作者tokyo291 (工口工口)
看板Statistics
标题[程式] 用数值分析方法找方程式的根
时间Fri Jun 7 01:02:25 2013
[软体程式类别]:
R
[程式问题]:
想算zero-inflated poisson 的MLE
以解出两个变数的理论MLE值
想将样本的资料代入去求MLE,利用数值分析的方法去求MLE方程式的解
发现解的数值很奇怪
[软体熟悉度]:
高(1年以上)
[问题叙述]:
iid
X_i----> f(x_i)= (1-π)+π*exp(-θ) ,x_i=0
exp(-θ)*θ^(x_i)
π*----------------- ,x_i=1,2,3,...
x_i!
for i=1,...,n
X_i服从ZIP(π,θ),π和θ的MLE分别为
_ θhat n-n_0
x=------------- ; πhat=------------------
1-exp(-θhat) n*(1-exp(-θhat))
其中n_0:x_i for i=1,...,n 中0的个数
_
x :(不为0的x_i的总合)/(n-n_0)
_
所以我就先生成一笔ZIP的data去得到x 然後先去解
_ θhat
x=------------- 这个方程式得到θhat,再去求πhat
1-exp(-θhat)
我利用数值分析里面的Newton-Raphson Method去解方程式的根
得到的答案很奇怪
[程式范例]:
y<-NULL
n<-2000
p<-c(0.3,2)
从ZIP(π=0.3,θ=2生成)
for(i in 1:n){
y[i]<-rbinom(1,1,p[1])
if(y[i]==0){y[i<-y[i]]}
else{y[i]<-rpois(1,p[2])}
}
table(y)
num.0<-as.vector(as.matrix(table(y))[1,])
_
mean.num0<-sum(y)/num.0
此为x
eq<-function(x){x/(1-exp(-x))-mean.num0}
θhat的方程式
eq1<-function(x){(1-exp(-x)-x*exp(-x))/((1-exp(-x))^2)}
上式的一阶微分
int<-3
起始值设定
xt<-NULL;xt[1]<-int
for(i in 1:100){
if(i==1){
xt[i+1]<-xt[i]-eq(xt[i])/eq1(xt[i])
}
else{
if(abs(xt[i]-xt[i-1])>10^(-6)){
xt[i+1]<-xt[i]-eq(xt[i])/eq1(xt[i])
}
else{break}
}
}
因为使用误差小於10^(-6)作为停止条件
为了避免xt[1]-xt[0]的情况所以多了一个if回圈
得到的答案是-0.5612140
感觉很奇怪
请问是哪里出了问题呢?
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 218.164.76.23
※ 编辑: tokyo291 来自: 218.164.76.23 (06/07 01:02)
※ 编辑: tokyo291 来自: 218.164.76.23 (06/07 01:06)
1F：→ Wush978:我用curve(eq1, from=-100, to=100) 06/07 22:37
2F：→ Wush978:看起来没有根耶 06/07 22:37
3F：→ Wush978:有兴趣欢迎到R_Language板讨论 06/07 22:37