开发: C++知识库 Java知识库 JavaScript Python PHP知识库人工智能区块链大数据移动开发嵌入式开发工具数据结构与算法开发测试游戏开发网络协议系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑笔记本显卡显示器固态硬盘硬盘耳机手机 iphone vivo oppo 小米华为单反装机图拉丁

-> 人工智能 -> 基于R和Python的极大似然估计的牛顿法实现 -> 正文阅读

[人工智能]基于R和Python的极大似然估计的牛顿法实现

0.前言

1.理论基础

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)$

?????? 当 $\gamma=1,\theta=0$ ,此时的Cauchy 分布称为标准Cauchy 分布

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

?????? Cauchy 分布的最特别的性质是其期望及高阶矩都不存在,自然也就无法对参数进行矩估计。但Cauchy分布的cdf具有很好的性质,可以利用一组样本的分位点来对参数进行点估计。

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

$\because F(\theta|\gamma,\theta)=\frac{1}{2}, F(\theta+\gamma|\gamma,\theta)=\frac{3}{4}$

$\therefore \hat{\theta}=X_{\frac{1}{2}},\hat{\theta}+\hat{\gamma}=X_{\frac{3}{4}}$

$\therefore \hat{\gamma}=X_{\frac{3}{4}}-X_{\frac{1}{2}},\hat{\theta}=X_{\frac{1}{2}}$

2.1理论基础

? ? ? ? 使用极大似然方法估计 $\lambda=\begin{pmatrix} \gamma\\ \theta \end{pmatrix}$ ,已知样本 $X=\{x_1,...,x_n\}$

$x_i \overset{i.i.d} {\sim}C(\gamma,\theta)$

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

$L(\lambda)=\prod_{i=1}^{n}\frac{1}{\pi}\frac{\gamma}{({x_i-\theta})^2+\gamma^2}$

$l(\lambda)=\sum_{i=1}^{n}ln(\frac{1}{\pi}\frac{\gamma}{({x_i-\theta})^2+\gamma^2})=-nln(\pi)+nln(\gamma)-\sum_{i=1}^{n}ln((x_i-\theta)^2+\gamma^2)$

$\dot{\l}(\lambda)=\bigl(\begin{smallmatrix} \frac{n}{\gamma}-\sum_{i=1}^{n}\frac{2\gamma}{(x_i-\theta)^2+\gamma^2}\\ \sum_{i=1}^{n}\frac{2(x_i-\theta)}{(x_i-\theta)^2+\gamma^2} \end{smallmatrix}\bigr)= \begin{pmatrix} 0\\0 \end{pmatrix}$

$\ddot{\l}(\lambda)=\begin{bmatrix} -\frac{n}{\gamma^2}-\sum_{i=1}^{n}\frac{2((x_i-\theta)^2-\gamma^2)}{((x_i-\theta)^2+\gamma^2)^2}&-\sum_{i=1}^{n}\frac{4\gamma(x_i-\theta)}{((x_i-\theta)^2+\gamma^2)^2} \\ -\sum_{i=1}^{n}\frac{4\gamma(x_i-\theta)}{((x_i-\theta)^2+\gamma^2)^2}&\sum_{i=1}^{n}\frac{2((x_i-\theta)^2-\gamma^2)}{((x_i-\theta)^2+\gamma^2)^2} \end{bmatrix}$

? ? ? ? 故迭代公式为

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

2.2算法

Step0:给定 $\hat{\theta}^{(0)}=\begin{pmatrix} quantile(X,0.75)-median(X)\\ median(X) \end{pmatrix}$ , $k=0$ ,停止精度 $\varepsilon=10^{-5}$

Step1:计算 $\dot{l}(\hat{\theta}^{(k)})$ ,如果

$||\dot{l}(\hat{\theta}^{(k)})||_2\leq \varepsilon$

?????? 则终止,否则进行下一步

Step2:计算 $\ddot{l}(\hat{\theta}^{(k)})$ ,令

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

Step3:计算 $\hat{\theta}^{(k)}=\hat{\theta}^{(k)}+d^{(k)},k=k+1$ ,?跳转Step1

2.2.1R语言实现

# 设立牛顿算法框架
Newtons = function(fun, x, theta,ep = 1e-5, it_max = 100){
  index = 0 ;k = 1
  while (k <= it_max){
    theta1 = theta
    obj = fun(x,theta)
    theta = theta - solve(obj$J,obj$f)
    norm = sqrt((theta - theta1) %*% (theta - theta1))
    if (norm < ep){
      index = 1; break
    }
    k = k + 1
  }
  obj = fun(x,theta)
  list(root = theta, it = k, index = index, FunVal = obj$f)
}
# 计算hessian矩阵和一阶导数
funs = function(x,theta){
  n = length(x)
  temp_1 = (x-theta[2])^2+theta[1]^2
  temp_2 = (x-theta[2])^2-theta[1]^2
  f = c(n/theta[1]-sum(2*theta[1]/temp_1), 
         sum((2*(x-theta[2]))/temp_1))
  J = matrix(c(-n/theta[1]^2-sum(2*temp_2/temp_1^2), -sum((4*theta[1]*(x-theta[2]))/temp_1^2),
                -sum((4*theta[1]*(x-theta[2]))/temp_1^2),sum(2*temp_2/temp_1^2)), nrow = 2, byrow = T)
  list(f = f, J = J)
}

