MATLAB最大似然估计中的数值问题

Question

我正在使用mle 和mlecov 来估计标量噪声信号n 的均值和方差，该信号假定为正态分布，具有以下均值和标准差模型：

mean(x,y) = @(x,y) k(1)+k(2)*x+k(3)*x.^2+k(4)*y+k(5)*y.^2;
sd(x,y)  = @(x,y) k(6)+k(7)*x+k(8)*x.^2+k(9)*y+k(10)*y.^2;

其中x 在 [0,3] 区间内，y 在 [0,pi/2] 区间内（因此，缩放似乎不会立即成为问题）。用于 MLE 的 n、x 和 y 值的样本有 10981 个样本。以下是一些定性显示样本的图表：

图 1. 噪声样本的直方图。

图 2. 噪声样本分别与 x 和 y 样本的散点图。

我的目标是计算k(i) 模型参数i=1,...,10 的最大似然估计值，以及它们的标准差kSE(i)（由渐近协方差矩阵输出的对角元素的平方根给出mlecov)。

对于最大似然估计，我最小化负对数似然：

我还为 MATLAB 提供了 mle 和 mlecov 使用的负对数似然 L(k(1),...,k(10)) 的分析梯度，这样梯度的数值近似值有望不会对我将要描述的数值问题产生影响。

数值问题

为了演示这个问题，我提出了三个场景。

场景 1。我直接在样本数据上运行mle 和mlecov。这将输出以下类似 Stata 的摘要：

-----------------------------------------------------------------------------
Coeffs   |      Val.     Std. Err.       z      P>|z|    [95% Conf. Interval]
---------+-------------------------------------------------------------------
   k1    |   -0.0153      0.0014     -11.27     0.000     -0.0179    -0.0126 
   k2    |    0.0075      0.0016       4.79     0.000      0.0045     0.0106 
   k3    |    0.0045      0.0006       7.44     0.000      0.0033     0.0056 
   k4    |    0.0131      0.0023       5.57     0.000      0.0085     0.0177 
   k5    |   -0.0101      0.0012      -8.45     0.000     -0.0125    -0.0078 
   k6    |    0.0114      0.0011      10.25     0.000      0.0092     0.0135 
   k7    |    0.0244      0.0011      21.86     0.000      0.0222     0.0266 
   k8    |   -0.0001      0.0004      -0.34     0.732     -0.0010     0.0007 
   k9    |   -0.0190      0.0018     -10.48     0.000     -0.0225    -0.0154 
   k10   |    0.0057      0.0009       6.32     0.000      0.0039     0.0074 
-----------------------------------------------------------------------------

“瓦尔”。列对应于k(i) 估计值和“Std. Err”。列对应于kSE(i)。 “P>|z|”列给出了原假设 k(i)==0 的单个系数 Wald 检验的 p 值（如果此 p 值为 <0.05，我们拒绝原假设并因此得出结论，系数 k(i) 在 95 % 水平）。

请注意，为了计算 k(i) 估计值的渐近协方差矩阵，mlecov 计算 L(k(1),...,k(10)) 的 Hessian H - 我提供了一个解析梯度。 H 的条件号为cond(H)=2.7437e3。 mlecov 函数对 Hessian 进行 Cholesky 因式分解，得到上三角矩阵 R 和 cond(R)=52.38。

场景 2。我将所有样本乘以0.1，从而在样本数据n*0.1、x*0.1 和y*0.1 上运行mle 和mlecov。这将输出以下摘要：

-----------------------------------------------------------------------------
Coeffs   |      Val.     Std. Err.       z      P>|z|    [95% Conf. Interval]
---------+-------------------------------------------------------------------
   k1    |   -0.0010      0.0001      -7.39     0.000     -0.0013    -0.0008 
   k2    |    0.0063      0.0016       3.97     0.000      0.0032     0.0093 
   k3    |    0.0494      0.0060       8.21     0.000      0.0376     0.0611 
   k4    |    0.0023      0.0024       0.95     0.340     -0.0024     0.0070 
   k5    |   -0.0462      0.0123      -3.75     0.000     -0.0704    -0.0221 
   k6    |    0.0014      0.0001      12.30     0.000      0.0012     0.0016 
   k7    |    0.0220      0.0011      20.86     0.000      0.0200     0.0241 
   k8    |    0.0078      0.0042       1.87     0.062     -0.0004     0.0160 
   k9    |   -0.0228      0.0020     -11.27     0.000     -0.0267    -0.0188 
   k10   |    0.0747      0.0097       7.70     0.000      0.0557     0.0937 
