目錄

0.前言

1.理論基礎

2.Cauchy分布的極大似然估計

2.1理論基礎

2.2算法

2.2.1R語言實現

2.2.2Python語言實現

3.Gamma 分布的極大似然估計

3.1理論基礎

3.2算法

3.2.1R語言實現

3.2.2Python語言實現

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.理論基礎

? ? ? ? 本節給出牛頓法求分布的極大似然參數估計的一般理論。

????????如果隨機變量?獨立同分布于,且已知一組樣本? ,為了估計該分布的參數，可以使用極大似然估計的方法。

?????? 首先寫出樣本的似然函數

?????? 對? 進行對數化處理,得到對數似然函數

?????? 則求解未知參數等價于求解以下等式方程組

?????? 不妨假設收斂解為? ,將在? 的鄰域內展開成泰勒級數得

?????? 這樣就得到一個迭代關系式

?????? 如果是連續的，并且待求的零點是孤立的，那么在零點周圍存在一個區域，只要初始值位于這個鄰近區域內，那么牛頓法必定收斂。 并且，如果不為0, 那么牛頓法將具有平方收斂的性能。 粗略地說，這意味著每迭代一次，牛頓法結果的有效數字將增加一倍。

2.Cauchy分布的極大似然估計

????????如果隨機變量服從柯西分布,記為 ,其中為最大值一半處的一半寬度的尺度參數(scale parameter ),為定義分布峰值位置的位置參數(location parameter )

?????? 當 ,此時的Cauchy 分布稱為標準Cauchy 分布

?????? Cauchy 分布的最特別的性質是其期望及高階矩都不存在,自然也就無法對參數進行矩估計。但Cauchy分布的cdf具有很好的性質,可以利用一組樣本的分位點來對參數進行點估計。

2.1理論基礎

? ? ? ? 使用極大似然方法估計,已知樣本

? ? ? ? 故迭代公式為

2.2算法

Step0:給定,,停止精度

Step1:計算 ,如果

?????? 則終止,否則進行下一步

Step2:計算 ,令

Step3:計算,?跳轉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 npimport scipy.stats as stnp.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 分布的極大似然估計

????????如果隨機變量 服從Gamma分布,記為,其中為形狀參數(shape parameter),β 為尺度參數(scale parameter ),λ 為位置參數(location parameter)

? ? ? ? 為了給出Gamma分布三個參數的矩估計，現考慮分布的一二三階原點矩(求解的技巧在于湊Gamma函數)

? ? ? ? 以一階原點矩的證明為例，將x拆分為

故

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

? ? ? ? ? ? ? ? ? ? ? ? ? ?

? ? ? ? ? ? ? ?

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??

?①? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

? ? ? ? 對于二階三階原點矩，分別將拆分為

? ? ? ? 易得

? ? ? ? 又

②

③

? ? ? ? 由①②③可得

3.1理論基礎

? ? ? ? 已知樣本,參數

? ? ? ? 故迭代公式為

3.2算法

Step0:給定,,停止精度

Step1:計算 ,如果

?????? 則終止,否則進行下一步

Step2:計算 ,令

Step3:計算,?跳轉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)}# 抽取隨機gammaset.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 npimport scipy.stats as stimport scipy.special as spnp.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分布三參數的極大似然估計過程中,如果使用牛頓法,很容易出現矩陣不正定的情況從而無法得到正確解,這個時候可以使用擬牛頓法或者修正的牛頓法。