songlin
这是关于一个最小二乘求参数的问题,我自己编的程序求出来的结果跟真值差得很远,但是又觉得思路没错,麻烦高手帮我看看程序,看看怎么改进才能得到好结果:
模型是:y=beta*X+g(t)+e
g(t)由于是无限维的,所以就用
Gm<-function(t,b0,b1,b2)
+ {((b1-b0)/(0.5)*t+b0*0.5/0.5)*I1(t)+((b2-b1)/0.5*t-(b2*0.5-b1*1)/0.5)*I2(t)}
替代。
然后就是由县委的问题了,参数是beta,b0,b1,b2
求的是目标函数最小,求一阶偏导后得到 f 函数,二阶偏导后得到df,然后就是利用牛顿迭代法做的
真值beta=3,b0=0,b1=0.5,b2=2
下边是我自己编的程序。
希望高手帮我看看,万分感谢
newt.result<-function(x){
t<-runif(500,0,1)
X<-runif(500,0,2)
Fe_new<-runif(500,0,1)
Fe<-function(Fe_new){for(i in 1:500){if(Fe_new<=0.01) y<-0.01
else if(Fe_new>=0.99) y<-0.99
else y<-Fe_new
};
return(y)}
Fe<-Fe(Fe_new)
e<-log(Fe/(1-Fe))
Fe1<-runif(500,0,1)
e1<-log(Fe1/(1-Fe1))
Z1<-e1+11/3
Z<-function(Z1){for(i in 1:500){if(Z1<0) y<-0
else y<-Z1}
return(y)}
Z<-Z(Z1)
y<-runif(500,0,1)
d<-function(X,t,e,Z){for(i in 1:500){if(3*X[i]+2*t[i]*t[i]+e[i]<=Z[i]) y[i]<-1 else y[i]<-0 };return(y)}
delta<-d(X,t,e,Z)
I1<-function(t){for(i in 1:500){if(t[i]<0.5) y[i]<-1 else y[i]<-0 };return(y)}
I2<-function(t){for(i in 1:500){if(t[i]>=0.5) y[i]<-1 else y[i]<-0 };return(y)}
Gm<-function(t,b0,b1,b2)
{((b1-b0)/0.5*t+b0*0.5/0.5)*I1(t)+((b2-b1)/0.5*t-(b2*0.5-b1*1)/0.5)*I2(t)}
f<-function(delta,X,Z,t,beta,b0,b1,b2)
{f1<-1/500*2*(((-delta+1/(1+exp(-Z+X*beta+Gm(t,b0,b1,b2))))*(exp(-Z+X*beta+Gm(t,b0,b1,b2))/(1+exp(-Z+X*beta+Gm(t,b0,b1,b2)))^2))%*%X)
f2<-1/500*2*(((-delta+1/(1+exp(-Z+X*beta+Gm(t,b0,b1,b2))))*(exp(-Z+X*beta+Gm(t,b0,b1,b2))/(1+exp(-Z+X*beta+Gm(t,b0,b1,b2)))^2))%*%((-t/0.5+1)*I1(t)))
f3<-1/500*2*(((-delta+1/(1+exp(-Z+X*beta+Gm(t,b0,b1,b2))))*(exp(-Z+X*beta+Gm(t,b0,b1,b2))/(1+exp(-Z+X*beta+Gm(t,b0,b1,b2)))^2))%*%(t/0.5*I1(t)+(1-t)/0.5*I2(t)))
f4<-1/500*2*(((-delta+1/(1+exp(-Z+X*beta+Gm(t,b0,b1,b2))))*(exp(-Z+X*beta+Gm(t,b0,b1,b2))/(1+exp(-Z+X*beta+Gm(t,b0,b1,b2)))^2))%*%((t-0.5)/0.5*I2(t)))
return(c(f1,f2,f3,f4))}
df<-function(delta,X,Z,t,beta,b0,b1,b2)
{u<-1/500*(2*(-delta+1/(1+exp(-Z+X*beta+Gm(t,b0,b1,b2))))*((exp(-Z+X*beta+Gm(t,b0,b1,b2))*(1-exp(-Z+X*beta+Gm(t,b0,b1,b2))))/(1+exp(-Z+X*beta+Gm(t,b0,b1,b2)))^3)+2*(exp(-Z+X*beta+Gm(t,b0,b1,b2)))^2/((1+exp(-Z+X*beta+Gm(t,b0,b1,b2)))^2)^2)
G0d<-(-2*t+1)*I1(t)
G1d<-2*t*I1(t)+2*(1-t)*I2(t)
G2d<-2*(t-0.5)*I2(t)
df11<-sum(u*X*X)
df12<-sum(u*X*G0d)
df13<-sum(u*X*G1d)
df14<-sum(u*X*G2d)
df21<-sum(u*X*G0d)
df22<-sum(u*G0d*G0d)
df23<-sum(u*G0d*G1d)
df24<-sum(u*G0d*G2d)
df31<-sum(u*X*G1d)
df32<-sum(u*G0d*G1d)
df33<-sum(u*G1d*G1d)
df34<-sum(u*G1d*G2d)
df41<-sum(u*X*G2d)
df42<-sum(u*G0d*G2d)
df43<-sum(u*G1d*G2d)
df44<-sum(u*G2d*G2d)
return(matrix(c(df11,df12,df13,df14,df21,df22,df23,df24,df31,df32,df33,df34,df41,df42,df43,df44),nrow=4,byrow=T))}
newt.raph<-function(delta,X,Z,t,beta,b0,b1,b2,tol=0.0001,maxit=100,details=F)
{current.beta<-beta
current.b0<-b0
current.b1<-b1
current.b2<-b2
change<-1
iter<-0
while(change > tol &&iter <= maxit)
{
cov.temp<-df(delta,X,Z,t,current.beta,current.b0,current.b1,current.b2)
beta<-current.beta - (solve(cov.temp)%*%f(delta,X,Z,t,current.beta,current.b0,current.b1,current.b2))[1]
b0<-current.b0 -(solve(cov.temp)%*%f(delta,X,Z,t,current.beta,current.b0,current.b1,current.b2))[2]
b1<-current.b1- (solve(cov.temp)%*%f(delta,X,Z,t,current.beta,current.b0,current.b1,current.b2))[3]
b2<-current.b2- (solve(cov.temp)%*%f(delta,X,Z,t,current.beta,current.b0,current.b1,current.b2))[4]
change<-sum((c(beta-current.beta,b0-current.b0,b1-current.b1,b2-current.b2))^2)
current.beta<-beta
current.b0<-b0
current.b1<-b1
current.b2<-b2
iter<-iter + 1
}
return(c(current.beta,current.b0,current.b1,current.b2))}
return(newt.raph(delta,X,Z,t,2.5,0,0,1))}
newt.result(x)
sss<-function(x)
{x<-1
s2<-0
ss<-c(0,0,0,0)
iter<-0
while(iter<500)
{sx<-newt.result(x)
s2<-s2+sx[1]*sx[1]
ss<-ss+sx
x<-x+1
iter<-iter+1
}
return(c((ss/500)[1],s2/500-((ss/500)^2)[1]))}
sss(x)