-----------------------------------------------------------------------------

p 值已更改。另外，现在cond(H)=9.3831e5 (!!!) 和cond(R)=968.6616。请注意，当我从均值和标准差模型中删除二阶项（x.^2 和 y.^2）时，不再存在此问题（即 p 值保持不变，k(i) 值除外，除了常数项k(1) 和k(6) 仅按0.1 缩放）。这是否表明存在数字问题？

场景 3。我还决定尝试将n、x 和y 缩放到区间[-1,1]，方法是将它们的样本除以最大元素（即n(i)=n(i)/max(abs(n))、x(i)=x(i)/max(abs(x)) 和y(i)=y(i)/max(abs(y))）。在此缩放样本上运行 mle 和 mlecov 会输出以下摘要：

-----------------------------------------------------------------------------
Coeffs   |      Val.     Std. Err.       z      P>|z|    [95% Conf. Interval]
---------+-------------------------------------------------------------------
   k1    |   -0.0347      0.0041      -8.40     0.000     -0.0428    -0.0266 
   k2    |    0.1193      0.0141       8.46     0.000      0.0917     0.1470 
   k3    |    0.0482      0.0164       2.94     0.003      0.0160     0.0803 
   k4    |   -0.0002      0.0120      -0.02     0.987     -0.0238     0.0234 
   k5    |   -0.0305      0.0103      -2.96     0.003     -0.0506    -0.0103 
   k6    |    0.0557      0.0035      16.11     0.000      0.0489     0.0624 
   k7    |    0.1131      0.0107      10.60     0.000      0.0922     0.1341 
   k8    |    0.1164      0.0128       9.13     0.000      0.0914     0.1414 
   k9    |   -0.1132      0.0094     -11.99     0.000     -0.1317    -0.0947 
   k10   |    0.0583      0.0079       7.37     0.000      0.0428     0.0738 
-----------------------------------------------------------------------------

p 值又变了！现在cond(H)=4.7550e3（高于场景 1（未缩放）但低于场景 2（一切乘以 0.1））。此外，cond(R)=68.9565，仅略高于场景 1。

我的问题

对我而言，这三个分析的预期行为是 k(i) 和 kSE(i) 会发生变化，但 p 值将保持不变 - 换句话说，缩放数据不应使任何模型系数更多或统计学意义较小。这与上述情况相反，每次 p 值都在变化！

请帮助我调试这个数字问题 - 或解释这是否实际上是预期的行为，我误解了一些东西。感谢您阅读这篇长文并提供帮助 - 我试图在此处封装所有相关问题的详细信息。

Answer 1

首先，我假设您正在控制抽样的随机种子，因此在所有情况下都是相同的。

已经解决了，我认为这可能与您尝试解决的优化问题有关。 我有第一手经验，当目标函数不是凸的时，微小的数值变化（在我的例子中，将对数似然函数缩放一个因子，或者等效地：添加所有数据点的副本）会改变你的结果。

我会尝试在所有参数中推导出对数似然函数的解析梯度。 这应该让您了解优化问题是否是凸的。 如果它不是凸的，则需要做一些事情来确保获得真正的 MLE。

• 优化函数 1000 次并选择对数似然最高的估计值
• 更改优化器的容差和步数
• 尝试其他优化器，例如信任区域搜索或粒子群

我会先模拟这个问题的一个更简单的版本，然后逐渐建立它，看看这个行为是从哪里开始发生的。例如，仅从 1 个参数表示均值，1 个参数表示噪声，然后看看 p 值会发生什么。