0
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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