2
$\begingroup$

Suppose that we have a distribution $X = aY$ where the distribution of $Y$ is Poisson{$\lambda$}, what is the maximum likelihood estimator of $a$?

By the method of moments we can estimate $a$ as $\sigma^2_X /\mu_X $. This is because $\mu_X = a \lambda$ and $\sigma^2_X = a^2 \lambda$.

The log-likelihood function is I think:

$$ l(\lambda,a|x_j) = \log \prod_{j=1}^{n} e^{-\lambda} \frac{\lambda^{x_j/a}} {(x_j/a)!} $$

$$ =-n\lambda + {\frac{1}{a}}\log{\lambda}\sum_{j=1}^n{x_j}-\sum_{j=1}^n \log{\left(\frac{x_j}{a}!\right)} $$

However I am having trouble solving this for $a$ after differentiating and setting to zero.

The motivation for this question is that I am trying to learn how to do maximum likelihood for this distribution. I know that $a$ can be estimated by finidng the greatest common divisor of the differences of the observed values of $X$ but I am not interested in finding $a$ that way. Especially because I have some real-world motivation for this question and actual observations would be distributed as $X + \epsilon$ where $\epsilon \sim(0,\sigma_\epsilon^2$).

$\endgroup$
5
  • 3
    $\begingroup$ Unless you have only a handful of observations you can, with extremely high confidence, estimate $a$ with perfect accuracy as the greatest common divisor of the differences among all the values. Upon dividing all values by $a,$ you will be in a standard Poisson estimation setting. If the noise affects the values as well as the counts, then we need more information about the nature of that noise: it might need to be incorporated into the probability model. $\endgroup$
    – whuber
    Commented Dec 1, 2022 at 0:25
  • $\begingroup$ Hrm, suppose there is sufficient noise such that say $a=0.1$, $\lambda=5$ can't be easily distinguished from $a=0.2$, $\lambda=2.5$ by this method. I am thinking normally distributed noise on the values so that the data looks continuous. Basically I am looking for some way to do the parameter estimation based on the shape of the distribution, not on any properties stemming from its discrete-ness. $\endgroup$
    – Ben S.
    Commented Dec 1, 2022 at 1:41
  • 2
    $\begingroup$ Those two scaled Poisson cases can easily be distinguished with high reliability in all but extremely small samples. Can you clarify what you mean by 'the data are too noisy for that' -- either they're integers scaled by a constant, or they aren't. If they are, whubers comment applies (where can this "noise" come in?). If they aren't scaled integers, don't use that model at all. The shape of the distribution is discrete (the Poisson literally has spikes at every integer), so using its discreteness is simply making use of the facts of its actual shape. $\endgroup$
    – Glen_b
    Commented Dec 1, 2022 at 2:00
  • $\begingroup$ I guess I was just hoping there is a way to leverage the intuition that every combination of $a$ and $\lambda$ creates a unique shape to the distribution. $\endgroup$
    – Ben S.
    Commented Dec 1, 2022 at 2:28
  • 1
    $\begingroup$ The method and the results will vary with the ratio $a/\sigma_\epsilon.$ When that ratio exceeds $1$ or so, $a$ will usually be clearly and uniquely estimable with no error. Otherwise, how to proceed might depend on $\lambda:$ large $\lambda$ can use accurate approximations. Thus, it would help to know what you expect $a,$ $\sigma_\epsilon,$ and $\lambda$ to be (within reasonable ranges) as well as how much data you typically have. $\endgroup$
    – whuber
    Commented Dec 1, 2022 at 20:34

2 Answers 2

4
$\begingroup$

Gaussian iid errors: If the errors are independent and normally distributed, the likelihood can be expressed as \begin{align} L(\lambda,\sigma,a) &=\prod_{i=1}^n f_{Y_i}(y_i), \\&=\prod_{i=1}^n\sum_{x=0}^\infty f_{Y_i|X_i}(y_i|x)P(X_i=x), \\&=\prod_{i=1}^n\sum_{x=0}^\infty \frac1{\sqrt{2\pi}\sigma}\exp\left(-\frac{(y_i-ax)^2}{2\sigma^2}\right)\frac{e^{-\lambda} \lambda^x}{x!} \end{align} via the law of total probability. In practice only a finite number of terms of the inner infinite sum needs to be computed.

This likelihood may have multiple optima but for a given $a$, $L(\lambda,\sigma,a)$ can only have a single optimum for $\lambda$ and $\sigma$. This suggest working with the profile likelihood for $a$, $$ L_p(a)=\sup_{\lambda,\sigma}L(\lambda,\sigma,a) $$ or its logarithm computed and shown below (lower panels) for simulated data (upper panels).

Unless the standard deviation $\sigma$ of the errors is too large, the range of likely values for $a$ as expected is very narrow centered around the true value (left panels, green vertical line), but for small $n$ and large $\sigma$, several values of $a$ may be plausible and the model may become nearly unidentifiable (if $\sigma$ is an unknown parameter in addition to $a$ and $\lambda$) as judged by $L_p(a)$ being nearly flat for small $a$ (right panels).

