# Functions for P matrices in Gompertz and Weibull models # Ardo, UCL 2016 #################################### # Gompertz model: # #################################### Pmatrix.Gomp<-function(t1,t2,lambda,xi,nnodes,h){ # xi cannot be zero: if( any(abs(xi)1){P[1,2]<-0}} #P[1,1:2]<-P[1,1:2]/som} # Rest of P-matrix: P[1,3]<-1-P[1,1]-P[1,2] P[2,2]<- exp(-lambda[3]/xi[3]*(exp(xi[3]*t2)-exp(xi[3]*t1))) P[2,3]<-1-P[2,2] P } #################################### # Weibull model: # #################################### Pmatrix.Weib<-function(t1,t2,lambda,tau,nnodes,h){ # Construct P-matrix: P<-diag(3) P[1,1]<-exp(-lambda[1]*(t2^tau[1]-t1^tau[1])-lambda[2]*(t2^tau[2]-t1^tau[2])) # Integrand for P[1,2]: integrand<-function(u){ p11<-exp(-lambda[1]*(u^tau[1]-t1^tau[1])-lambda[2]*(u^tau[2]-t1^tau[2])) q12<-lambda[1]*tau[1]*u^(tau[1]-1) p22<-exp(-lambda[3]*(t2^tau[3]-u^tau[3])) p11*q12*p22 } # Approx using composite Simpson's rule (code from Wikipedia): S <- integrand(t1) for(i in seq(1,nnodes,by=2)){ x <- t1+ h*i; S = S+ 4 * integrand(x)} for(i in seq(2,nnodes-1,by=2)){x <- t1 + h*i; S = S+ 2 * integrand(x)} S<-S+integrand(t2) P[1,2]<- h * S / 3 # Check approximation: #if(is.na(P[1,2])){P[1,2]<-0}else{if((P[1,1]+P[1,2])>1){P[1,2]<-0}} #P[1,1:2]<-P[1,1:2]/som} # Rest of P-matrix: P[1,3]<-1-P[1,1]-P[1,2] P[2,2]<- exp(-lambda[3]*(t2^tau[3]-t1^tau[3])) P[2,3]<-1-P[2,2] P }