Q: Fitting Poisson Distribution To Censored Data Repeat Experiments I am working on analyzing some data that come from the following experimental setup:
Imagine doing repeated trials of an experiment where a lever is pulled and a random number shows up drawn from a poisson distribution with parameter $\lambda$ (here representing the expected number/trial). Unfortunately you are really bad at reading numbers and so you are only able to ascertain if the number was 0 or $\ge 1$. What you are interested in though is deriving $\lambda$ from this censored data set.
For a single fixed value of $\lambda$ the current literature approach to solving this problem is to do the following: run the trials a large number of time and count how many times $\ge 1$ comes up and ratio that to the number of trials let's call this $\mu$. $\mu$ is then set equal to 1 - $p_0$ (the probability of getting 0) because it represents the sum of the probabilities of having gotten 1, 2, 3, etc and solved for $\lambda$ by the definition of the poisson distribution: $\lambda = -ln(1-\mu)$.
My question is this: is there a way to solve this problem if $\lambda$ is free to vary from trial to trial and I only care about getting the average value of $\lambda$?
As a note I realize the literature procedure doesn't work if $\lambda$ is large and so the experiments are usually configured to avoid that. 

A: Perhaps to belabor the obvious: The success of estimating $E(\lambda)$ depends on the variability
of $\lambda.$ 
I tried a simulation in R choosing $\lambda \sim \mathsf{Gamma}(6, rate = 2),$ so that $E(\lambda) = 3.$ Then generating Poisson values with
such values of $\lambda$ gives a sample of $Y_i$s with mean near 3, but
the distribution of $Y$ is not $\mathsf{Pois}(3),$ and in particular
$P(Y=0) \ne e^{-3}.$
set.seed(1234)
m = 10^6
y = rpois(m, rgamma(m, 6, 2))
mean(y);  var(y)
[1] 3.001025
[1] 4.505838
mean(y==0)
[1] 0.088019
dpois(0, 3)
[1] 0.04978707

hist(y, prob=T, br=(-1:21)+.5, ylim=c(0,.25), col="skyblue2")
  points(0:21, dpois(0:21, 3), col="red")


By contrast, if $\lambda \sim \mathsf{Unif}(2.9,3.0),$ so that
$E(\lambda) = 3,$ with very little variability, then it seems
possible to estimate $E(\lambda)$ from $P(Y = 0).$
set.seed(1235)
m = 10^6
y = rpois(m, runif(m, 2.9, 3.1))
mean(y);  var(y)
[1] 2.998847
[1] 3.004599
mean(y==0)
[1] 0.049865
dpois(0, 3)
[1] 0.04978707