# 抽取1000个样本
set.seed(80)
sample = rcauchy(1000,scale = 2,location = 2)

# 计算cauchy分布参数的分位数估计作为初始值
gamma_ = quantile(sample,0.75) - median(sample)
theta_ = median(sample)
theta = c(gamma_,theta_)

Newtons(funs,sample,theta)

2.2.2Python语言实现

import numpy as np
import scipy.stats as st

np.random.seed(12)
sample = st.cauchy.rvs(loc=1,scale = 1,size = 100) # scipy.stats.cauchy.rvs(loc=0, scale=1, size=1, random_state=None)
gamma_ = np.quantile(sample,0.75) - np.median(sample)
theta_ = np.median(sample)
theta = np.array([gamma_,theta_])

def Newtons(fun,x,theta,ep=1e-5,it_max = 100):
    index = 0
    k = 1
    while k <= it_max:
        theta1 = theta
        obj = fun(x,theta)
        theta = theta - np.dot(np.linalg.inv(obj[1]),obj[0])
        norm = np.sqrt(np.sum((theta - theta1)**2))
        if norm < ep:
            index = 1
            break
        k = k+1
    obj = fun(x,theta)
    print('gamma,theta估计值为%.3f,%.3f'%(theta[0],theta[1]))
def funs(x,theta):
    n = len(x)
    temp_1 = (x - theta[1]) ** 2 + theta[0] ** 2
    temp_2 = (x - theta[1]) ** 2 - theta[0] ** 2
    f = np.array([n/theta[0]-np.sum(2*theta[0]/temp_1),
                  np.sum((2*(x-theta[1]))/temp_1)])
    j = np.array([[-n/theta[0]**2-sum(2*temp_2/temp_1**2),-np.sum((4*theta[0]*(x-theta[1]))/temp_1**2)],
                  [-np.sum((4*theta[0]*(x-theta[1]))/temp_1**2),sum(2*temp_2/temp_1**2)]])
    return [f,j]

Newtons(funs,sample,theta)

3.Gamma 分布的极大似然估计

????????如果随机变量 $x$ 服从Gamma分布,记为 $x\sim{G(\alpha,\beta,\lambda)}$ ,其中 $\alpha$ 为形状参数(shape parameter),β 为尺度参数(scale parameter ),λ 为位置参数(location parameter)

$f(x|\alpha,\beta,\lambda)=\frac{1}{\beta^{\alpha}\Gamma(\alpha)}(x-\lambda)^{\alpha-1}e^{-\frac{x-\lambda}{\beta}},x\in(\lambda,\infty)$

? ? ? ? 为了给出Gamma分布三个参数的矩估计，现考虑分布的一二三阶原点矩(求解的技巧在于凑Gamma函数)

? ? ? ? 以一阶原点矩的证明为例，将x拆分为

$x=(x-\lambda)+\lambda$

故

$E(x)=\int_{\lambda}^{\infty}x\frac{1}{\beta^{\alpha}\Gamma(\alpha)}(x-\lambda)^{\alpha-1}e^{-\frac{x-\lambda}{\beta}}dx$ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

$=\frac{1}{\beta^{\alpha}\Gamma(\alpha)}\int_{\lambda}^{\infty}(x-\lambda+\lambda)(x-\lambda)^{\alpha-1}e^{-\frac{x-\lambda}{\beta}}dx$ ? ? ? ? ? ? ? ? ? ? ? ? ? ?

$=\frac{1}{\beta^{\alpha}\Gamma(\alpha)}\int_{\lambda}^{\infty}(x-\lambda)^{\alpha}e^{-\frac{x-\lambda}{\beta}}+\lambda(x-\lambda)^{\alpha-1}e^{-\frac{x-\lambda}{\beta}} dx$ ? ? ? ? ? ? ? ?

$=\frac{1}{\beta^{\alpha}\Gamma(\alpha)}\int_{\lambda}^{\infty}[\beta^{\alpha}(\frac{x-\lambda}{\beta})^{\alpha}e^{-\frac{x-\lambda}{\beta}}+\lambda(\frac{x-\lambda}{\beta})^{\alpha-1}e^{-\frac{x-\lambda}{\beta}}] d\frac{x-\lambda}{\beta}$

