2
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ If your question is answered, please consider upvoting the helpful responses and accepting one of them, not erasing your question to thank people for answering (then future readers cannot benefit). Upvoting and accepting answers is how we express gratitude here. $\endgroup$
    – Dave
    Sep 10, 2021 at 0:11

2 Answers 2

3
$\begingroup$

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$$

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – helpneeded
    Sep 6, 2021 at 5:28
  • $\begingroup$ 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). $\endgroup$
    – Glen_b
    Sep 6, 2021 at 5:33
3
$\begingroup$

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...
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – helpneeded
    Sep 6, 2021 at 5:49
  • $\begingroup$ @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. $\endgroup$ Sep 6, 2021 at 13:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.