Estimating the parameters of a Poisson times a constant 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$).

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