$=\frac{1}{\beta^{\alpha}\Gamma(\alpha)}\int_{0}^{\infty}[\beta^{\alpha}(t)^{\alpha}e^{-t}+\lambda(t)^{\alpha-1}e^{-t}] dt$ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??

$=\frac{1}{\beta^{\alpha}\Gamma(\alpha)}(\beta^{\alpha+1}\Gamma(\alpha+1)+\beta^{\alpha}\lambda\Gamma(\alpha))$ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??

$=\alpha\beta+\lambda$ ?①? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

? ? ? ? 对于二阶三阶原点矩，分别将 $x^2,x^3$ 拆分为

$x^2=(x-\lambda)^2+2\lambda(x-\lambda)+\lambda^2$

$x^3=(x-\lambda)^3+3\lambda(x-\lambda)^2+3\lambda^2(x-\lambda)+\lambda^3$

? ? ? ? 易得

$E(x^2)=\alpha^2\beta^2+\alpha\beta^2+2\alpha\beta\lambda+\lambda^2$

$E(x^3)=\alpha^3\beta^3+3\alpha^2\beta^3+2\alpha\beta^3+3\alpha^2\beta^2\lambda+3\alpha\beta^2\lambda+3\alpha\beta\lambda^2+\lambda^3$

? ? ? ? 又

$Var(x)=E[(x-E(x))^2]=E(x^2)-E^2(x)=\alpha\beta^2$ ②

$m_3=E[(x-E(x))^3]=E(x^3)-3E(x^2)E(x)+2E^3(x)=2\alpha\beta^3$ ③

? ? ? ? 由①②③可得

$\therefore\hat{\alpha}=\frac{4(Var(x))^3}{m_3^2},\hat{\beta}=\frac{m_3}{2Var(x)},\hat{\lambda}=mean(x)-2\frac{(Var(x))^2}{m_3}$

3.1理论基础

? ? ? ? 已知样本 $X=\{x_1,...,x_n\}$ ,参数 $\theta=\begin{bmatrix} \alpha\\ \beta \\ \lambda \end{bmatrix}$

$f(x|\theta)=\frac{1}{\beta^{\alpha}\Gamma(\alpha)}(x-\lambda)^{\alpha-1}e^{-\frac{x-\lambda}{\beta}}$

$L(\theta)=\prod_{i=1}^{n}\frac{1}{\beta^{\alpha}\Gamma(\alpha)}(x_i-\lambda)^{\alpha-1}e^{-\frac{x_i-\lambda}{\beta}}$

$l(\theta)=-n\alpha ln(\beta)-n ln(\Gamma(\alpha))+(\alpha-1)\sum_{i=1}^{n}ln(x_i-\lambda)-\frac{1}{\beta}\sum_{i=1}^{n}(x_i-\lambda)$

