I am trying to code a steepest ascent algorithm to optimize parameters used
in a survivor function type problem. My unknown parameters (alpha, Beta0,
and Beta1) for which I have been able to optimize using Newton's method. I
keep getting an error because my alpha becomes negative and I can't
calculate the likelihood.
Here is my log likelihood I am optimizing (in LaTex):
l=\sum _{ i=1 }^{ n }{ \left( { w }_{ i }log\left\{ { t }_{ i }^{ \alpha
}exp\left\{ { \beta }_{ o }+{ \delta }_{ i }{ \beta }_{ 1 } \right\}
\right\} \right) - } { t }_{ i }^{ \alpha }exp\left\{ { \beta }_{ o }+{
\delta }_{ i }{ \beta }_{ 1 } \right\} +\sum _{ i=1 }^{ n }{ { w }_{ i
}log\left\{ \alpha { t }^{ -1 } \right\} } \\
Here is the code I have so far towards the bottom is the data I have been
given a text version of the data:
#Length of remission
lordata=read.table('~/Dropbox/CSUF/M534/Rick-Math 534/Midterm
1/Table2_2GH.txt',header=T)
attach(lordata)
#This is my loglikelihood function
#w = 0 or 1 depending on whether they were censored or not
#t is the suvival time
#d is 0 if they were the control group and 1 if they were in the treatment
group
#theta is a vector of paramters (a,Bo, and B1 over which to maximize upon
lordata[4,3]
like<-function(theta,data,gcomp=FALSE,hesscomp=FALSE)
{
a = theta[1]
Bo=theta[2]
B1=theta[3]
d=data[,1] #treatment
w=data[,2] #censored
t=data[,3] #survival time
mu = (t^a)*exp(Bo+d*B1)
l = sum(w*log(mu)-mu+w*log(a/t))
if(gcomp==TRUE) {
grad=matrix(0,3,1)
grad[1]=sum(w*log(t)-mu*log(t)+w/a) #alpha
grad[2]=sum(w-mu) #Bo
grad[3]=sum(d*w-d*mu) #B
}
if(hesscomp==TRUE){
#Recall mu = (t^a)*exp(Bo+d*B1)
hess = matrix(0,3,3)
hess[1,1] = sum(-(mu*log(t)^2)-w/(a^2)) #dada
hess[1,2] = hess[2,1] = sum(-log(t)*mu) #dadB0
hess[2,2] = sum(-mu) #dBodBo
hess[1,3] = hess[3,1] = sum(-d*log(t)*mu) #dadB1
hess[2,3] = hess[3,2] = sum(-d*mu) #dBodB1
hess[3,3] = sum(-mu*d^2) #dB1dB1
}
list(l=l,grad=if(gcomp) grad,hess=if(hesscomp)hess)
}
#Steepest ascent algorithm----------------------------------
ascent2 <- function(data,maxit,theta) {
for (it in 1:maxit){
a =like(theta,data,gcomp=T,hesscomp=F)
theta1 = theta + a$grad # Steepest Ascent
atmp = like(theta,data,gcomp=F,hesscomp=F)
halve = 0;
print(c(it, halve, a$l))
print(c(it, halve, atmp$l))
while (theta[1]<=0)
{
halve = halve + 1
theta1 = theta + a$grad/2^(halve)
}
atmp = like(theta1,data,gcomp = FALSE)
while (atmp$l < a$l & halve <= 20){
halve = halve+1
theta1 = theta + a$grad/2^(halve) # Steepest Ascent
atmp = like(theta1,data,gcomp=F,hesscomp=F)
print(c(it, halve,atmp$l))
}
if (halve >= 20) print('Step-halving failed after 20 halvings')
theta = theta1
print(sprintf('it = %2.0f a = %12.12f B0 = %12.12f B1 = %12.12f
l=%12.12f',
it,theta[1],theta[2],theta[3],atmp$l))
print ('-----------------------------------------')
}
}
#Data in Text form
Treatment Censored Time
1 0 6
1 1 6
1 1 6
1 1 6
1 1 7
1 0 9
1 0 10
1 1 10
1 0 11
1 1 13
1 1 16
1 0 17
1 0 19
1 0 20
1 1 22
1 1 23
1 0 25
1 0 32
1 0 32
1 0 34
1 0 35
0 1 1
0 1 1
0 1 2
0 1 2
0 1 3
0 1 4
0 1 4
0 1 5
0 1 5
0 1 8
0 1 8
0 1 8
0 1 8
0 1 11
0 1 11
0 1 12
0 1 12
0 1 15
0 1 17
0 1 22
0 1 23
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
