在Statlect上搜索概率和统计术语
统计列克特
指数 > 的基本原理 统计 > 最大似然 > 算法

最大似然-MATLAB示例

通过 博士

在演讲中 最大 可能性-算法 我们已经解释了如何计算最大值 通过数值方法对参数的似然估计。在本讲座中,我们 提供真人在线斗地主完整的示例,说明如何使用MATLAB进行操作。

在接下来的内容中,我将假定您有权访问MATLAB安装 同时包含统计信息和优化工具箱(安装 的 八度 -免费软件 与MATLAB非常相似-应该提供相同的功能)。

目录

数据

我们有100个独立抽奖的样本 标准学生t分布n 自由程度。参数 n 未知,我们希望通过最大可能性进行估算。

数据(100个观测值)存储在MATLAB文件中 数据垫,您需要下载。

参数化

注意参数 n 必须严格为正,即它必须属于区间 [eq1]. 因此,我们需要解决的优化问题才能进行估算 n 是真人在线斗地主约束优化问题。如演讲中所述 最大似然 - 算法,最好尽可能避免约束问题。 在这种情况下,可能是因为 n 可以很容易地重新设置参数 如 [eq2]哪里 $ heta $ 是我们的新参数,对此没有任何限制,因为它可以 间隔中的任何值 [eq3].

对数似然函数编码

似然函数被编码为真人在线斗地主例程,该例程以的值作为输入 参数和数据,并返回输出的值 对数似然及其符号已更改。代码如下。

function val=log_lik(theta,data)
n=exp(theta);
val=-sum(log(tpdf(data,n)));

的 name 的 the function 是 log_lik. It takes 如 arguments the parameter theta 和 the vector 的 observations data.

的 function tpdf (which 是 part 的 the 统计istics 工具箱)计算标准学生t的概率密度函数 distribution. In particular, tpdf(data,n) returns a 密度向量(向量中每个观测值的真人在线斗地主密度 data),在假设度数下 的 freedom 是 equal 至 n. 通过 taking the natural logarithm 与 log function 和 summing over all 向量的项,我们获得样本的对数似然。其他 words, 与 command sum(log(tpdf(data,df))) we 计算 对数似然[eq4]哪里 $ x_ {j} $ 是 an observation (a component 的 the vector data), $ T $ 是个 sample size (the dimension 的 the vector data) 和 [eq5] 是个 可能性 密度函数$ x_ {j} $, 给定参数等于 $ heta $.

最后,我们在前面加上减号来更改对数似然的符号 因为我们将使用的优化例程执行 默认情况下最小化,而我们想最大化对数似然 功能。这样获得的值存储在变量中 val,由函数返回。

优化算法

本小节中提供的代码运行优化算法,其中 为了从数字上找到参数的最大似然估计。

load data
options=optimset('Display','off','MaxIter',10000,'TolX',10^-30,'TolFun',10^-30);
theta0=1;
[theta,fval,exitflag,output,grad,hessian]=fminunc('log_lik',theta0,options,data);
exp(theta)

We first load the data file 与 command load data.

然后,我们设置优化算法的一些选项。选项 Display 是 set 至 的f, which 表示优化算法将以静默方式运行,而不会显示 output 的 each iteration. 的 option MaxIter 是 set 至 10000,这表示该算法将执行 maximum 的 10,000 iterations. 的 option TolX 是 set 至 10^-30, which means that the termination 至lerance 在参数上将是 $ 10 ^ {-30} $. Also the option TolFun 是 set 至 10^-30,表示终止公差为 要优化的功能的价值将是 $ 10 ^ {-30} $.

参数的初始值(我们的初始猜测)设置为等于1 与 command theta0=1.

最后,基于MATLAB导数的优化算法 fminunc 叫做。它作为参数:

该函数作为输出返回:

的 last line 的 code (exp(theta)) displays the 参数估计 $ n,$, 等于 至[eq6]哪里 $ widehat {heta} $ 是的最大似然估计 $ heta $ (stored in the variable theta)。 If the code 是 正确执行,显示值应为3.9389,这是我们的观点 感兴趣参数的估计。这对应于 1.3709为 $ heta $.

您还可以显示渐变的值(通过输入 grad)。从数字上接近零的数字 应该显示,这意味着最优的一阶条件 (即一阶导数等于零)近似满足。

多次开始

