Mean of a truncated Poisson random variable

Question

Let $Y$ be a Poisson random variable such that $E[Y]= \lambda$ for any $C\geq 1$ let truncation of $Y$ at $C$ be denoted by $X$ that is $X= min(Y,C)$.

How can I estimate the mean of $Y$ when I have only $X$?

I think I should give some background story to make my question more clear:

Let $Y$ be real demand for a specific product and the seller can only sell what she has in her warehouse which is $C$. Now the seller wants to know the mean of the true demand for her product but she has only sales data $X$ which has been historically censored by stock $C$. So $X= min(Y,C)$ might be more clear now.

Edit:

As @BruceET answer suggests an edit, I will answer his questions and I will post also my approach in this edit. (Thank you for your answer @AlgebraicsAnonymous).

I am aware of the difference between Truncation and Censoring but the paper that I am trying to apply on my data does not make this distinction (see section 5.1.1 in Amjad and Shah (2017) originaly published on SIGMETRICS)

The paper suggests two steps to estimate Y. At the firs steps they use "Universal Singular Value Thresholding" algorithm to estimate X because of missing values in the data Ref: Chatterjee (2015). Here is my code (I use soft thresholding):

usvd <- function(X,C,eta){

  # step 1 ------------------------------------------------------------------

  Y <- X

  Y[is.na(Y)]<- 0


  # step 2 ------------------------------------------------------------------

  A <- svd(Y)



  # step 3 ------------------------------------------------------------------


  P <- (1-sum(is.na(X))/length(X))

  # step 4 ------------------------------------------------------------------


  threshold <- function(N, eta) {

    my_nulls <- vector("numeric", length = length(N))

    S <- sign(N) * pmax(abs(N) - eta, my_nulls)

    S

  } 


  # step 5 ------------------------------------------------------------------


  W <- 1/P * A$u%*%diag(threshold(A$d,eta))%*%t(A$v)


   # step 6 ------------------------------------------------------------------



  M <- matrix(data = NA,ncol = ncol(Y),nrow = nrow(Y))

  for (i in 1:length(W)) {

    M[[i]] <- if(W[[i]]< 0){

      M[[i]] <- 0

    } else{if(W[[i]]>C[[i]]){

      M[[i]]<- C[[i]]

    }else{
      M[[i]]<-W[[i]]
    }

    }
  }

  return(M)

}

At the second step they use $M= [m_{ij}]_{i<=N,j<=T} $ where $m_{ij}=E[X_{ij}]$ and $C$ to get an estimate of $Y$ by using bisection method. That is where I need the mean of $Y$ see section 5.1.1 in Amjad and Shah (2017)

Here is my code:

lambda=colMeans(M)

real_Y <- function(M,C,lambda){

M <- M  
lambda <- lambda 

C <- C  
func <- function(C){

    #lambda <- colMeans(M)
   lambda=lambda
    M=M

    qk <- function(C) {C - ppois(C, lambda=lambda)}



    lambda * qk(C-1) / qk(C)-M

  }


  bisection <- function(f, a, b, n = 1000, tol = 1e-7) {
    # If the signs of the function at the evaluated points, a and b, stop the function and return message.
    if (!(f(a) < 0) && (f(b) > 0)) {
      stop('signs of f(a) and f(b) differ')
    } else if ((f(a) > 0) && (f(b) < 0)) {
      stop('signs of f(a) and f(b) differ')
    }

    for (i in 1:n) {
      c <- (a + b) / 2 # Calculate midpoint

      # If the function equals 0 at the midpoint or the midpoint is below the desired tolerance, stop the 
      # function and return the root.
      if ((f(c) == 0) || ((b - a) / 2) < tol) {
        return(c)
      }

      # If another iteration is required, 
      # check the signs of the function at the points c and a and reassign
      # a or b accordingly as the midpoint to be used in the next iteration.
      ifelse(sign(f(c)) == sign(f(a)), 
             a <- c,
             b <- c)
    }
    # If the max number of iterations is reached and no root has been found, 
    # return message and end function.
    print('Too many iterations')
  }

# print the result --------------------------------------------------------


 print(bisection(func,M,C))

}

