5
$\begingroup$

I need to implement an algorithm in R that finds the MLE for $\lambda$ and $k$ (the truncation level: only values larger than k are allowed) for a truncated Poisson distribution. I seem to not be able to find this anywhere on google, and I'm struggling to derive the MLE in the first place.

$\endgroup$
3
  • $\begingroup$ Sorry to be thick but $k$ (I assume to be the truncation level) is unknown? $\endgroup$
    – JimB
    Oct 16, 2015 at 23:45
  • $\begingroup$ yes the truncation level is unknown $\endgroup$
    – Kathi
    Oct 16, 2015 at 23:53
  • $\begingroup$ What is your starting point to derive left-truncation log likelihood? Could you please clarify? $\endgroup$
    – user311529
    Feb 16, 2021 at 12:12

1 Answer 1

6
$\begingroup$

Update: Now that the question is more specific, I've added the log likelihood and R code for left-trunction.

Right-truncation. For $Y_i$ distributed as a truncated Poisson with parameters $\lambda$ and truncation parameter $k$, the log likelihood of $n$ random samples is given by

$$\log(\lambda) \sum_{i=1}^n y_i -n \lambda -n \log \left({\sum_{j=0}^k {{e^{-\lambda}\lambda^j}\over{j!}}}\right) $$

This is maximized for any value of $\lambda$ when $\hat{k}=\max(y_1, y_2,\ldots, y_n)$ (which is the smallest value that $k$ can take). The maximum likelihood estimate of $\lambda$ requires an iterative procedure. Moore (1952, Biometrika) provides a good and easily calculated initial estimate of $\hat{\lambda}$:

$$\hat{\lambda}_0 ={{\sum_{i=1}^n y_i}\over{\sum_{i=1}^n I(y_i<\hat{k})}}$$

where $I(\cdot)$ is the indicator function taking on 1 if $y_i<\hat{k}$ and 0 otherwise. (So it's a mean adjusted upwards a bit.)

Some not so elegant R code follows:

# MLE estimation of a truncated Poisson with unknown truncation level

# Objective is to find estimates of lambda (underlying Poisson mean) and
# k (unknown trunction value).

# Generate some samples from a known truncated Poisson distribution
  lambda <- 10  # Poisson mean
  k <- 8        # Truncation level
  n <- 100      # Sample size
  y <- rpois(n*4,lambda) # Oversample to be more certain of getting n samples
  # Keep just the first n legitimate observations
  y <- y[y <= k][1:n]

# MLE for k and initial estimate for lambda
# From Moore (1952)  Biometrika
# "The estimation of the Poisson parameter from a truncated distribution"
  khat <- max(y)
  lambda0 <- sum(y)/length(y[y < khat])

# Define log of the likelihood function
  logL <- function(lambda,ysum,n,k) {
    log(lambda)*ysum - n*lambda - n*ppois(k,lambda,log.p=TRUE)
  }

# Find maximum likelihood estimate of lambda  
  trPoisson <- optim(lambda0, logL, ysum=sum(y), n=length(y), k=khat, 
    method="BFGS", control=list(fnscale=-1))

# Show results
  cat("MLE of truncation value =",khat,"\n")
  # MLE of truncation value = 8 
  cat("MLE of lambda =",trPoisson$par,"\n")
  # MLE of lambda = 9.765054 

Left-truncation. For left-truncation the log likelihood is the following:

$$\log(\lambda) \sum_{i=1}^n y_i -n \lambda -n \log \left(1-{\sum_{j=0}^{k-1} {{e^{-\lambda}\lambda^j}\over{j!}}}\right) $$

with the sum on the far right being zero when $k<0$. Below is some R code:

# MLE estimation of a left-truncated Poisson with unknown truncation level

# Objective is to find estimates of lambda (underlying Poisson mean) and
# k (unknown trunction value).

# Generate some samples from a known truncated Poisson distribution
  lambda <- 10   # Poisson mean
  k <- 7         # Truncation level
  n <- 100       # Sample size
  y <- rpois(n*30,lambda) # Oversample to be more certain of getting n samples
  y <- y[y >= k][1:n]

# MLE for k and initial estimate for lambda
  khat <- min(y)
  lambda0 <- mean(y)  # Probably any starting value greater than zero will work

# Define log of the likelihood function
  logL <- function(lambda,ysum,n,k) {
       log(lambda)*ysum - n*lambda - n*log(1-ppois(k-1,lambda))
  }

# Find maximum likelihood estimate of lambda  
  trPoisson <- optim(lambda0, logL, ysum=sum(y), n=length(y), k=khat, 
    method="L-BFGS-B", lower=0, upper=Inf, control=list(fnscale=-1))

# Show results
  cat("MLE of truncation value =",khat,"\n")
  cat("MLE of lambda =",trPoisson$par,"\n")
$\endgroup$
9
  • $\begingroup$ khat <- max(y) lambda0 <- sum(y)/length(y[y < k]) the k here should be a khat no?Also: the random sample you create is this theoretically allowed? You basically just sample from a poisson and kick oyt the values not allowed. How can I proof that this is correct? $\endgroup$
    – Kathi
    Oct 17, 2015 at 11:22
  • $\begingroup$ If I define $k$ to be the lower bound truncation level (i.e. only values above $k$ are allowed, then using khat as defined above is not very useful. Would I use khat = min{y1,y2,..yn} then? $\endgroup$
    – Kathi
    Oct 17, 2015 at 11:39
  • $\begingroup$ Do you want truncation from the left or the right? Do you want truncation or censoring? $\endgroup$
    – JimB
    Oct 17, 2015 at 15:13
  • $\begingroup$ I want that if k = 5 for instance only values larger than 5 i.e 6,7,8,9... (i.e. only values to the right)are acceptable. So a lower truncation. I am not aware what censoring is... $\endgroup$
    – Kathi
    Oct 17, 2015 at 15:22
  • $\begingroup$ You are correct. The R code using k to assign a value to lambda0 should be khat. I'll fix that. $\endgroup$
    – JimB
    Oct 17, 2015 at 15:26

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.