id.prob<-function (tt) { library(mvtnorm) #============================ Makeham<-function(tt) { a2=0.030386513 a3=0.006688287 b3=0.039047537 t<-tt-20 h.t<-a2+a3*exp(b3*t) S.t<-exp(-a2*t+a3/b3*(1-exp(b3*t))) return(S.t*h.t) } #=========================== trans<-function (age) { indic=c(2,2) lage=log(age) xpars=c(2.1862908,7.5528077,8.5806697,2.3319461,8.8959507,9.1380187,0.3720293) beta<-c(0,0) alpha1<-c(0,0) alpha2<-c(0,0) beta[1]<-xpars[1] alpha1[1]<-xpars[2] alpha2[1]<-xpars[3] beta[2]<-xpars[4] alpha1[2]<-xpars[5] alpha2[2]<-xpars[6] corr<-matrix(rep(0,4),nc=2) corr[1,1]=1 corr[2,2]=1 corr[1,2]=xpars[7] corr[2,1]=corr[1,2]
LLK<-0 L<-c(0,0) R<-c(0,0)
for(j in 1:2){ if(indic[j]==1){ L[j]<--Inf R[j]<-alpha1[j]-beta[j]*lage} if(indic[j]==2){ L[j]<-alpha1[j]-beta[j]*lage R[j]<-alpha2[j]-beta[j]*lage} if(indic[j]==3){ L[j]<-alpha2[j]-beta[j]*lage R[j]<-Inf} }
LLK<-pmvnorm(lower=L,upper=R,corr=corr)[1]
return(LLK) }
#===========================
tot<-Makeham(tt)*trans(tt) return(tot) }
And then:
> id.prob(20) [1] 0.0001417763 > id.prob(21) [1] 0.00018078 > id.prob(20:21) [1] 0.0001417763 0.0001375794 Warning messages: 1: number of items to replace is not a multiple of replacement length 2: number of items to replace is not a multiple of replacement length 3: number of items to replace is not a multiple of replacement length 4: number of items to replace is not a multiple of replacement length
Ultimately I need to be able to integrate this function, which isn't going to happen (correctly) if I've messed-up writing the function. Thanks for any advice on this.
Regards, Lyle W. Konigsberg Professor of Anthropology University of Tennessee Knoxville, TN [EMAIL PROTECTED]
______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html