Hey Sky
2010-Sep-15 19:37 UTC
[R] optim with BFGS--what may lead to this, a strange thing happened
Dear R Users
on a self-written function for calculating maximum likelihood probability (plz
check function code at the bottom of this message), one value, wden, suddenly
jump to zero. detail info as following:
w[11]=2.14
lnw =2.37 2.90 3.76 ...
regw =1.96 1.77 1.82 ....
wden=0.182 0.178 0.179...
w[11]=2.14
lnw=2.37 2.90 3.76 ...
regw =1.96 1.77 1.82 ....
wden=0.182 0.179 0.180...
w[11]=-0.51
lnw=2.37 2.90 3.76 ...
regw=57.92 57.50 57.60 ...
wden=0 0 0 0 ....1 0 0 0
which I use cat() to get (by the way, how to write from the result from cat()
directly to a file?).
the lnw is coming from data. and w[11] is one of the estimated parameters; wden
and regw is estimated equation. wrk is a 0/1 dummy
wden equation: wden<-ifelse(wrk==1,dnorm((lnw-regw)/w[5])/w[5],1)
when I use simulated data, it works fine. but when using real data, it does not
work. when I did not add the following into the function,
func<-ifelse(func<.Machine$double.xmin,0.000001,func)
it will iterate for a while and give a non-finite finite difference in par[54]
report. but when I add it, very soon it will give the upper result. what may
lead to this and how to fix it?
thanks for any suggestion in advance.
Nan
from Montreal
# the main function
mymatrix<-function(par,data) {
# define the parameter matrix used in following part
vbar2<-matrix(0,n,nt)
vbar3<-matrix(0,n,nt)
v8 <-matrix(0,n,nt)
regw<-matrix(0,n,nt)
wden<-matrix(0,n,nt)
lia<-matrix(0,n,nt)
ccl<-matrix(1,n,ns)
eta<-c(0,0)
# setup the parts for loglikelihood
q1<-exp(par[1])
pr1<-q1/(1+q1)
pr2<-1-pr1
eta[2]<-par[2]
a<-par[3:6]
b<-par[7:11]
w<-par[12:npar]
for(m in 1:ns){
regw<-exp(w[1]+w[2]*eta[m])*actr+exp(w[3]+ w[4]*eta[m])*acwrk
vbar2=a[1]+ eta[m]+regw*a[2]+acwrk*a[3]+actr*a[4]
vbar3=b[1]+b[2]*eta[m]+regw*b[3]+acwrk*b[4]+actr*b[5]
v8=1+exp(vbar2)+exp(vbar3)
lia<-ifelse(home==1,1/v8,
ifelse(wrk==1,exp(vbar2)/v8,
ifelse(tr==1,exp(vbar3)/v8,1)))
wden<-ifelse(wrk==1,dnorm((lnw-regw)/w[5])/w[5],1)
ccl[,m]<-lia[,1]*lia[,2]*lia[,3]*lia[,4]*lia[,5]*lia[,6]*lia[,7]*lia[,8]*
wden[,1]*wden[,2]*wden[,3]*wden[,4]*wden[,5]*wden[,6]*wden[,7]*wden[,8]
}
func<-pr1*ccl[,1]+pr2*ccl[,2]
#
func<-ifelse(func<.Machine$double.xmin,0.000001,func)
#
f<-sum(log(func))
return(-f)
}
#*******************************************
mydata<-read.table("C:/Documents and Settings/lciqss/Desktop/R-2.11.1/my
data_wage_eq.txt", head=F)
nt<<-8 # number of periods
ns<<-2
n<<-50 # number of people
npar<<-16 # number of parameters
tr<-as.matrix(mydata[,1:nt])
wrk<-as.matrix(mydata[,(nt+1):(2*nt)])
home<-as.matrix(mydata[,(2*nt+1):(3*nt)])
actr<-as.matrix(mydata[,(3*nt+1):(4*nt)])
acwrk<-as.matrix(mydata[,(4*nt+1):(5*nt)])
lnw<-as.numeric(mydata[,(5*nt+1)])
guess<-rep(0,times=npar)
guess[npar]<-1.0
system.time(r1<-optim(guess,mymatrix,data=mydata, hessian=F))
guess<-r1$par
system.time(rmat<-optim(guess,mymatrix,data=mydata,
method="BFGS",hessian=T,
control=list(trace=T, maxit=1000)))
rmat$par
std.err<-sqrt(diag(solve(rmat$hessian)))
result<-cbind(rmat$par,std.err,rmat$par/std.err)
colnames(result)<-c("paramters","std,err","t
test")
#print(result)
Hey Sky
2010-Sep-15 19:53 UTC
[R] optim with BFGS--what may lead to this, a strange thing happened
sorry, a typo. the following code is an example from my computer and the real
one just has more variables.
the w[11] is the w[5] in the wden equation. the value of wden is calculated
based on the ifelse condition. it may equal to 1 when one is working; if not, it
will equal to the value from -- dnorm((lnw-regw)/w[5])/w[5]
----- Original Message ----
From: Hey Sky <heyskywalker at yahoo.com>
To: R <r-help at r-project.org>
Sent: Wed, September 15, 2010 3:37:23 PM
Subject: optim with BFGS--what may lead to this, a strange thing happened
Dear R Users
on a self-written function for calculating maximum likelihood probability (plz
check function code at the bottom of this message), one value, wden, suddenly
jump to zero. detail info as following:
w[11]=2.14
lnw =2.37 2.90 3.76 ...
regw =1.96 1.77 1.82 ....
wden=0.182 0.178 0.179...
w[11]=2.14
lnw=2.37 2.90 3.76 ...
regw =1.96 1.77 1.82 ....
wden=0.182 0.179 0.180...
w[11]=-0.51
lnw=2.37 2.90 3.76 ...
regw=57.92 57.50 57.60 ...
wden=0 0 0 0 ....1 0 0 0
which I use cat() to get (by the way, how to write from the result from cat()
directly to a file?).
the lnw is coming from data. and w[11] is one of the estimated parameters; wden
and regw is estimated equation. wrk is a 0/1 dummy
wden equation: wden<-ifelse(wrk==1,dnorm((lnw-regw)/w[5])/w[5],1)
when I use simulated data, it works fine. but when using real data, it does not
work. when I did not add the following into the function,
func<-ifelse(func<.Machine$double.xmin,0.000001,func)
it will iterate for a while and give a non-finite finite difference in par[54]
report. but when I add it, very soon it will give the upper result. what may
lead to this and how to fix it?
thanks for any suggestion in advance.
Nan
from Montreal
# the main function
mymatrix<-function(par,data) {
# define the parameter matrix used in following part
vbar2<-matrix(0,n,nt)
vbar3<-matrix(0,n,nt)
v8 <-matrix(0,n,nt)
regw<-matrix(0,n,nt)
wden<-matrix(0,n,nt)
lia<-matrix(0,n,nt)
ccl<-matrix(1,n,ns)
eta<-c(0,0)
# setup the parts for loglikelihood
q1<-exp(par[1])
pr1<-q1/(1+q1)
pr2<-1-pr1
eta[2]<-par[2]
a<-par[3:6]
b<-par[7:11]
w<-par[12:npar]
for(m in 1:ns){
regw<-exp(w[1]+w[2]*eta[m])*actr+exp(w[3]+ w[4]*eta[m])*acwrk
vbar2=a[1]+ eta[m]+regw*a[2]+acwrk*a[3]+actr*a[4]
vbar3=b[1]+b[2]*eta[m]+regw*b[3]+acwrk*b[4]+actr*b[5]
v8=1+exp(vbar2)+exp(vbar3)
lia<-ifelse(home==1,1/v8,
ifelse(wrk==1,exp(vbar2)/v8,
ifelse(tr==1,exp(vbar3)/v8,1)))
wden<-ifelse(wrk==1,dnorm((lnw-regw)/w[5])/w[5],1)
ccl[,m]<-lia[,1]*lia[,2]*lia[,3]*lia[,4]*lia[,5]*lia[,6]*lia[,7]*lia[,8]*
wden[,1]*wden[,2]*wden[,3]*wden[,4]*wden[,5]*wden[,6]*wden[,7]*wden[,8]
}
func<-pr1*ccl[,1]+pr2*ccl[,2]
#
func<-ifelse(func<.Machine$double.xmin,0.000001,func)
#
f<-sum(log(func))
return(-f)
}
#*******************************************
mydata<-read.table("C:/Documents and Settings/lciqss/Desktop/R-2.11.1/my
data_wage_eq.txt", head=F)
nt<<-8 # number of periods
ns<<-2
n<<-50 # number of people
npar<<-16 # number of parameters
tr<-as.matrix(mydata[,1:nt])
wrk<-as.matrix(mydata[,(nt+1):(2*nt)])
home<-as.matrix(mydata[,(2*nt+1):(3*nt)])
actr<-as.matrix(mydata[,(3*nt+1):(4*nt)])
acwrk<-as.matrix(mydata[,(4*nt+1):(5*nt)])
lnw<-as.numeric(mydata[,(5*nt+1)])
guess<-rep(0,times=npar)
guess[npar]<-1.0
system.time(r1<-optim(guess,mymatrix,data=mydata, hessian=F))
guess<-r1$par
system.time(rmat<-optim(guess,mymatrix,data=mydata,
method="BFGS",hessian=T,
control=list(trace=T, maxit=1000)))
rmat$par
std.err<-sqrt(diag(solve(rmat$hessian)))
result<-cbind(rmat$par,std.err,rmat$par/std.err)
colnames(result)<-c("paramters","std,err","t
test")
#print(result)