Ref: for bisection method

About lost sales: I think the lost sales caused by stock level. And the data has been polluted all the time because of the low-stock so we do not want to change level of $C$ before having an estimate of true demand $Y$.

Answer 1

Note that if $Y$ is a one-dimensional random variable (ie. treating the one-dimensional case first), then for $k \in \mathbb N$ we have $$ P(\min\{Y, C\} = k) = \begin{cases} 0 & k > C \\ \sum_{j \ge C} e^{-\lambda} \frac{\lambda^j}{j!} & k = C \\ e^{-\lambda} \frac{\lambda^k}{k!} & k < C \end{cases}; $$ for simplification, we assume $C \in \mathbb Z$ here. Thus, in order to get the distribution (and hence the expactation) of $Y$, it is necessary and sufficient that each entry of the matrix $C$ is bigger than two. For, you may then read of the coefficient $\lambda_{i,j}$ (ie. the coefficient of the Poisson distribution in the $(i,j)$-th entry) from the equation $$ \alpha_{i,j} = e^{\lambda_{i,j}} \lambda_{i,j}, $$ where $\alpha_{i,j} = P(\min\{Y_{i,j}, C_{i,j}\} = 1)$; for this you do need the Lambert W-function.

Answer 2

Comment: It is important to distinguish between truncation and censoring. Roughly speaking, the distinction is that in truncation, one has no knowledge of subjects with values beyond $C$ or even how many there were, whereas there are various kinds of censoring (see link) in which one knows the true sample size, but for some reason does not see values beyond $C.$

You have not given quite enough information for us to know whether the sales data correspond to a truncated sample or a censored sample. You use both of the words truncation and consoring in your question.

---

Suppose that $Y \sim \mathsf{Pois}(\lambda = 10)$ and that $C = 15.$ Then $P(Y \le 15) = .9513,$ and the seller is losing very little.

ppois(15, 10)
[1] 0.9512596

While $E(Y) = 10,$ you have $E(X) \approx 9.64,$ according to the simulation below.

y = rpois(10^6, 10);  x = y[y <= 15];  mean(x)
[1] 9.637339

With the probabilities $p_i = P(Y=i) = e^{-10}10^i/i!,$ for $i = 0,1,2, \dots ,$ the mean $E(Y) = \sum_{i=0}^\infty ip_i = 10$ can be found very accurately by summing only the first 101 terms:

i = 0:100;  sum(i*dpois(i, 10))
[1] 10

For the random variable $X,$ we use slightly inflated probabilities $p_j^* = p_j/0.9512596$ and sum only the terms for $j = 0, 1, \dots, C=15.$ Then $E(X) = \sum_{j=1}^{15} jp_j^* \approx 9.635,$ which is the same result we got from simulation above.

j = 0:15;  sum(j*dpois(j, 10)/ppois(15, 10))
[1] 9.635031

The average of the daily sales in the records can be used to estimate $E(X).$

---

If this is a textbook exercise, please edit your Question, telling us the topic of the relevant chapter, stating your problem more explicitly, and showing what approaches you have tried.

If this is a real-life inventory problem then there may be two approaches: (a) Try to understand what triggers lost sales. Does the client place orders too large to fill? Are orders filled partially in case of shortage. Or, maybe the sellers just put up a sign saying, "Out of Widgets. Closed for Today. Sorry for Any Inconvenience." In short, better communication with customers may provide a better answer than you can get by statistical methods. (b) Although the mean of an (unrestricted) sample is the best way to estimate $E(Y) = \lambda,$ depending on $\lambda$ and $C$ there other methods to estimate $\lambda$ from $X$-values. With better understanding of the exact situation, perhaps one of us can suggest a practical procedure for estimation.