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

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?

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