pdf("rpd15_1ar.pdf") grimes <- matrix(c( 1, 0.45, 0.34170, 2, 0.45, -0.00438, 3, 0.45, 0.82531, 4, 1.30, 1.77967, 5, 1.30, 0.95384, 6, 1.30, 0.64080, 7, 2.40, 1.75136, 8, 2.40, 1.27497, 9, 2.40, 1.17332, 10, 4.00, 3.12273, 11, 4.00, 2.60958, 12, 4.00, 2.57429, 13, 6.10, 3.17881, 14, 6.10, 3.00782, 15, 6.10, 2.67061, 16, 8.05, 3.05959, 17, 8.05, 3.94321, 18, 8.05, 3.43726, 19, 11.15, 4.80735, 20, 11.15, 3.35583, 21, 11.15, 2.78309, 22, 13.15, 5.13825, 23, 13.15, 4.70274, 24, 13.15, 4.25702, 25, 15.00, 3.60404, 26, 15.00, 4.15029, 27, 15.00, 3.42484),byrow=T,ncol=3) t <- as.matrix(grimes[,2]) y <- as.matrix(grimes[,3]) theta <- matrix(c(5,2,1),ncol=1) n <- nrow(y) p <- nrow(theta) difftheta <- 100 while(difftheta > .000001) { yhat <- theta[1,] * (1-exp(-(t/theta[2,])^theta[3,])) F1 <- 1-exp(-(t/theta[2,])^theta[3,]) F2 <- -(theta[1,]*theta[3,]/theta[2,]) * ((t/theta[2,])^theta[3,])*exp(-(t/theta[2,])^theta[3,]) F3 <- theta[1,] * log(t/theta[2,])*((t/theta[2,])^theta[3,])*exp(-(t/theta[2,])^theta[3,]) F <- cbind(F1,F2,F3) theta1 <- theta + solve(t(F) %*% F) %*% t(F) %*% (y-yhat) difftheta <- t(theta1-theta) %*% (theta1-theta) theta <- theta1 } yhat <- theta[1,] * (1-exp(-(t/theta[2,])^theta[3,])) F1 <- 1-exp(-(t/theta[2,])^theta[3,]) F2 <- -(theta[1,]*theta[3,]/theta[2,]) * ((t/theta[2,])^theta[3,])*exp(-(t/theta[2,])^theta[3,]) F3 <- theta[1,] * log(t/theta[2,])*((t/theta[2,])^theta[3,])*exp(-(t/theta[2,])^theta[3,]) F <- cbind(F1,F2,F3) s2 <- (t(y-yhat) %*% (y-yhat))/(n-p) vtheta <- s2[1]*solve(t(F) %*% F) setheta <- sqrt(diag(vtheta)) ttheta <- theta/setheta ptheta <- 2*(1-pt(abs(ttheta),n-p)) thetalo <- theta-setheta*qt(.975,n-p) thetahi <- theta+setheta*qt(.975,n-p) print(cbind(theta, setheta, ttheta, ptheta, thetalo, thetahi)) ybar <- matrix(rep(1/n,n*n),ncol=n) %*% y sstotalu <- t(y) %*% y sstotalc <- t(y-ybar) %*% (y-ybar) ssresidual <- t(y-yhat) %*% (y-yhat) ssmodel <- sstotalu-ssresidual dftotalu <- n dftotalc <- n-1 dfresidual <- n-p dfmodel <- p msmodel <- ssmodel/dfmodel msresidual <- ssresidual/dfresidual print(cbind(dfmodel, ssmodel, msmodel)) print(cbind(dfresidual, ssresidual, msresidual)) print(cbind(dftotalu, sstotalu)) print(cbind(dftotalc, sstotalc)) # Proportional Response estimates ht <- matrix(c(1,5,15),ncol=1) httheta <- 1-exp(-(ht/theta[2])^theta[3]) dhb1 <- 0*ht dhb2 <- -(theta[3]/theta[2])*((ht/theta[2])^theta[3])*exp(-(ht/theta[2])^theta[3]) dhb3 <- (log(ht/theta[2]))*((ht/theta[2])^theta[3])*exp(-(ht/theta[2])^theta[3]) dhb <- cbind(dhb1,dhb2,dhb3) s2h <- dhb %*% vtheta %*% t(dhb) seh <- sqrt(diag(s2h)) seh <- as.matrix(seh) ciht <- qt(.975,n-p)*seh ciht htthetalo <- httheta - ciht htthetahi <- httheta + ciht print(cbind(ht,httheta,htthetalo,htthetahi)) dev.off()