簡述

研究了下計算公式，簡化了一下，用r語言實現(xiàn)了。

算法解釋

• 牛頓迭代法
  $$x_{k+1}=x_k?\frac{f(x_k)}{f^{\prime }(x_k)}$$
  求解的方程是
  $$f(x)=0$$

• 通過極大似然估計，構(gòu)造對數(shù)似然方程，之后再關(guān)于$\alpha$和$\beta$ 求偏導(dǎo)數(shù)。之后，得到關(guān)于$\alpha$的非線性對數(shù)似然方程。然后，$\beta$可以用$\alpha$表示。

• 再進一步的簡化（去掉無關(guān)的項，再整理相關(guān)的項）
  得到需要求解的方程為
  $$log(\alpha )+digamma(\alpha )=0$$

• 再求一下這個方程左邊的導(dǎo)數(shù)
  $$\frac{1}{\alpha }+trigamma(\alpha )$$

• 初始值，使用通過矩估計得到的參數(shù)。

代碼部分

除掉前面的讀取數(shù)據(jù)加一行的空格，不就是小于20行咩

• TIMES 的是迭代次數(shù)
• 1e-5 表示的是最小的變動精度

講真，我用這個精度，我算了4次迭代，就得到正確結(jié)果了。

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)

總結(jié)

以上是生活随笔為你收集整理的【R语言-20行代码】牛顿迭代法求伽马函数极大似然估计法的参数估计的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。