I have a problem with four-parameter estimation. I have derived the variances for the estimated parameters using Monte Carlo simulations (numerical ones) and theoretical ones using the inverse of the Fisher information matrix.
The theoretical and the numerical ones have differences. The numerical ones can have even smaller variances for some parameters than the theoretical ones. I suspect the reason is the biased nature of the estimator. I wanted to compute the biased CRBs so that they could make a fair comparison. On top of that, I need to show some parameter sweeps with these variances.
The biased CRBs can be found by the following formula.
$$ \text{CRB}(\theta) = \text{diag}\{ \left[\mathbf{1} + \nabla_\theta \mathbb{B}(\theta) \right] \mathbf{I}^{-1}(\theta) \left[\mathbf{1} + \nabla_\theta \mathbb{B}(\theta) \right]^\text{T} \}, $$
where $\mathbf{1}$ is the identity matrix, $\mathbf{I}$ is the Fisher information matrix. $\mathbb{B}(\theta)$ is the bias vector, but $\nabla_\theta \mathbb{B}(\theta)$ is the bias-gradient matrix.
The bias can come from factors such as insufficient sample size in the integrals involved in the log-likelihood ($N$), constrained parameter space, optimizer errors, etc. For my problem, I think I have all of them, especially insufficient $N$ and constrained space. A functional form of $\nabla_\theta \mathbb{B}(\theta)$ is almost impossible. However, I can see that there are some Bayesian inference techniques to compute the biased-CRBs, such as in [1]. It goes beyond my understanding and is a bit complicated.
I was thinking that I could probably use the numerical biases to get the bias gradients with Monte Carlo simulations. However, it is computationally so expensive. As the gradients are with respect to the parameters and I want to do a few parameter sweeps, I need to perform Monte Carlo simulations for $M^4 \times M_c$ computations of the optimizer where $M$ is the number of values in one parameter sweep, and $M_c$ is the number of Monte Carlo runs. It also takes some time to prepare the simulated measurements for each case. I code with MATLAB, and I can only parallelize one of these sweeps.
Are there smarter ways to obtain the biased-CRBs conceptually or with programming tactics? I would appreciate any sort of help.
[1] C. L. Matson and A. Haji, “Biased Cramér-Rao lower bound calculations for inequality-constrained estimators,” J. Opt. Soc. Am. A, vol. 23, no. 11, pp. 2702–2713, 2006, doi: 10.1364/josaa.23.002702.