我们不知道 fminunc 算法,以及我们的似然函数的属性,可以保证 数值收敛才能真正解决问题。由于我们有 前面已经解释过,当没有理论保证时,数值 可以实现收敛, 多次开始 方法是 通常遵循:数值优化算法运行几次, 包含参数的随机起始值,以及 算法(或大多数算法)得出相同的建议解决方案,然后 这可以作为证据证明拟议的解决方案是真人在线斗地主很好的近似值 真正的解决方案。

以下代码是上述代码的多次启动变体。

load data
stream=RandStream('mt19937ar','Seed',1);
RandStream.setDefaultStream(stream);
options=optimset('Display','off','MaxIter',10000,'TolX',10^-30,'TolFun',10^-30);
for j=1:10
theta0=1+10*randn(1,1);
[theta(j),fval(j),exitflag(j),output(j),grad(j),hessian(j)]=...
fminunc('log_lik',theta0,options,data);
end
theta

的 two commands involving RandStream initialize the 随机数生成器。您不必担心这两个命令。他们 只是为了确保如果您在计算机上运行此代码,您将 得到与我得到的结果完全相同的结果。

重复优化十次,每次都是初始值 theta0 是随机生成的。特别是,我们添加到 our initial guess (theta0=1) a random number extracted 从真人在线斗地主 标准正态分布 (randn(1,1)), multiplied 通过 10.

输出存储在向量中,向量的条目由变量索引 j。特别是提出的十个解决方案是 stored in the vector theta, 和 the ten values 的 the 与这些解相对应的对数似然函数存储在 vector val.

如果在运行脚本后通过键入以下内容打印建议的解决方案 theta,您将看到它们都等于真人在线斗地主 我们之前发现(1.3709),除了真人在线斗地主-第二个-等于 至12.8117。因此,我们遇到了一些麻烦。这证实了我们 怀疑只满足一次跑步是不明智的 算法。

通过 printing the vectors valgrad,我们可以很容易地看到12.8117不能是 解,因为它对应于目标函数的值 高于对应于1.3709的值,并且因为梯度不是 等于零。

如果我们多次运行该算法(例如,增加 repetitions 至 1000 in the for cycle), we will see 还提出了其他解决方案,但这并不是真正的解决方案, 由于相同的原因(即非零导数, 目标函数)。

因此,我们可以将多次启动算法的结果作为证明 1.3709是真人在线斗地主声音解决方案。同时,这使我们想起了真人在线斗地主事实, 如果我们只运行一次优化算法,那确实是真人在线斗地主坏主意 不确定其收敛性。

估计方差

产生对方差的估计 $ widehat {heta} $, 我们可以使用讲座中介绍的方法之一 最大 似然-协方差矩阵估计。特别是因为 fminunc 提供了粗麻布的数值估计 矩阵,我们可以使用基于此估计的方法。

请记住,最大似然估计器的分布 $ widehat {heta} $ 可以由均值等于的多元正态分布近似 真实参数 $ heta _ {0} $ 和协方差矩阵等于 [eq7]哪里[eq8]是 渐近协方差矩阵的估计 [eq9] 表示二阶导数矩阵。换句话说, 的方差 $ widehat {heta} $[eq10]

在我们的例子中,只有真人在线斗地主参数,所以矩阵实际上是真人在线斗地主标量 我们可以 写[eq11]的 value stored 通过 fminunc in the variable hessian 是...的二阶导数 log_lik,这是(请参见 以上)[eq12]通过 线性,这是相同的 如 [eq13]因此, 我们需要做的所有事情以获得对方差的估计 $ widehat {heta} $ 是 至 计算 reciprocal 的 hessian, 与 command hessian^-1 (if you do this after running the 多次启动算法,仅取向量的一项 hessian,对应于估算值1.3709)。如果 一切都正确完成,估计值应等于0.0970。这个 可以转化为对方差的估计 $ widehat {n} $ 增量法,乘以 的估计方差 $ widehat {heta} $ 通过 [eq14]. 结果,我们得到估计的方差 $ widehat {n} $ 是1.5046。

MATLAB文件

本讲座中介绍的所有MATLAB代码都存储在 压缩文件,你 可以下载。

如何引用

请引用为:

Taboga, Marco (2017). "最大似然-MATLAB示例", 列克特ures on 可能性 theory 和 mathematical 统计, Third edition. Kindle Direct Publishing. Online appendix. //www.junruiqiche.com/fundamentals-of-statistics/maximum-likelihood-matlab-example.

这本书

该网站上提供的大多数学习材料现在都以传统教科书格式提供。