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.

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

