MLE of $f(x\mid\theta) = \theta x^{\theta−1}e^{−x^{\theta}}I_{(0,\infty)}(x)$

Question

$X_1$, $X_2$, and $X_3$ are i.i.d random variables having pdf

$$f(x\mid\theta) = \theta x^{\theta−1}e^{−x^{\theta}}I_{(0,\infty)}(x)$$ where $\theta \gt0$. Using the observed values of these random variables given below, give the mle of $\theta$.

$x_1= 0.60$, $x_2= 5.17$, $x_3= 0.23$

My Attempt:

We have

$$\begin{align*} L(\theta\mid x)=\prod_{i=1}^3 \theta x_i^{\theta−1}e^{−x_i^{\theta}} &=\theta^3 \left(\prod_{i=1}^3 x_i\right)^{\theta-1} \mathsf{exp}\left(-\left(\sum x_i^{\theta}\right)\right) \end{align*}$$

Upon taking the taking the natural log of both sides, we get

$$\mathcal{L}(\theta\mid x)=3\cdot\text{log}(\theta)+(\theta-1)\left(\sum \text{log } x_i\right)-\sum\left(x_i^{\theta}\right)$$

Then we have

$$0=\frac{d}{d\theta}\left(\mathcal{L}(\theta\mid x)\right)=\frac{3}{\theta}+\sum\text{log }x_i-\sum\left(\text{log }x_i\cdot x_i^{\theta}\right)$$

Is this valid thus far? If so, how can I proceed from here to find which value of $\theta$ makes this function $0$? It's not immediately clear to me how to do so.

Edit:

Using Newton's Method to obtain a mle we have

$$\hat{\theta}_{j+1} = \hat{\theta}_j - \frac{g(\hat{\theta}_j)}{g'(\hat{\theta}_j)}$$

where $g(\theta)=\frac{d}{d\theta}\left(\mathcal{L}(\theta\mid x)\right)$.

With an initial value of $\theta=1$, I get that this converges to approximately $0.67715$ after only a few iterations.

Check:

$$\frac{3}{0.67715}+\sum\text{log }x_i-\sum\left(\text{log }x_i\cdot x_i^{0.67715}\right)\approx 0$$

Is it necessary to show that we have a global maximum still?

Answer 1

If you're not sure whether or not your answer is correct, a useful check is to plot a graph of the log-likelihood function and see if your purported MLE looks visually to give the maximising value. I will do this below, but I include the mathematics for deriving the MLE in the general case.

---

MLE in the general case: For IID data from this distribution, you have log-likelihood:

$$\ell_\mathbf{x}(\theta) = n \ln \theta + (\theta-1) \sum_{i=1}^n \ln x_i - \sum_{i=1}^n x_i^\theta \quad \quad \text{for } \theta>0.$$

The corresponding score function is:

$$s_\mathbf{x}(\theta) = \frac{d\ell_\mathbf{x}}{d\theta}(\theta) = \frac{n}{\theta} + \sum_{i=1}^n (1-x_i^\theta) \ln x_i,$$

and the observed information is:

$$I_\mathbf{x}(\theta) = - \frac{d^2\ell_\mathbf{x}}{d\theta^2}(\theta) = \frac{n}{\theta^2} + \sum_{i=1}^n x_i^\theta (\ln x_i)^2 > 0.$$

We can see from the positive observed information that the log-likelihood is strictly concave, which means that the MLE will be at the unique critical point (unless the score is monotone, in which case the maximum is approached at the boundary of the parameter range, and there is no MLE). The critical point is given implicitly by solving the score equation:

$$0 = s_\mathbf{x}(\hat{\theta}) = \frac{n}{\hat{\theta}} + \sum_{i=1}^n (1 - x_i^\hat{\theta}) \ln x_i.$$

There is no closed-form expression for the MLE in this case, so we need to find it numerically.

---

Iterative algorithm for MLE: Applying Newton's method, with your chosen starting-point, gives:

$$\hat{\theta}_0 = 1 \quad \quad \quad \hat{\theta}_{k+1} = \hat{\theta}_{k} + \frac{s_\mathbf{x}(\hat{\theta})}{I_\mathbf{x}(\hat{\theta})} = \hat{\theta}_{k} + \frac{n \hat{\theta}_k + \hat{\theta}_k^2 \sum_{i=1}^n (1 - x_i^{\hat{\theta}_k}) \ln x_i}{n + \hat{\theta}_k^2 \sum_{i=1}^n x_i^{\hat{\theta}_k} (\ln x_i)^2}.$$

(Note: The starting point you have chosen is a reasonable one. With some calculus, it is possible to show that $\mathbb{E}(X) = \Gamma(1 + 1/\theta)$, so we could approximate $\bar{x} \approx \Gamma(1 + 1/\theta)$ as a starting point for the iteration. However, the problem is that this already requires numerical solution, so it is not a great starting point. The value you have chosen is reasonable, and the iteration should converge quite rapidly in any case.) We can implement this iteration algorithm in the following R code:

#Create function to find the MLE via iteration
#The input m is the number of iterations to perform (default is five iterations)
MLE_ITERATION <- function(x, m = 5) {
    n  <- length(x);
    T  <- 1;
    theta <- rep(T, m+1);
    for (k in 1:m) {
        NUMERATOR   <- n*theta[k] + theta[k]^2 * sum((1-x^theta[k])*log(x));
        DENOMINATOR <- n + theta[k]^2 * sum(x^theta[k]*(log(x))^2);
        theta[k+1]  <- theta[k] + NUMERATOR/DENOMINATOR; }
    theta; }    

---

Application to your data set: You have the data vector $\mathbf{x} = (0.60, 5.17, 0.23)$. With $m=10$ iterations (which is more than you need) you get the MLE $\hat{\theta} = 0.6771516$. The log-likelihood and MLE are shown here:

Here is the R code used to generate the MLE and the plot:

#Enter your data
x <- c(0.60, 5.17, 0.23);

#Choose number of iterations
m <- 10;

#Generate the iterations, and display the last value
THETA_ITER <- MLE_ITERATION(x, m);
THETA_ITER[m+1];

[1] 0.6771516

#Generate vectorised log-likelihood function
LOG_LIKE <- function(theta) {
    LL <- rep(0, length(theta)); 
    for (i in 1:length(theta)) {
        LL[i] <- length(x)*log(theta[i]) + (theta[i]-1)*sum(log(x)) - sum(x^theta[i]); }
    LL }

DATA <- data.frame(Theta    = 1:200/100,
                   Log_Like = LOG_LIKE(1:200/100)); 

#Plot the log-likelihood function with MLE
library(ggplot2);
ggplot(data = DATA, aes(x = Theta, y = Log_Like)) + 
    geom_line(size = 1.2) +
    geom_vline(xintercept = THETA_ITER[m+1],
               size = 1.2, linetype = 'dashed', colour = 'red') +
    theme(plot.title    = element_text(hjust = 0.5, face = 'bold'),
          plot.subtitle = element_text(hjust = 0.5)) +
    ggtitle('Plot of Log-Likelihood Function') +
    labs(subtitle = '(Red line shows MLE - obtained via iteration)') +
    xlab(expression(theta)) + ylab('Log-Likelihood');