解釋模型如何工作是數據科學中最為基本最為關鍵的問題之一。當你建立了一個模型之後，它給了你預期的結果，但是它背後的過程是什麼呢？作為一個數據科學家，你需要對這個經常被問到的問題做出解答。

例如，假設您建立了一個預測公司股價的模型。您注意到夜深人靜的時候，股票價格上漲得很快。背後可能有多種原因，找出可能性最大的原因便是最大似然估計的意義所在。這一概念常被用于經濟學、MRIs、衛星成像等領域。

來源：YouTube

在這篇文章中，我們将研究最大似然估計(以下簡稱MLE)是如何工作的，以及它如何用于确定具有任何分布的模型的系數。理解MLE将涉及到概率和數學，但我将嘗試通過例子使它更通俗易懂。

. 注：如前所述，本文假設您已經了解概率論的基本知識。您可以通過閱讀這篇文章來澄清一些基本概念：每個數據科學專業人員都應該知道的概率分布常識。

• 為什麼要使用最大似然估計(MLE)？
• 通過一個實例了解MLE
• 進一步了解技術細節分布參數似然對數似然最大似然估計
• 利用MLE确定模型系數
• R語言的MLE實現

假設我們想預測活動門票的銷售情況。數據的直方圖和密度如下。

你将如何為這個變量建模？該變量不是正态分布的，而且是不對稱的，因此不符合線性回歸的假設。一種常用的方法是對變量進行對數、平方根（sqrt）、倒數等轉換，使轉換後的變量服從正态分布，并進行線性回歸建模。

讓我們試試這些轉換，看看結果如何：

對數轉換：

平方根轉換：

倒數轉換：

所有這些都不接近正态分布，那麼我們應該如何對這些數據進行建模，才能不違背模型的基本假設？如何利用正态分布以外的其他分布來建模這些數據呢？如果我們使用了不同的分布，又将如何來估計系數？

這便是最大似然估計(MLE)的主要優勢。

在研究統計和概率時，你肯定遇到過諸如x＞100的概率，因為x服從正态分布，平均值為50，标準差為10。在這些問題中，我們已經知道分布（在這種情況下是正态分布）及其參數（均值和标準差），但在實際生活問題中，這些參數是未知的，并且必須從數據中估計出來。MLE可以幫助我們确定給定數據的分布參數。

讓我們用一個例子來加深理解：假設我們用數據來表示班級中學生的體重(以kg為單位)。數據如下圖所示(還提供了用于生成數據圖的R代碼)：

圖 1

x = as.data.frame(rnorm(50,50,10))

ggplot(x, aes(x = x)) geom_dotplot()

這似乎遵循正态分布。但是我們如何得到這個分布的均值和标準差呢？一種方法是直接計算給定數據的平均值和标準差，分别為49.8公斤和11.37公斤。這些值能很好地表示給定的數據，但還不能最好地描述總體情況。

我們可以使用MLE來獲得更穩健的參數估計。因此，MLE可以定義為從樣本數據中估計總體參數(如均值和方差、泊松率(Lambda)等)的方法，從而使獲得觀測數據的概率(可能性)最大化。

為了加深對MLE的理解，嘗試猜測下列哪一項會使觀察上述數據的概率最大化？

1. Mean = 100, SD = 10 均值=100，标準差=10
2. Mean = 50, SD = 10 均值=50 ，标準差=10

顯然，如果均值為100，我們就不太可能觀察到上述數據分布圖形。

知道MLE能做什麼之後，我們就可以深入了解什麼是真正的似然估計，以及如何對它最大化。首先，我們從快速回顧分布參數開始。

首先，來了解一下分布參數。維基百科對這個詞的定義如下：“它是一個概率分布的量化指數”，可以将它視為樣本總數的數值特征或一個統計模型。通過下面的圖表來理解它：

鐘形曲線的寬度和高度由兩個參數決定-均值和方差。這就是正态分布的分布參數。同樣，泊松分布由一個參數lambda控制，即事件在時間或空間間隔内發生的次數

大多數分布都有一個或兩個參數，但有些分布可以有多達4個參數，比如4參數β分布。

