Maximum likelihood using simulated annealing

Question

I have a sequence of random variable $Y_i=X_i\eta_{i}$ where $X_i \sim \mathcal{P}(\lambda)$ and are independent, $\eta _i \sim \mathcal{N}(0,1)$ also, and a set of values $(y_1,...,y_n)$. I would like to find the best parameter $\lambda$ that maximizes the likelihood function defined as :
$\mathcal{L}_{\lambda}(y_1,...y_n)=\prod\limits_{i=1}^{n}\mathcal{L}_{\lambda}(y_i)=\prod\limits_{i=1}^{n}\int\mathcal{L}_{\lambda}(y_i /x_i)\mathcal{L}_{\lambda}(x_i)dx_i=\prod\limits_{i=1}^{n}\int\frac{e^\frac{-y_i^{2}}{2x_i^{2}}}{\sqrt{2\pi}x_i}\lambda e^{-\lambda x_i}dx_i=\prod\limits_{i=1}^{n}\mathbb{E}(\frac{e^\frac{-y_i^{2}}{2X_i^{2}}}{\sqrt{2\pi}X_i})=\mathbb{E}(\prod\limits_{i=1}^{n}\frac{e^\frac{-y_i^{2}}{2X_i^{2}}}{\sqrt{2\pi}X_i})$
The first problem that i am encountering is that when i plot the expected value above using the approximation $\sum\frac{ \frac{e^\frac{-y_i^{2}}{2X_i^{2}}}{\sqrt{2\pi}X_i}}{N}$ with respect to different values of $\lambda $, i find that it is a decreasing function and that a value of $\lambda$ close to $0$ is a solution to my problem. This seems strange because i am asked to find the best parameter $\lambda$ using the simulated annealing algorithm. Also, the calculations i am doing to evaluate the likelihood function for a certain $(y_1,...,y_n)$ are really heavy, which makes the simulated annealing algorithm obsolete. Can anyone help me? thanks

Answer 1

Here is an approach using Mathematica.

First construct the density function for $Y_i$:

PDF[TransformedDistribution[x z, {x \[Distributed] ExponentialDistribution[\[Lambda]], 
   z \[Distributed] NormalDistribution[0, 1]}], y]

$$\begin{array}{cc} \{ & \begin{array}{cc} \frac{\lambda G_{0,3}^{3,0}\left(\frac{y^2 \lambda ^2}{8}| \begin{array}{c} 0,0,\frac{1}{2} \\ \end{array} \right)}{2 \sqrt{2} \pi } & y\neq 0 \\ \frac{\lambda G_{0,3}^{3,0}\left(\frac{y^2 \lambda ^2}{8}| \begin{array}{c} 0,0,\frac{1}{2} \\ \end{array} \right)}{\sqrt{2} \pi } & y=0 \\ \end{array} \\ \end{array}$$

Now in the following we assume that no value of $Y_i$ is exactly zero:

(* Generate some data *)
SeedRandom[12345];
\[Lambda]0 = 3;
n = 100;
x = RandomVariate[ExponentialDistribution[\[Lambda]0], n];
z = RandomVariate[NormalDistribution[], n];
y = x z;

(* Log likelihood *)
logL = n Log[\[Lambda]] - (3 n/2) Log[2] - n Log[\[Pi]] + 
   Sum[Log[MeijerG[{{}, {}}, {{0, 0, 1/2}, {}}, (
      y[[i]]^2 \[Lambda]^2)/8]], {i, n}];

(* Method of moments estimate *)
\[Lambda]Initial = Sqrt[2]/StandardDeviation[y]
(* 3.159259259137787 *)

(* Maximum likelihood estimate *)
mle = FindMaximum[{logL, \[Lambda] > 0}, {\[Lambda], \[Lambda]Initial}]
(* {33.75465530200371,{\[Lambda] -> 2.7946001175877813}} *)
(* Alternative (but slower) *)
(* mle=NMaximize[{logL,\[Lambda]>0},\[Lambda],Method -> "SimulatedAnnealing"] *)

(* Check to see if first partial derivative is approximately zero *)
dlogL\[Lambda] = D[logL, \[Lambda]] /. mle[[2]]
(* 1.000782830828939 * 10^(-6) *)

(* Estimate of standard error *)
se = Sqrt[-1/(D[logL, {\[Lambda], 2}]) /. mle[[2]]]
(* 0.39489891632730545 *)

(* Another check:  plot the log likelihood *)
Plot[logL, {\[Lambda], 1, 4}, Frame -> True, FrameLabel -> {"\[Lambda]", "Log(likelihood)"}]