How to compute MLE for Pareto distribution data with right censored observation by using R programming?

Question

I am doing survival analysis and writing codes to compute MLE for several distributions. Yet, I get stuck while writing for Pareto distribution with right censored observation.

For complete/uncensored data, it can be dealt with by using the following coding provided by Macro in this post: How do I fit a set of data to a Pareto distribution in R?

    pareto.MLE <- function(x)
    {
      n <- length(x)
      m <- min(x)
      a <- n/sum(log(x)-log(m))
      return( c(m,a) ) 
    }

    # example. 
    library(VGAM)
    set.seed(1)
    z = rpareto(1000, 2, 6) 
    pareto.MLE(z)
    [1] 1.000014 5.065213

So, how about right censored data? I tried to follow the steps the previous post shows for right censored data by using its survival function where S(t)=(b/t)^a, but I still couldn't manage to get a solution, as I got n*log(m)-sum(log(t))=0, there's no more a after the differentiation w.r.t. a.

What should I do?

Answer 1

For a random sample $x_1, x_2,\ldots,x_n$ the likelihood with parameters $m$ and $a$ with the censoring done at $x_c$ is given by

$$\left(\left(\frac{m}{x_c}\right)^{\alpha }\right)^{n_c} \prod _{i=1}^{n-n_c} a m^a x_i^{-a-1}$$

where $n_c$ is the number of censored observations. The maximum likelihood estimator of $m$ is still $\min x_i$ and the maximum likelihood estimator of $a$ becomes

$$\frac{n-n_c}{\sum _i^{n} \log (x_i)-n \log (\hat{m})}$$

Some R code follows:

# Function to create a random sample from a right-censored Pareto distribution
r.censored.pareto <- function(n, m, a, xc) {
  # n is sample size
  # m and a are the Pareto parameters
  # xc is the censoring value
  x <- m*(1 - runif(n))^(-1/a)
  x[x >= xc] <- xc
  x
}

# Function to return maximum likelihood estimates of m and a
censored.pareto.MLE <- function(x, xc) {
  # x is the random sample from the censored Pareto distribution
  # xc is the censoring value

  # Count the number of censored observations
  nc <- length(x[x >= xc])
  m.mle <- min(x)
  a.mle <- (length(x) - nc)/(sum(log(x)) - length(x)*log(m.mle))
  return(c(m.mle, a.mle))
}

# Generate random sample from a right-censored Pareto distribution
  set.seed(12345)
  x <- r.censored.pareto(1000, 5, 5, 8)

# Find maximum likelihood estimates of m and a
  censored.pareto.MLE(x, 8)
# [1] 5.001137 4.864037