# Numerical issues in Maximum Likelihood Estimation

I am estimating the following noise model:

\begin{equation} n \sim \mathcal N\left(\mu(x,y),\sigma(x,y)^2\right)\in\mathbb R \end{equation}

\begin{equation} \begin{cases} \mu(x,y) = k_1+k_2x+k_3x^2+k_4y+k_5y^2 \\ \sigma(x,y) = k_6+k_7x+k_8x^2+k_9y+k_{10}y^2 \\ \end{cases} \end{equation}

where $x\in[0,3]$ and $y\in[0,\pi/2]$ (thus, scaling does not immediately seem to be an issue). The sample size $\{\hat n_i, \hat x_i, \hat y_i\}_{i=1}^N$ has size $N=10981$. Here are some graphs to show the sample qualitatively: Figure 1. Histogram of the noise samples $\{\hat n_i\}_{i=1}^N$. Figure 2. Scatter plot of the noise samples $\{\hat n_i\}_{i=1}^N$ vs. the $\{\hat x_i\}_{i=1}^N$ and $\{\hat y_i\}_{i=1}^N$ samples respectively.

I use MATLAB's mle and mlecov functions in order to extract the maximum likelihood estimates $\tilde k_j$ of the $k_j$, $j=1,...,10$, model parameters as well as their standard deviation, $\tilde k_{j,\text{SE}}$ (given by the square root of the diagonal elements of the asymptotic covariance matrix output by mlecov).

For the maximum likelihood estimation, I minimize the negative log likelihood $L(k_1,...,k_{10})$:

\begin{equation} \tilde k_{j=1,\ldots,10} = \mathrm{arg}\min_{k_1,\ldots,k_{10}} L(k_1,...,k_{10}) := \mathrm{arg}\min_{k_1,\ldots,k_{10}}\left[ -\ln\left(\prod_{i=1}^N \mathrm{Pr}(\hat n_i | \hat x_i, \hat y_i, k_1,\ldots,k_{10}) \right) \right] \end{equation}

where

\begin{equation} \mathrm{Pr}(\hat n_{i} | \hat x_i, \hat y_i, k_1,\ldots,k_{10}) = \frac{1}{\sqrt{2\pi\sigma(\hat x_i,\hat y_i)^2}}\exp\left( -\frac{(\hat n_{i}-\mu(\hat x_i,\hat y_i))^2}{2\sigma(\hat x_i,\hat y_i)^2} \right). \end{equation}

I also give MATLAB the analytical gradient of $L(k_1,...,k_{10})$, used by mle and mlecov such that numerical approximations of the gradient hopefully do not contribute to the numerical issue I am about to describe.

Numerical Issue

To demonstrate the issue, I present three scenarios.

Scenario 1. I directly run mle and mlecov on the sample data $\{\hat n_i, \hat x_i, \hat y_i\}_{i=1}^N$. This outputs the following Stata-like summary:

-----------------------------------------------------------------------------
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
-----------------------------------------------------------------------------


The "Val." column corresponds to $\tilde k_i$ and the "Std. Err." column corresponds to $\tilde k_{j,\text{SE}}$. The "P>|z|" column gives the p-value for a single coefficient Wald test of the null hypothesis $\tilde k_i=0$ (if this p-value is $<0.05$, we reject the null hypothesis and thus conclude that the coefficient $k_i$ may be significant at the 95% level).

Note that to compute the asymptotic covariance matrix of the $\tilde k_i$ estimates, mlecov computes the Hessian $H$ of $L(k_1,...,k_{10})$ - which I provide an analytic gradient for. The condition number of $H$ is $\text{cond}(H)=2.7437\cdot 10^3$. The mlecov function does a Cholesky factorization of the Hessian, which gives the upper-triangular matrix $R$ with $\text{cond}(R)=52.38$.

Scenario 2. I multiply all samples by $0.1$ and thus run mle and mlecov on the sample data $\{\hat n_i\cdot 0.1, \hat x_i\cdot 0.1, \hat y_i\cdot 0.1\}_{i=1}^N$. This outputs the following summary:

-----------------------------------------------------------------------------
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
-----------------------------------------------------------------------------


The p-values have changed. Also, now $\text{cond}(H)=9.3831\cdot 10^5$ (!!!) and $\text{cond}(R)=968.6616$. Note that when I remove the second order terms ($x^2$ and $y^2$) from the $\mu$ and $\sigma$ models, there is no longer this problem (i.e. the p-values stay the same and the $\widetilde k_i$ values, except for the constant terms $k_1$ and $k_6$, are simply scaled by $0.1$). Does this indicate a numerical issue?

Scenario 3. I decided to also try scaling $n$, $x$ and $y$ to the interval $[-1,1]$ by dividing the samples by their largest element (i.e. $\hat n_i\leftarrow \frac{\hat n_i}{\max_{j=1,...,N}(|\hat n_j|)}$, $\hat x_i\leftarrow \frac{\hat x_i}{\max_{j=1,...,N}(|\hat x_j|)}$, $\hat y_i\leftarrow \frac{\hat y_i}{\max_{j=1,...,N}(|\hat y_j|)}$). Running mle and mlecov on this scaled sample outputs the following summary:

-----------------------------------------------------------------------------
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
-----------------------------------------------------------------------------


The p-values have changed again! Now $\text{cond}(H)=4.7550\cdot 10^3$ (higher than Scenario 1 (unscaled) but lower than Scenario 2 (everything multiplied by 0.1)). Also, $\text{cond}(R)=68.9565$, which is only slightly higher than for Scenario 1.

My problem

The expected behavior across the three analyses, for me, is that $\tilde k_i$ and $\tilde k_{i,SE}$ would change but the p-values would remain the same - in other words, scaling the data should not make any model coefficient more or less statistically significant. This is contrary to the above scenarios, where the p-values change each time!

Please help me to debug this numerical issue - or explain whether this is in fact the expected behavior and I have misunderstood something. Thank you for reading this long post and helping - I tried to encapsulate all relevant problem details here.

• I haven't used MATLAB but does it allow for increasing the precision of the calculations? Is the data available? – JimB Apr 17 '18 at 3:41