從圖2和圖3中我們可以看到，給定一組分布參數，一些數據值比其他數據的概率更大。從圖1中，我們已經看到，當平均值是50而不是100時，給定的數據更有可能發生。然而，在現實中，我們已經觀察到了這些數據。因此，我們面臨着一個逆向問題：給定觀測數據和一個感興趣的模型，我們需要在所有概率密度中找到一個最有可能産生數據的概率密度函數/概率質量函數(f(x_\θ)。

為解決這一逆向問題，我們通過逆轉f(x=θ)中數據向量x和(分布)參數向量θ來定義似然函數，即，

L(θ;x) = f(x| θ)

在MLE中，可以假定我們有一個似然函數L(θ；x)，其中θ是分布參數向量，x是觀測集。我們感興趣的是尋找具有給定觀測值(x值)的最大可能性的θ值。

如果假設觀測集(Xi)是獨立的同分布随機變量，概率分布為f0(其中f0=正态分布，例如圖1)，那麼手頭的數學問題就變得簡單了。似然函數可以簡化為：

為了求這個函數的極大值/極小值，我們可以取此函數w.r.tθ的導數，并将其設為0(因為零斜率表示最大值或極小值)。因為我們這裡有乘積，所以需要應用鍊式規則，這對乘積來說是相當麻煩的。一個聰明的訣竅是取似然函數的對數，并使其最大化。這将乘積轉換為加法，并且由于對數是一個嚴格遞增的函數，因此不會影響θ的結果值。所以我們有：

為找到對數似然函數LL（Th；x）的極大值，我們可以：

• 取LL(θ；x)函數w.r.tθ的一階導數，并将其等價于0，
• 取LL(θ；x)函數w.r.tθ的二階導數，并确認其為負值。

在許多情況下，微積分對最大化似然估計沒有直接幫助，但最大值仍然可以很容易地識别出來。在尋找最大對數似然值的參數值時，沒有任何東西比一階導數等于零具有更為 “優先”或特殊的位置。當需要估計一些參數時，它僅僅是一個方便的工具而已。

在通常情況下，能對函數求最大值的argmax的方法都可能适合于尋找log似然函數的最大值。這是一個無約束的非線性優化問題.我們尋求一種優化算法以下列方式工作：

1. 從任意起始點可靠地收斂到局部最小化
2. 速度盡可能快

使用優化技術來最大化似然是非常常見的；可以有多種方法來實現(比如：牛頓法、Fisher評分法、各種基于共轭梯度的方法、最陡下降法、Nder-Mead型(單純形)方法、BFGS法和多種其他技術)。

結果表明，當模型假設為高斯時，MLE估計等價于一般的最小二乘法。

你可以參考這裡來證明它。

現在讓我們看看如何使用MLE來确定預測模型的系數。

假設我們有一個n個觀測量y1，y2，…，yn的樣本，它們可以被看作是獨立的泊松随機變量：Yi ∼ P(µi)。此外，假設我們希望讓均值i(方差也同樣！)依賴于變量xi組成的向量。我們可以構成如下簡單的線性模型-

其中θ是模型系數的向量。這個模型的缺點是右邊的線性預測器可以假定為任何實際值，而表示期望計數的左側泊松均值必須是非負的。這個問題的一個簡單的解決辦法是用線性模型來模拟平均值的對數。因此，我們考慮了一個具有對數鍊接對數的廣義線性模型，它可以寫成如下所示：

我們的目的是利用MLE找到θ：

現在，泊松分布如下：

利用在上一節中學到的對數似然概念來求θ。取上述方程的對數，忽略包含log(y！)的常數，我們得到的對數似然函數是：

其中，µi依賴協變量xi和θ系數的向量。我們可以用變元µi = exp(xi’θ)代替，求解方程，得到最大似然值θ。得到了θ向量之後，我們就可以通過乘以xi和θ向量來預測平均值的期望值。

在本節中，我們将采用一個真實的數據集，運用前面學到的概念，來解決一個問題。您可以從此鍊接下載數據集。數據集中的示例如下：

售票日期計票

25-08-2012 00:00 8

25-08-2012 01:00 2

25-08-2012 02:00 6

25-08-2012 03:00 2

25-08-2012 04:00 2

25-08-2012 05:00 2

它有從2012年8月25日到2014年9月25日每小時售出的門票數量(約18K記錄)。我們的目的是預測每小時售出的門票數量。這是本文第一節中讨論的同一個數據集。

這個問題可以用回歸、時間序列等技術來解決。在這裡，我們将使用我們已經學習到的統計建模技術，用R語言實現。

首先，分析一下數據。在統計建模中，我們更關心的是目标變量是如何分布的。讓我們看一看計數的分布：

hist(Y$Count, breaks = 50,probability = T ,main = "Histogram of Count Variable")

lines(density(Y$Count), col="red", lwd=2)

這可以看作是泊松分布，或者我們甚至可以嘗試拟合指數分布。

由于手頭的變量是票的計數，泊松分布是一個更合适的模型。指數分布通常用于模拟事件之間的時間間隔。

讓我們來計算一下這兩年售出的票的數量：

看上去随着時間的推移，門票的銷售有了很大的增長。為了将問題簡單化，讓我們僅将時間作為一個因子來建模，其中時間定義為2012年8月25日以來幾個星期。我們可以把它寫成：

其中，µ（售出的票數）是泊松分布的平均值，而θ0和θ1是我們需要估計的系數。

組合方程1和2，我們得到如下的對數似然函數：

我們可以使用Rstats 4包中的mle() 函數來估計系數θ0和θ1。它需要下列主要參數：

1. 需要最小化的負似然函數：這與我們剛剛導出的函數相同，但前面有一個負号[最大對數似然度等于最小化負對數似然]。
2. 系數向量的起點：這是系數的初始預測值。結果可以根據這些值而變化，因為函數可能達到局部最小值。因此，通過運行不同起點的函數來驗證結果是很好的辦法。
3. BFGS是默認的對似然函數進行優化的方法。

在我們的示例中，負對數似然函數的編碼如下：

nll <- function(theta0,theta1) {

x <- Y$age[-idx]

y <- Y$Count[-idx]

mu = exp(theta0 x*theta1)

-sum(y*(log(mu)) - mu)

}

我将數據分為訓練集和測試集，以便客觀地評價模型的性能。idx是測試集中行的指數。

set.seed(200)

idx <- createDataPartition(Y$Count, p=0.25,list=FALSE)

接下來，調用mle函數來獲得參數：

est <- stats4::mle(minuslog=nll, start=list(theta0=2,theta1=0))

summary(est)

Maximum likelihood estimation

Call:

stats4::mle(minuslogl = nll, start = list(theta0 = 2, theta1 = 0))

Coefficients:

Estimate Std. Error

theta0 2.68280754 2.548367e-03

theta1 0.03264451 2.998218e-05

-2 log L: -16594396

我們得到了系數的估計值，利用RMSE作為獲得測試集結果的評估度量：

pred.ts <- (exp(coef(est)['theta0'] Y$age[idx]*coef(est)['theta1'] ))

rmse(pred.ts, Y$Count[idx])

86.95227

現在，讓我們看看我們的模型和标準線性模型(正态分布的誤差)的對比，本模型是用對數計數建模。

lm.fit <- lm(log(Count)~age, data=Y[-idx,])

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 1.9112992 0.0110972 172.2 <2e-16 ***

age 0.0414107 0.0001768 234.3 <2e-16 ***

pred.lm <- predict(lm.fit, Y[idx,])

rmse(exp(pred.lm), Y$Count[idx])

93.77393

可以看出來，标準線性模型的RMSE比我們的泊松分布模型要高。讓我們比較這兩個模型在樣本上的殘差圖，看看這些模型在不同的區域中表現如何：

與常規線性回歸相比，泊松回歸的誤差更接近于零。

在Python中，也可以通過使用scipy.optimize.minimize()函數來實現目标函數的最小化，對初始值的估計同BFGS、L-BFGS等參數和方法類似。

在R語言中，使用stats包中的glm 函數建模更加簡單。它支持泊松，伽瑪，二項分布，Quasi,，逆高斯，拟二項分布，拟泊松分布等等。對于上面所示的示例，可以使用以下命令直接獲得系數：

glm(Count ~ age, family = "poisson", data = Y)

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) 2.669 2.218e-03 1203 <2e-16 ***

age 0.03278 2.612e-05 1255 <2e-16 ***

在Python中也可以使用pymc.glm()函數，并設置為pm.glm.Familes.Poisson()系列。

對上述例子的一種思考是，參數空間中是否存在比标準線性模型估計更好的系數。正态分布是缺省分布，也是最廣泛使用的分布形式，但如果采用其它更為正确的分布，則可以得到更好的結果。最大似然估計是一種可以用于估計分布參數而不考慮所使用的分布的技術。因此，下次當您手頭有建模問題時，首先看看數據的分布情況，看看有沒有比正态分布更有意義的分布！

詳細的代碼和數據在我的Gizub存儲庫中可以找到。有關使用年齡變量的數據讀取、格式化和建模的示例，請參閱“ModelingSingleVariables.R”文件。此外，我還使用了多個變量進行建模，它保存于“ModelingMultipleVariables.R”文件中。

更多精彩资讯请关注tft每日頭條，我们将持续为您更新最新资讯!