基于R和Python的极大似然估计的牛顿法实现

最新推荐文章于 2023-09-13 13:51:09 发布

张乃晟

最新推荐文章于 2023-09-13 13:51:09 发布

阅读量6.9k

点赞数 12

文章标签： python 概率论算法

本文链接：https://blog.csdn.net/zns972630879/article/details/120399944

版权

0.前言

最近在学习Theory and Method of Statistics（统计理论方法）,使用的教材是由Bradley Efron 、Trevor Hastie共同编写的Computer Age Statistical Inference: Algorithms, Evidence, and Data Science(《计算机时代的统计推断：算法、演化和数据科学》)。书中第四章讲述的Fisherian Inference and Maximum Likelihood Estimation（费雪推断和极大似然估计）,其中提到现实应用中极大似然估计并没有那么容易求解,比如Cauchy分布和Gamma分布。

如果极大似然估计方法没有显式解,可以考虑用数值计算的方法求解(如牛顿法);更进一步，如果二阶导不存在或Hessian矩阵非正定,可以使用拟牛顿法;再复杂一些,可以使用MM算法(EM是MM的特例) 。本文以牛顿法为例,给出求解 Cauchy分布、Gamma分布的极大似然估计参数的理论并使用R和Python实现。

1.理论基础

本节给出牛顿法求分布的极大似然参数估计的一般理论。

如果随机变量 $x$ 独立同分布于 $f(x|\theta)$ ,且已知一组样本 $X=\{x_1,...,x_n\}$ ,为了估计该分布的参数，可以使用极大似然估计的方法。

首先写出样本的似然函数

$L(\theta)=f(x_1,...,x_n|\theta)=\prod_{i=1}^{n}f(x_i|\theta)$

对 $L(\theta)$ 进行对数化处理,得到对数似然函数 $l(\theta)$

$l(\theta)=ln(L(\theta))=\sum_{i=1}^{n}ln(f(x_i,\theta))$

则求解未知参数 $\theta$ 等价于求解以下等式方程组

$\dot{l}(\hat{\theta})=\frac{\partial{l(\hat{\theta})}}{\partial{\hat{\theta}}}=\vec{0}$

不妨假设收敛解为 $\hat{\theta}$ ,将 $\dot{l}(\theta)$ 在 $\hat{\theta}$ 的邻域内展开成泰勒级数得

$\dot{l}(\theta)=\dot{l}(\hat{\theta})+\ddot{l}(\hat{\theta})(\theta-\hat{\theta})+o(||\theta-\hat{\theta}||_2^2)$

$\dot{l}(\theta)=\dot{l}(\hat{\theta})+\ddot{l}(\hat{\theta})(\theta-\hat{\theta}),\theta\rightarrow \hat{\theta}$

$\dot{l}(\hat{\theta})=0\Rightarrow \theta=\hat{\theta}-(\ddot{l}(\hat{\theta}))^-\dot{l}(\hat{\theta})$

这样就得到一个迭代关系式

$\hat{\theta}^{(k+1)}=\hat{\theta}^{(k)}-(\ddot{l}(\hat{\theta}^{(k)}))^-\dot{l}(\hat{\theta}^{(k)})$

如果是连续的，并且待求的零点是孤立的，那么在零点周围存在一个区域，只要初始值位于这个邻近区域内，那么牛顿法必定收敛。并且，如果不为0, 那么牛顿法将具有平方收敛的性能。粗略地说，这意味着每迭代一次，牛顿法结果的有效数字将增加一倍。

2.Cauchy分布的极大似然估计

如果随机变量 $x$ 服从柯西分布,记为 $x\sim{C(\gamma,\theta)}$ ,其中 $\gamma$ 为最大值一半处的一半宽度的尺度参数(scale parameter ), $\theta$ 为定义分布峰值位置的位置参数(location parameter )

$f(x|\gamma,\theta)=\frac{1}{\pi\gamma}\frac{1}{(1+(\frac{x-\theta}{\gamma})^2)}=\frac{1}{\pi}\frac{\gamma}{({x-\theta})^2+\gamma^2},x\in(-\infty,\infty)$