# Numerical MLE for Rayleigh distribution

I am given a rayleigh distribution described by, =$$f\left(x|\theta\right)\:=\:\frac{x}{\theta ^2}e^{-\left(\frac{x^2}{2\theta ^2}\right)}$$

I need to find a numerical estimate of the MLE of $$\theta^2$$ using the Newton-Raphson method. I have gone through a found the derivatives of the log likelihood function as follows: $$\frac{\partial }{\partial \theta}l = 1/\theta^3 \sum_{i=1}^n (x_i^2) -2n/\theta$$ and $$\frac{\partial^2 }{\partial \theta^2}l = -3/\theta^4 \sum_{i=1}^n (x_i^2) +2n/\theta^2$$

I have attempted to plug these into the Newton-Raphson algorithm but have been unable to come up with a reasonable solution. So I am just hoping for some clarification as to whether my errors have come from incorrect derivatives of the log likelihood function or from incorrect application of the Newton-Raphson method.

– Dave
Sep 10, 2021 at 0:11

I think your derivatives look okay. We don't have a basis to judge your implementation though. What steps have you taken to check it?

However, if we're trying to estimate $$\theta^2$$, why not do the derivatives directly in terms of that, and maybe save a step, since you're then iterating directly to an estimate of the desired quantity.

I think it goes like this, but you might want to double check my work:

$${l(\theta^2)} = \sum_i \left[\log(x_i) - \log(\theta^2) - \frac{1}{\theta^2} x_i^2/2\right]$$

$$\frac{\partial l}{\partial \theta^2}= -n/\theta^2 + \frac{1}{(\theta^2)^2}\sum_i x_i^2/2$$

$$\frac{\partial^{2} l}{\partial (\theta^2)^2} = +n/(\theta^2)^2 - \frac{1}{(\theta^2)^3}\sum_i x_i^2$$

• Thanks for that, you were right in that my implementation fell apart as I was taking derivatives with respect to theta rather than theta squared. Sep 6, 2021 at 5:28
• That should still get you to an estimate of $\theta$, and you can then get $\widehat{\theta^2}$ by squaring $\hat\theta$. It's just an extra step (well I guess it's also more involved getting an estimate of the standard error as well). Sep 6, 2021 at 5:33

The density is

$$f(x \vert \theta) = \dfrac{x}{s} \exp\left( -\dfrac{x^2}{2s}\right)$$

Here $$s = \theta^2$$. We will take derivatives with respect to $$s$$ for simplicity. The first derivative of the log likelihood is

$$\dfrac{\partial \ell}{\partial s} = \sum_i \dfrac{2s - x_i^2}{2s^2}$$

and the second derivative is

$$\dfrac{\partial ^2 \ell }{\partial s^2} = \sum_i \dfrac{x^2-s}{s^3}$$

These are fairly straight forward to code in python

import numpy as np

def dll(s, x):

return -np.divide(np.power(x, 2) - 2*s, 2*np.power(s, 2))

def ddll(s, x):

return -np.divide(s-np.power(x, 2), np.power(s, 3))



Let's generate some data

from scipy.stats import rayleigh
sample = rayleigh(loc=0,scale = 2.5).rvs(10_000)



and implement newton's method for 1000 iterations

s = 1.0

for i in range(1000):

s -= dll(s, sample).sum()/ddll(s, sample).sum()


Results will depend on the random seed, but my sample gives me

s
>>>6.177422...


Let's compare to the maximum likelihood estimate

1/(2*sample.size)*np.power(sample,2).sum()
>>>6.177422...


Our implementation is within acceptable precision of the MLE obtained via algebra.

We can further compare with scipy's implementation of newton's method

from scipy.optimize import newton
newton(lambda s: dll(s, sample).sum(), x0=1.0)
>>> 6.177422...

• Hi, would you be able to clarify how you got the first and second derivatives of the log-likelihood function into the form you have? Sep 6, 2021 at 5:49
• @helpneeded I used computer algebra to get the derivatives. They are probably not exactly what you would get were you to apply the chain/product/quotient rules for differentiation. Remember that $\ell(s;x) = \sum \log(f(x\vert s))$. Just be patient and take derivatives of that. Sep 6, 2021 at 13:55