Rounding error: A similar approach can be taken if the errors are due to rounding off to a certain (known) number of decimal places. This would lead to a likelihood being piecewise constant in intervals of possible values of $a$ and zero elsewhere. Such rounding errors would not be independent so it seems questionable if they can be dealt with using the above likelihood based on iid normal errors.

profileLogLikelihood<- function(a, y) {
  logLik <- function(par) {
    lambda <- exp(par[1])
    sigma <- exp(par[2])
    
    res <- 0
    # Compute range of x values contributing something within the machine precision
    # based on a Gaussian approximation in x of the terms in the sum
    prec <- c(a/sigma^2, 1/lambda)
    sd <- sqrt(1/sum(prec))
    nsd <- 1.5*sqrt(-2*log(.Machine$double.eps)) # 
    
    for (i in 1:length(y)) {
      xhat <- sum(prec*c(y[i]/a, lambda))/sum(prec)
      x <- max(0, floor(xhat - nsd*sd)):ceiling(xhat + nsd*sd)

      # log of each term in inner sum
      logs <- dnorm(y[i], a*x, sigma, log=TRUE) + dpois(x, lambda, log=TRUE)
      # log sum exp trick 
      maxlogs <- max(logs)
      logsumexp <- maxlogs + log(sum(exp(logs - maxlogs)))
      # add contribution to the total log likelihood
      res <- res + logsumexp
    }
    res
  }
  # profile out lambda and sigma
  lambdahat <- mean(y)/a
  optim(c(log(lambdahat),0), logLik, control=list(fnscale=-1,reltol=1e-12,maxit=1000))$value
}
profileLogLikelihood <- Vectorize(profileLogLikelihood, vectorize.args = "a")

example <- function(a, n, sd, seed) {
  set.seed(seed)
  x <- rpois(n, 2)
  y <- a*x + rnorm(n,sd=sd)
  plot(density(y,bw = sd),main="")
  points(y,rep(0,n),pch="+")
  abline(v=a*(0:50),col="grey")
  curve(profileLogLikelihood(a, y=y), xname="a", 0.1, 1, n=500, ylab="profile logLik(a)")
  abline(v=a, col="green")
}
par(mfcol=c(2,2))
example(a=.4, n=50, sd=.05, seed=1)
example(a=.4, n=10, sd=.1, seed=2)

Created on 2022-12-01 with reprex v2.0.2

$\endgroup$
3
  • $\begingroup$ Wow thank you, this is super interesting! The plots are quite fascinating. I look forward to digging into your code, I have never used optim(). $\endgroup$
    – Ben S.
    Commented Dec 2, 2022 at 14:53
  • $\begingroup$ Can you tell, if there was no noise, would there be an ML estimator of $a$ that is different from var/mean? $\endgroup$
    – Ben S.
    Commented Dec 2, 2022 at 14:56
  • 1
    $\begingroup$ Yes, as noted by @whuber, the greatest common divisor will usually be the MLE of $a$. An exception would be that the $X_i$ happen to take only even integer values in which case the MLE of $a$ may be half of the greatest common divisor of the $Y_i$'s but outcomes like that are as you can imagine very unlikely to happen except for small sample sizes $n$. $\endgroup$ Commented Dec 2, 2022 at 15:18
0
$\begingroup$

EDIT: this was an answer to the original form of the question.

A simple answer was staring me in the face although I am not sure if this is likely not the maximum likelihood estimator of $a$ (will happily select an answer that can). If $Y\sim\textrm{Poisson}(\lambda) $ then $\mu_Y=\lambda$ and $\sigma^2_Y = \lambda$. If $X=aY$ then $\mu_X=a\lambda$ and $\sigma^2_X = a^2\lambda$. Hence $\sigma^2_X/\mu_X$ estimates $a$.

Example in R:

> Y = rpois(1000,0.7)
> X = 0.2*Y
> var(X)/mean(X)
[1] 0.2109116
$\endgroup$
4
  • 1
    $\begingroup$ This is the method of moments estimator of 𝑎 but it is certainly not the maximum likelihood estimator and probably not very efficient for small sample sizes. I guess what you call noise in the data involves at least some error due to rounding off the a certain number of decimal places in which case it is straightforward to write down the likelihood. Or is there anything more to it than that? $\endgroup$ Commented Dec 1, 2022 at 14:20
  • $\begingroup$ I would think of the noise as adding $\epsilon \sim N(0,\sigma^2$) to each value but rounding might be an OK model too. I have been trying to write down the likelihood but have trouble solving it for $a$. $\endgroup$
    – Ben S.
    Commented Dec 1, 2022 at 14:34
  • 1
    $\begingroup$ @whuber I think the question is clear enough given information given in the comments, so I'm voting to reopen. I was about to post an answer. $\endgroup$ Commented Dec 1, 2022 at 14:39
  • $\begingroup$ Eagerly awaiting! I don't think the question should be closed, even if one doesn't see the point to the question it's a valid question that can be answered (as my own var/mean answer shows) $\endgroup$
    – Ben S.
    Commented Dec 1, 2022 at 15:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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