【R語言-20行程式碼】牛頓迭代法求伽馬函式極大似然估計法的引數估計
阿新 • • 發佈:2018-11-03
簡述
研究了下計算公式,簡化了一下,用r語言實現了。
演算法解釋
-
牛頓迭代法
求解的方程是
-
通過極大似然估計,構造對數似然方程,之後再關於 和 求偏導數。之後,得到關於 的非線性對數似然方程。然後, 可以用 表示。
-
再進一步的簡化(去掉無關的項,再整理相關的項)
得到需要求解的方程為
-
再求一下這個方程左邊的導數
-
初始值,使用通過矩估計得到的引數。
程式碼部分
除掉前面的讀取資料加一行的空格,不就是小於20行咩
- TIMES 的是迭代次數
- 1e-5 表示的是最小的變動精度
講真,我用這個精度,我算了4次迭代,就得到正確結果了。
library(xlsx)
ray = read.xlsx('D:/Code/R/Data in Excel/Chapter 8/gamma-arrivals.xls',1)
mean_ray = mean(ray[,1])
var_ray = var(ray[,1])
alpha = mean_ray**2 / var_ray
origin_X = alpha
f = function(a){
log(a)+ digamma(a)
}
ff = function(a){
1. / a + trigamma(a)
}
TIMES = 10
for (i in 1:TIMES){
x = origin_X -f(origin_X) / (ff(origin_X))
if (abs(x - origin_X) < 1e-5) {
print(i)
break
} else {
origin_X = x
}
}
print(x)