Maximum Likelihood Fit of Matrix: How can I find the Covariance Matrix?

Question

I have a 2d Data Matrix $X_{i,j}$ with a corresponding error matrix $E_{i,j}$. I now have a model $X_{i,j} = a_i + b_j + u_{i,j}$ with $u_{i,j}$ being distributed as $N(0, E_{i,j}^2 + \delta_j^2)$. I want to determine the parameters $a_i, b_j$ and $\delta_j$. Currently I am minimizing the following log likelihood function: $$\text{logL}_{i,j} = 1/2 \text{log}(2\pi(E_{i,j}^2 + \delta_j^2) - 1/2(X_{i,j} - a_i - b_j)^2 / (E_{i,j}^2 + \delta_j^2)$$ with logL = -$\sum\limits_{i,j}\text{logL}_{i,j}$. I am minimizing with the constraint that $\sum\limits_j b_j= 0$ (edit: fixed expression). Therefore I can cannot use the scipy methods which would give me directly an estimation for the inverse hessian (edit: I am currently using SLSQP). The model works very well but I have no access to the fit uncertainties. My question is: How can I determine the covariance matrix in my problem. edit:correctet equation typo

Answer 1

This is more of an extended comment rather than an answer.

First, I wonder if you really want just a common $a$ rather than $a_1, a_2,\ldots$. Otherwise you have a case where the number of parameters increases with sample size which can make obtaining good estimators more difficult.

If the model is $X_{i,j}=a+b_j+u_{i,j}$ with the $u_{i,j}$'s being independent and $u_{i,j}\sim N(0,E_{i,j}^2+\delta_j^2)$, then the log likelihood is

$$\log(L)=-\sum _{i=1}^I \sum _{j=1}^J \left(\frac{(X_{i,j}-a-b_j)^2}{E_{i,j}^2+\delta_j^2}+\log \left(2 \pi \left(E_{i,j}^2+\delta_j^2\right)\right)\right)/2$$

The constraint $\sum_{j=1}^{J} b_j=0$ can be implemented by letting $b_1=-\sum_{j=2}^{J} b_j$.

Rather than listing all of the second order partial derivatives, it's much easier to use a computer algebra system (CAS). The example below uses Mathematica.

(* Generate some data *)
nI = 60;
nJ = 2;
ee = ConstantArray[1, {nI, nJ}];
aa = 0;
bb = ConstantArray[0, nJ];
δδ = ConstantArray[2, nJ];
SeedRandom[12345]; x = 
 Table[aa + bb[[j]] + 
   RandomVariate[
     NormalDistribution[
      0, (ee[[i, j]]^2 + δδ[[j]]^2)^(1/2)], 1][[1]], {i,
    nI}, {j, nJ}];

(* Incorporate the restriction and create the log of the likelihood *)


b[1] = -Sum[b[j], {j, 2, nJ}]
logL = Sum[LogLikelihood[NormalDistribution[0, (ee[[i, j]]^2 + δ[j]^2)^(1/2)], 
  {x[[i, j]] - a - b[j]}], {i, nI}, {j, nJ}];

(* Maximum likelihood estimates *)
mle = FindMaximum[{logL, And  @@ Table[δ[j] > 0, {j, nJ}]}, 
  Join[{{a, 0}}, Table[{b[j], 0}, {j, 2, nJ}], Table[{δ[j], 1}, {j, nJ}]]]
(* {-268.3853361029104`,{a -> -0.7177507444852211`,b[2] -> 0.31930855913062345`,
     δ[1] -> 2.0864419431898744`,δ[2] -> 1.9791700109759636`}} *)

(* Approximate parameter covariance matrix *)
(cov = -Inverse[(D[logL, {Join[{a}, Table[b[j], {j, 2, nJ}], Table[δ[j], {j, 1, nJ}]], 2}]) /. mle[[2]]]) // MatrixForm

$$\left( \begin{array}{cccc} 0.0427931 & -0.00181719 & \text{4.3471562574947714$\grave{ }$*${}^{\wedge}$-18} & 0. \\ -0.00181719 & 0.0427931 & -\text{4.347156257494772$\grave{ }$*${}^{\wedge}$-18} & 0. \\ \text{4.347156257494772$\grave{ }$*${}^{\wedge}$-18} & -\text{4.347156257494772$\grave{ }$*${}^{\wedge}$-18} & 0.0548579 & 0. \\ 0. & 0. & 0. & 0.0514367 \\ \end{array} \right)$$

(* Approximate standard errors for parameters *)
Thread[Join[{aSE}, Table[bSE[j], {j, 2, nJ}], Table[δSE[j], {j, nJ}]] -> Diagonal[cov]^0.5]
(* {aSE -> 0.206865031629407`,bSE[2] -> 0.206865031629407`, δSE[1] -> 0.23421773987043146`,
    δSE[2] -> 0. 22679660389187373`} *)

A bootstrap as suggested by @MarkLStone is a good choice that generally gives good results with a wider variety of distributions but does require a bit more programming and more details to implement than what is given in the referenced answer.