Hello,
The likelihood includes two parameters to be estimated: lambda
(=beta0+beta1*x) and alpha. The algorithm for the estimation is as
following:
1) with alpha=0, estimate lambda (estimate beta0 and beta1 via GLM)
2) with lambda, estimate alpha via ML estimation
3) with updataed alpha, replicate 1) and 2) until alpha is converged to a
value
I coded 1) and 2) (it works), but faced some problems with the loop. It
produced some errors. Is there someone to help me to figure this problem for
the loop. Here are two codes: [I] one is for 1) and 2) (working), and [II]
another one is for 1)-3) (making errors).
Thank you in advance,
Jin
[I]for 1) and 2) (working)
% data set
alpha<-1
verpi<-c(5^alpha,10^alpha-5^alpha,14^alpha-10^alpha,18^alpha-14^alpha)
r<-c(1,0,0,1)
k<-c(3,2,2,2)
x<-c(0.5,0.5,1.0,1.0)
% estimate lambda (lambda=beta0+beta1*x)
GLM_results <- glm(r/k ~ x, family=binomial(link='cloglog'),
offset=log(verpi),weights=k)
beta0<-GLM_results$coefficients[[1]]
beta1]<-GLM_results$coefficients[[2]]
lambda1<-beta0+beta1*x[1]
lambda2<-beta0+beta1*x[2]
% using lambda, estimate alpha (a=alpha) through ML estimation
L<-function(a){
s1_f1<-(exp(-lambda1*(0^a-0^a))-exp(-lambda1*(5^a-0^a)))
s2_f2<-(exp(-lambda1*(10^a)-lambda2*(14^a-10^a+14^a-14^a))
-exp(-lambda1*(10^a)-lambda2*(14^a-10^a+18^a-14^a)))
logl<-log(s1_f1*s2_f2)
return(-logl)}
optim(1,L)
alpha<-optim(1,L)$par
}
[II] for 1)-3) (making errors)
alpha[1]<-1
beta0=rep(0,10)
beta1=rep(0,10)
lambda1=rep(0,10)
lambda2=rep(0,10)
lambda3=rep(0,10)
lambda4=rep(0,10)
verpi=rep(0,10)
s1_f1=rep(0,10)
s2_f2=rep(0,10)
a=rep(0,10)
logl=rep(0,10)
beta0[1]=0
beta1[1]=0
lambda1[1]=0
lambda2[1]=0
lambda3[1]=0
lambda4[1]=0
verpi[1]=0
s1_f1[1]=0
s2_f2[1]=0
a[1]=1
logl[1]=0
for(i in 2:11){
verpi[i]<-c(5^alpha[i-1],10^alpha[i-1]-5^alpha[i-1],14^alpha[i-1]-10^alpha[i-1],18^alpha[i-1]-14^alpha[i-1])
r<-c(1,0,0,1)
k<-c(3,2,2,2)
x<-c(0.5,0.5,1.0,1.0)
GLM_results <- glm(r/k ~ x, family=binomial(link='cloglog'),
offset=log(verpi[i]),weights=k)
GLM_results
logLik(GLM_results[i])
beta0[i]<-GLM_results$coefficients[[1]]
beta1[i]<-GLM_results$coefficients[[2]]
beta0[i]
beta1[i]
lambda1[i]<-beta0[i]+beta1[i]*x[1]
lambda2[i]<-beta0[i]+beta1[i]*x[2]
lambda3[i]<-beta0[i]+beta1[i]*x[3]
lambda4[i]<-beta0[i]+beta1[i]*x[4]
lambda1[i]
lambda2[i]
lambda3[i]
lambda4[i]
L<-function(a){
s1_f1[i]<-(exp(-lambda1[i]*(0^a[i]-0^a[i]))-exp(-lambda1[i]*(5^a[i]-0^a[i])))
s2_f2[i]<-(exp(-lambda1[i]*(10^a[i])-lambda2[i]*(14^a[i]-10^a[i]+14^a[i]-14^a[i]))
-exp(-lambda1*(10^a[i])-lambda2[i]*(14^a[i]-10^a[i]+18^a[i]-14^a[i])))
logl[i]<-log(s1_f1[i]*s2_f2[i])
return(-logl[i])}
optim(1,L)
alpha[i]<-optim(1,L)$par
}
alpha