作者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