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]]

