【R语言-20行代码】牛顿迭代法求伽马函数极大似然估计法的参数估计

深藏阁楼爱情的钟 2022-04-21 22:30 374阅读 0赞

## 简述 ##

研究了下计算公式，简化了一下，用r语言实现了。

## 算法解释 ##

*  牛顿迭代法  
     x k + 1 = x k − f ( x k ) f ′ ( x k ) x\_\{k+1\} = x\_k - \\frac\{f(x\_k)\}\{f&\#x27;(x\_k)\} xk\+1=xk−f′(xk)f(xk)  
    求解的方程是  
     f ( x ) = 0 f(x) = 0 f(x)=0
 *  通过极大似然估计，构造对数似然方程，之后再关于 α \\alpha α和 β \\beta β 求偏导数。之后，得到关于 α \\alpha α的非线性对数似然方程。然后， β \\beta β可以用 α \\alpha α表示。
 *  再进一步的简化（去掉无关的项，再整理相关的项）  
    得到需要求解的方程为  
     l o g ( α ) + d i g a m m a ( α ) = 0 log(\\alpha) + digamma(\\alpha) = 0 log(α)\+digamma(α)=0
 *  再求一下这个方程左边的导数  
     1 α + t r i g a m m a ( α ) \\frac\{1\}\{\\alpha\} + trigamma(\\alpha) α1\+trigamma(α)
 *  初始值，使用通过矩估计得到的参数。

## 代码部分 ##

除掉前面的读取数据加一行的空格，不就是小于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)