R: How can I create expected values of a variable assuming Poisson distribution? Working in R, I have a dataset with values which I would like to compare to expected values if it was a Poisson distribution. Is there a way to do it?
Example:
n <- c(1,2,3,4,5,6,7,8,9,10)
people <- c(850, 200, 100, 60, 40, 20, 25, 10, 7, 60 )

df <- data.frame(n, people)

In my case, n is number of times event happened and people is number of people (10 means attending 10 or more events). I would like to see how many people would be in each category assuming poisson distribution.
I am completely clueless how to approach this.
 A: From your data and the context of your question, it appears that this is a problem where you are actually observing data from a censored and truncated Poisson distribution, where you don't observe people who go to zero events and you observe people with ten or more events in one category.  Assuming this distributional form with a maximum observation of $\dot{x}=10$ gives you the following probability mass function for an individual value:
$$p_X(x) = 
\begin{cases}
\frac{e^{-\lambda}}{1-e^{-\lambda}} \cdot \frac{\lambda^x}{x!} 
& & & \text{for }x=1,2,3,..., \dot{x}-1, \\[6pt]
1 - \frac{e^{-\lambda}}{1-e^{-\lambda}} \cdot \sum_{i=1}^{\dot{x}-1} \frac{\lambda^i}{i!}  
& & & \text{for }x=\dot{x}. \\[6pt]
\end{cases}$$
Suppose you observe $n$ IID data points $x_1,...,x_n$ from this distribution.  To facilitate analysis, let $\dot{n} = \sum_{i=1}^n \mathbb{I}(x_i=\dot{x})$ be the number of censored points and $\bar{x}_n = \sum_{i=1}^n x_i \mathbb{I}(x_i<\dot{x})/(n-\dot{n})$ be the sample mean of the non-censored points.  The log-likelihood function for this data is:
$$\begin{align}
\ell_\mathbf{x}(\lambda)
= \text{const} 
&+ (n-\dot{n}) \bigg[ \bar{x}_n \log (\lambda) - \lambda - \log (1-e^{-\lambda}) \bigg] \\[6pt]
&+ \dot{n} \log \bigg( 1 - \frac{e^{-\lambda}}{1-e^{-\lambda}} \cdot \sum_{i=1}^{\dot{x}} \frac{\lambda^i}{i!} \bigg).
\end{align}$$
This function can be maximised numerically to get the maximum likelihood estimator (MLE).  The statistic $(n, \dot{n}, \bar{x}_n)$ is a sufficient statistic in this distribution, so we can create a function to find the MLE that takes either the full dataset or this summary of the data.  It is useful to create a function to compute the MLE of the rate parameter for IID data from a censored Poisson distribution.  Here we give a relatively simple function for this task, with the optimisation performed on the parameter $p = \log(\lambda)$ for purposes of numerical stability.
dpois.ct <- function(x, xmax, lambda, log = FALSE) {
  
  #Check input
  if (!is.numeric(x))   stop('Input x should be a numeric vector')
  
  #Compute log-probabilities
  LOGPROBS <- rep(-Inf, length(x))
  for (i in 1:length(x)) {
    if (x[i] %in% 1:xmax) {
      LOGPROBS[i] <- dpois(x[i], lambda, log = TRUE) }
    if (x[i] == xmax) {
      LOGPROBS[i] <- ppois(xmax-1, lambda, lower.tail = FALSE, log = TRUE) } }
  LOGPROBS <- LOGPROBS - VGAM::log1mexp(lambda)
  
  #Return output
  if (log) { LOGPROBS } else { exp(LOGPROBS) } }


MLE.pois.ct <- function(x, xmax, ...) {
  
  #Set objective function and compute MLE
  NEGLOGLIKE <- function(p) { 
    LL <- dpois.ct(x, xmax, lambda = exp(p), log = TRUE)
    -sum(LL) }
  MLE <- exp(nlm(NEGLOGLIKE , p = log(mean(x)), ...)$estimate)
  names(MLE) <- 'MLE.rate'
  
  #Give output
  MLE }

We can implement this for your data to get the MLE and produce a corresponding barplot of the estimated distribution.  We first generate your data and use the MLE.pois.ct function to compute the MLE.  From the output below we see that a reasonable estimate of the rate parameter in your problem is $\hat{\lambda} = 1.876321$.  The barplot shows the estimated probabilities under the model (the blue bars) against the actual relative frequencies in your data (the black dots).  As you can see from the barplot, your data do not appear to follow a censored and truncated version of the Poisson distribution, so your model assumption seems unreasonable here.
#Generate the data vector
x <- rep(0, sum(people))
i <- 1
p <- 0
while (i <= length(people)) {
  x[(p+1):(p+people[i])] <- n[i]
  p <- p+people[i]
  i <- i+1 }

#Compute the sample mean of your data and the MLE
MLE.rate <- MLE.pois.ct(x, xmax = 10)
MLE.rate

MLE.rate 
1.876321

#Compute estimated probabilities in censored Poisson distribution
PROBS <- dpois.ct(n, xmax = 10, lambda = MLE.rate)
names(PROBS) <- n

#Barplot of estimated distribution
BARPLOT <- barplot(PROBS, col = 'blue', ylim = c(0,1),
           main = 'Estimated censored-truncated Poisson distribution',
           xlab = 'Number of Events', ylab = 'Estimated Probability')
points(x = BARPLOT , y = people/sum(people), pch = 16, cex = 1.2)