$\dot{l}(\theta)=\begin{pmatrix} -n ln(\beta)-n\frac{\Gamma'(\alpha)}{\Gamma(\alpha)}+\sum_{i=1}^{n}ln(x_i-\lambda)\\ -n\frac{\alpha}{\beta}+\frac{1}{\beta^2}\sum_{i=1}^{n}(x_i-\lambda) \\ (\alpha-1)\sum_{i=1}^{n}\frac{-1}{x_i-\lambda}+\frac{n}{\beta} \end{pmatrix}=\begin{pmatrix} 0\\ 0 \\ 0 \end{pmatrix}$

$\ddot{l}(\theta)=\begin{pmatrix} -n(\frac{\Gamma'(\alpha)}{\Gamma(\alpha)})' &-\frac{n}{\beta} &\sum_{i=1}^{n}\frac{-1}{x_i-\lambda} \\ -\frac{n}{\beta} &n\frac{\alpha}{\beta^2} -\frac{2}{\beta^3}\sum_{i=1}^{n}(x_i-\lambda) &-\frac{n}{\beta^2} \\ \sum_{i=1}^{n}\frac{-1}{x_i-\lambda} & -\frac{n}{\beta^2} & (\alpha-1)\sum_{i=1}^{n}\frac{-1}{(x_i-\lambda)^2} \end{pmatrix}$

? ? ? ? 故迭代公式为

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

3.2算法

Step0:给定 $\hat{\theta}^{(0)}=\begin{pmatrix} \frac{4(Var(x))^3}{m_3^2}\\ \frac{m_3}{2Var(x)}\\mean(x)-2\frac{(Var(x))^2}{m_3} \end{pmatrix}$ , $k=0$ ,停止精度 $\varepsilon=10^{-5}$

Step1:计算 $\dot{l}(\hat{\theta}^{(k)})$ ,如果

$||\dot{l}(\hat{\theta}^{(k)})||_2\leq \varepsilon$

?????? 则终止,否则进行下一步

Step2:计算 $\ddot{l}(\hat{\theta}^{(k)})$ ,令

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

Step3:计算 $\hat{\theta}^{(k)}=\hat{\theta}^{(k)}+d^{(k)},k=k+1$ ,?跳转Step1

3.2.1R语言实现

# 设立牛顿算法框架
Newtons = function(fun, x, theta,ep = 1e-5, it_max = 100){
  index = 0; k = 1
  while (k <= it_max ){
    theta1 = theta; obj = fun(x,theta)
    theta = theta - solve(obj$J,obj$f)
    norm = sqrt((theta - theta1) %*% (theta - theta1))
    if (norm < ep){
      index = 1; break
    }
    k = k + 1
  }
  obj = fun(x,theta)
  list(root = theta, it = k, index = index, FunVal = obj$f)
}
# 计算hessian矩阵和一阶导数
funs = function(x,theta){
  n = length(x)
  f = c(-n*log(theta[2])-n*digamma(theta[1])+sum(log(x-theta[3])), 
         -n*theta[1]/theta[2]+1/(theta[2]^2)*sum(x-theta[3]),
         (theta[1]-1)*sum(-1/(x-theta[3]))+n/theta[2])
  J = matrix(c(-n*trigamma(theta[1]), -n/theta[2], sum(-1/(x-theta[3])),
                -n/theta[2],n*theta[1]/(theta[2]^2)-2/(theta[2]^3)*sum(x-theta[3]),-n/(theta[2]^2),
                sum(-1/(x-theta[3])),-n/(theta[2]^2),(theta[1]-1)*sum(-1/(x-theta[3])^2)), nrow = 3, byrow = T)
  list(f = f, J = J)
}

# 抽取随机gamma
set.seed(80)
sample = rgamma(100,2,)

# 计算gamma分布参数的矩估计作为初始值
s_m = mean(sample)
s_v = var(sample)
m_3 = (sum((sample-s_m)^3))/length(sample)
alpha = (4*(s_v)^3)/(m_3^2)
beta = m_3/(2*s_v)
lambda = s_m-2*((s_v)^2/m_3)
theta = c(alpha,beta,lambda)

Newtons(funs,sample,theta)

3.2.2Python语言实现

import numpy as np
import scipy.stats as st
import scipy.special as sp

np.random.seed(12)
sample = st.gamma.rvs(2,scale = 1,size = 50) # scipy.stats.gamma.rvs(a, loc=0, scale=1, size=1, random_state=None)
s_m = sample.mean()
s_v =sample.var()
m_3 = np.mean(((sample - s_m)**3))
alpha = (4*(s_v)**3)/(m_3**2)
beta = m_3/(2*s_v)
lambda_ = s_m-2*((s_v)**2/m_3)
theta = np.array([alpha,beta,lambda_])

def Newtons(fun,x,theta,ep=1e-5,it_max = 100):
    index = 0
    k = 1
    while k <= it_max:
        theta1 = theta
        obj = fun(x,theta)
        theta = theta - np.dot(np.linalg.inv(obj[1]),obj[0])
        norm = np.sqrt(np.sum((theta - theta1)**2))
        if norm < ep:
            index = 1
            break
        k = k+1
    obj = fun(x,theta)
    print('alpha,beta,lambda估计值为%.3f,%.3f,%.3f'%(theta[0],theta[1],theta[2]))
def funs(x,theta):
    n = len(x)
    f = np.array([-n*np.log(theta[1])-n*sp.digamma(theta[0])+np.sum(np.log(x-theta[2])),
                  -n*theta[0]/theta[1]+1/(theta[1]**2)*np.sum(x-theta[2]),
                  (theta[0]-1)*np.sum(-1.0/(x-theta[2]))+n/theta[1]])
    j = np.array([[-n*sp.polygamma(1,theta[0]),-n/theta[1],np.sum(-1/(x-theta[2]))],
                  [-n/theta[1],n*theta[0]/(theta[1]**2)-2/(theta[1]**3)*np.sum(x-theta[2]),-n/(theta[1]**2)],
                  [np.sum(-1/(x-theta[2])),-n/(theta[1]**2),(theta[0]-1)*np.sum(-1/(x-theta[2])**2)]])
    return [f,j]

Newtons(funs,sample,theta)

注意：对Gamma分布三参数的极大似然估计过程中,如果使用牛顿法,很容易出现矩阵不正定的情况从而无法得到正确解,这个时候可以使用拟牛顿法或者修正的牛顿法。

人工智能最新文章

2022吴恩达机器学习课程——第二课（神经网

第十五章规则学习

FixMatch: Simplifying Semi-Supervised Le

数据挖掘Java——Kmeans算法的实现

大脑皮层的分割方法

【翻译】GPT-3是如何工作的

论文笔记:TEACHTEXT: CrossModal Generaliz

python从零学（六）

详解Python 3.x 导入(import)

【答读者问27】backtrader不支持最新版本的