怎么用r生成双参数指数分布的数据（一文读懂最大似然估计）

怎么用r生成双参数指数分布的数据（一文读懂最大似然估计）例如，假设您建立了一个预测公司股价的模型。您注意到夜深人静的时候，股票价格上涨得很快。背后可能有多种原因，找出可能性最大的原因便是最大似然估计的意义所在。这一概念常被用于经济学、MRIs、卫星成像等领域。 本文将研究MLE是如何工作的，以及它如何用于确定具有任何分布的模型的系数。简介解释模型如何工作是数据科学中最为基本最为关键的问题之一。当你建立了一个模型之后，它给了你预期的结果，但是它背后的过程是什么呢？作为一个数据科学家，你需要对这个经常被问到的问题做出解答。

作者：阿尼·辛格

翻译： 陈之炎

校对：丁楠雅

本文约4200字，建议阅读10 分钟。

本文将研究MLE是如何工作的，以及它如何用于确定具有任何分布的模型的系数。

简介

解释模型如何工作是数据科学中最为基本最为关键的问题之一。当你建立了一个模型之后，它给了你预期的结果，但是它背后的过程是什么呢？作为一个数据科学家，你需要对这个经常被问到的问题做出解答。

例如，假设您建立了一个预测公司股价的模型。您注意到夜深人静的时候，股票价格上涨得很快。背后可能有多种原因，找出可能性最大的原因便是最大似然估计的意义所在。这一概念常被用于经济学、MRIs、卫星成像等领域。

来源：YouTube

在这篇文章中，我们将研究最大似然估计(以下简称MLE)是如何工作的，以及它如何用于确定具有任何分布的模型的系数。理解MLE将涉及到概率和数学，但我将尝试通过例子使它更通俗易懂。

注：如前所述，本文假设您已经了解概率论的基本知识。您可以通过阅读这篇文章来澄清一些基本概念：《每个数据科学专业人员都应该知道的概率分布常识》：https://www.analyticsvidhya.com/blog/2017/09/6-probability-distributions-data-science/

目录

• 为什么要使用最大似然估计(MLE)？
• 通过一个实例了解MLE
• 进一步了解技术细节
• 分布参数
• 似然
• 对数似然
• 最大似然估计
• 利用MLE确定模型系数
• R语言的MLE实现

为什么要使用最大似然估计(MLE)？

假设我们想预测活动门票的销售情况。数据的直方图和密度如下。

你将如何为这个变量建模？该变量不是正态分布的，而且是不对称的，因此不符合线性回归的假设。一种常用的方法是对变量进行对数、平方根（sqrt）、倒数等转换，使转换后的变量服从正态分布，并进行线性回归建模。

让我们试试这些转换，看看结果如何：

• 对数转换：

• 平方根转换：

• 倒数转换：

所有这些都不接近正态分布，那么我们应该如何对这些数据进行建模，才能不违背模型的基本假设？如何利用正态分布以外的其他分布来建模这些数据呢？如果我们使用了不同的分布，又将如何来估计系数？

这便是最大似然估计(MLE)的主要优势。

举一个例子来加深对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. 均值=100，标准差=102. 均值=50 ，标准差=10

显然，如果均值为100，我们就不太可能观察到上述数据分布图形。

进一步了解技术细节

知道MLE能做什么之后，我们就可以深入了解什么是真正的似然估计，以及如何对它最大化。首先，我们从快速回顾分布参数开始。

• 分布参数

首先，来了解一下分布参数。维基百科对这个词的定义如下：“它是一个概率分布的量化指数”，可以将它视为样本总数的数值特征或一个统计模型。通过下面的图表来理解它：

图 2

钟形曲线的宽度和高度的两个参数决定均值和方差。这就是正态分布的分布参数。同样，泊松分布由一个参数lambda控制，即事件在时间或空间间隔内发生的次数。

图 3

大多数分布都有一个或两个参数，但有些分布可以有多达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似然函数的最大值。这是一个无约束的非线性优化问题。我们寻求一种优化算法以下列方式工作：

• 从任意起始点可靠地收敛到局部最小化
• 速度尽可能快

使用优化技术来最大化似然是非常常见的，可以有多种方法来实现(比如：牛顿法、Fisher评分法、各种基于共轭梯度的方法、最陡下降法、Nder-Mead型(单纯形)方法、BFGS法和多种其他技术)。

结果表明，当模型假设为高斯时，MLE估计等价于一般的最小二乘法。

你可以参考下面这篇文章来证明它：

link：http://people.math.gatech.edu/~ecroot/3225/maximum_likelihood.pdf

用MLE确定模型系数

现在让我们看看如何使用MLE来确定预测模型的系数。

假设我们有一个n个观测量y1，y2，…，yn的样本，它们可以被看作是独立的泊松随机变量：Yi ∼ P(µi)。此外，假设我们希望让均值i(方差也同样！)依赖于变量xi组成的向量。我们可以构成如下简单的线性模型：

中θ是模型系数的向量。这个模型的缺点是右边的线性预测器可以假定为任何实际值，而表示期望计数的左侧泊松均值必须是非负的。这个问题的一个简单的解决办法是用线性模型来模拟平均值的对数。因此，我们考虑了一个具有对数链接对数的广义线性模型，它可以写成如下所示：

我们的目的是利用MLE找到θ。

现在，泊松分布如下：

利用在上一节中学到的对数似然概念来求θ。取上述方程的对数，忽略包含log(y！)的常数，我们得到的对数似然函数是：

其中，µi依赖协变量xi和θ系数的向量。我们可以用变元µi = exp(xi’θ)代替，求解方程，得到最大似然值θ。得到了θ向量之后，我们就可以通过乘以xi和θ向量来预测平均值的期望值。

用R语言实现MLE

在本节中，我们将采用一个真实的数据集，运用前面学到的概念，来解决一个问题。您可以从此链接下载数据集：

https://s3-ap-south-1.amazonaws.com/av-blog-media/wp-content/uploads/2018/07/ Train_Tickets.csv

数据集中的示例如下：

售票日期 计票

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。它需要下列主要参数：

• 需要最小化的负似然函数：这与我们刚刚导出的函数相同，但前面有一个负号[最大对数似然度等于最小化负对数似然]。
• 系数向量的起点：这是系数的初始预测值。结果可以根据这些值而变化，因为函数可能达到局部最小值。因此，通过运行不同起点的函数来验证结果是很好的办法。
• 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 estimationCall:stats4::mle(minuslogl = nll start = list(theta0 = 2 theta1 = 0)) Coefficients: Estimate Std. Errortheta0 2.68280754 2.548367e-03theta1 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”文件中。

原文标题：

An Introductory Guide to Maximum Likelihood Estimation (with a case study in R)

原文链接：

https://www.analyticsvidhya.com/blog/2018/07/introductory-guide-maximum-likelihood-estimation-case-study-r/

译者简介

陈之炎，北京交通大学通信与控制工程专业毕业，获得工学硕士学位，历任长城计算机软件与系统公司工程师，大唐微电子公司工程师，现任北京吾译超群科技有限公司技术支持。目前从事智能化翻译教学系统的运营和维护，在人工智能深度学习和自然语言处理（NLP）方面积累有一定的经验。业余时间喜爱翻译创作，翻译作品主要有：IEC-ISO 7816、伊拉克石油工程项目、新财税主义宣言等等，其中中译英作品“新财税主义宣言”在GLOBAL TIMES正式发表。能够利用业余时间加入到THU 数据派平台的翻译志愿者小组，希望能和大家一起交流分享，共同进步。