5
$\begingroup$

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

enter image description here

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Stelios Mar 22 '15 at 18:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.