# RandomVariate, MSE and Cramer-Rao bound calculation

I'm using the following code to plot the variance of a moment estimate of a parameter of Rice distribution and its' Cramer-Rao bound.

1. First I generate $10^3$ vectors $\texttt{R}$ ($10^3$ elements in each) with Rice distribution with a fixed parameter $\alpha=10$ for a given parameter $\beta$. Then from this sample I calculate the bias and the variance of estmates: $\texttt{EstimateBias}$, $\texttt{EstimateVariance}$.
2. Then I form a table for different parameters $\beta$ (in log scale $10^{-\frac{\beta}{10}}$).
3. Then I evaluate Cramer-Rao bound as the inverse of the Fisher information $\texttt{FisherInformation}$. For this I take the second derivative of the probability density, invert and substitute the mean for the sample $x\to \texttt{Mean}[\texttt{RiceDistribution}[\alpha ,\beta ]]$ (i.e. taking expectation).
4. Then I plot them in log scale adding to the variance of the evaluated estimator the square of its bias (to get the MSE). Black - Cramer-Rao bound devided by the amount of samples, red - MSE of the estimator.

And what I get is actually wrong since the MSE is below the Cramer-Rao lower bound.

Is there any problem with the code or am I wrong with understanding?

RandomGenerator[α_,β_,Nstat1_,Nstat2_]:=RandomVariate[RiceDistribution[α,β],{Nstat1,Nstat2}];
Program[α_,β_,Nstat1_,Nstat2_]:=
Block[{},
R=RandomGenerator[α,10^(-(β/10)),Nstat1,Nstat2];
Estimates=ParallelTable[FindDistributionParameters[R[[i,All]],RiceDistribution[x,10^(-(β/10))],ParameterEstimator->"MethodOfMoments"],{i,1,Nstat2,1}][[All,1,2]];
EstimateBias=α-Mean[Estimates];
EstimateVariance=Variance[Estimates];
];

MoMEstimatesTable=Table[Progr[10,β,10^3,10^3];{EstimateBias,EstimateVariance},{β,-10,10,1}];

FisherInformation[α_,β_]:=
Evaluate[
-(FullSimplify[D[PDF[RiceDistribution[α,β],x],α,2}],
Assumptions->{x>0}]/.x->Mean[RiceDistribution[α,β]])^-1
];

FisherInformationTable=ParallelTable[FisherInformation[10,10^(-β/10)],{β,-10.,10,1}];

ListPlot[
{10 Log10[FisherInformationTable/1000],
10 Log10[MoMEstimatesTable[[All, 2]]+ MoMEstimatesTable[[All, 1]]^2]},
DataRange -> {-10., 10}, Joined -> True, Frame -> True, GridLines -> Automatic,
PlotStyle -> {{Black, Thick}, {Red, Thick}}
]


• Are you sure about the CRB formula? I believe it equals $\left(\mathbb{E}\left(\frac{\partial^2 \ln f(x;a)}{\partial a^2}\right)\right)^{-1}$. Also, using the rule x-> Mean[...] on an expression (function) does not provide the expectation of the function, rather the function evaluated at the mean of the random variable, which is not the same in general. – Stelios Mar 22 '15 at 18:17