Reverse Poisson?

[largeselection]

Not sure if anyone can help with this. Basically I'm trying to find out how to apply the reverse of the Poisson function in excel. So as of now I have poisson(x value, mean, true-cumulative) and that lets me get the probability for that occurence. Basically I want to know how I can get the minimum/maximum x value based on a given probability.

So if I have a list of data (700 rows) and I want to find out what the minimum starting value should be given a desired average and the fact that I want the lowest value to be at the 0.05% probability. So 0.05% = (x, 35, True) solve for x.

I know I can prob do this with solver, but I am trying to figure out a way to do this formulaicly without having to use the solver (as I may have to use this many times)

Thanks for any insight/links/help as always!

 

[shg]

Do a lookup against the CDF?

       --A--- --B--- -C-- -----------------------D------------------------
   1   lambda        Prob                                                 
   2       9           5%                                                 
   3     x     CDF                                                        
   4       0    0.0%      B4 and down: =POISSON(x, lambda, TRUE)          
   5       1    0.1%    3 C5: =INDEX(x, MATCH(Prob, PoissonCum))          
   6       2    0.6%   13 C6: {=INDEX(x, MATCH(Prob, 1 - PoissonCum, -1))}
   7       3    2.1%                                                      
   8       4    5.5%                                                      
   9       5   11.6%                                                      
  10       6   20.7%                                                      
  11       7   32.4%                                                      
  12       8   45.6%                                                      
  13       9   58.7%                                                      
  14      10   70.6%                                                      
  15      11   80.3%                                                      
  16      12   87.6%                                                      
  17      13   92.6%                                                      
  18      14   95.9%                                                      
  19      15   97.8%                                                      
  20      16   98.9%                                                      
  21      17   99.5%                                                      
  22      18   99.8%                                                      
  23      19   99.9%                                                      
  24      20  100.0%

 

[shg]

New solution:

Function PoissonInv(dVal As Double, dMean As Double) As Variant
    ' shg 2011
 
    ' For a Poisson process with mean dMean, _
      ' returns the largest integer such that the CDF <= dVal
 
    ' E.g., =POISSON(5, 10, TRUE) returns 0.0670859629
    ' PoissonInv(0.0670859629, 10) returns 5
 
    Dim iX          As Long
    Dim dCDF        As Double
 
    ' these variables are used to simplify this summation:
    '    dCDF = dCDF + Exp(-dMean) * dMean ^ iX / .Fact(iX)
    Dim dExpMean    As Double   ' =Exp(-dMean)
    Dim dFact       As Double   ' incremental factorial
    Dim dPowr       As Double   ' incremental power
 
    If dVal < 0 Or dVal >= 1 Then
        PoissonInv = CVErr(xlErrValue)
 
    ElseIf dVal > 0 Then
        dExpMean = Exp(-dMean)
        dFact = 1
        dPowr = 1
 
        Do While dCDF < dVal
            dCDF = dCDF + dExpMean * dPowr / dFact
            iX = iX + 1
            dFact = dFact * iX
            dPowr = dPowr * dMean
        Loop
 
        PoissonInv = iX - IIf(dCDF / dVal > 1.000000000001, 2, 1)
    End If
End Function

 

[jerestat]

There is a problem with this though. With the Poisson Distribution, x can never be negative as is it impossible to have a negative number of occurrences of something. Also, I'm not sure that I understand why there is a need to correct at the bottom.

Make one minor change:
PoissonInv = Application.WorksheetFunction.Max(iX - IIf(dCDF/dval > 1.00000000000001, 2, 1), 0) 'there can be zero occurrences but no less... +

As far as I can tell it is working now. I have tries looping through a number of different scenarios and it seems to be working well except for when lambda is very low and X is relatively high. Of course, if the lamda is low and the X is high then the pdf will be very close to zero and the CDF very close to 1. So this, would not be too much of a problem I don't think.

 

[shg]

Can you give an example that requires that change?

 

[jerestat]

Actually I think the problem was in the looping construct. It seems to be coming up with an error at the same test. I'm not sure whether or not it is the random number combination or the loop itself. (Does the seed reset itself, that may be why?)

Sub TestPoissonInv()    Dim topMean, topLambda As Integer, x, lambda As Double, i, j, k As Long, ws As Worksheet
    'topMean is the the highest value that the mean can be
    j = 2
    Set ws = ActiveSheet
    ws.Cells(1, 1) = "X"
    ws.Cells(1, 2) = "lambda"
    ws.Cells(1, 3) = "Poisson.Dist"
    ws.Cells(1, 4) = "PoissonInv"
    ws.Cells(1, 5) = "Check"
        For topLambda = 1 To 50
            For topMean = 1 To 100
                For i = 2 To 10 ^ 3
                    x = Application.WorksheetFunction.RoundDown(Rnd() * topMean, 0)
                    lambda = Rnd() * topLambda
                    ws.Cells(j, 1) = x
                    ws.Cells(j, 2) = lambda
                    ws.Cells(j, 3) = "=POISSON.DIST(" & x & "," & lambda & ", TRUE)"
                    ws.Cells(j, 4) = IIf(ws.Cells(j, 3) = 1, "N/A", "=POISSONINV(" & ws.Cells(j, 3) & "," & lambda & ")")
                    ws.Cells(j, 5) = (ws.Cells(j, 1) = ws.Cells(j, 4))
                    j = j + 1
                Next i
            Next topMean
        Next topLambda
End Sub

 

[shg]

That routine just got a cleanup and improved description:

Function PoissonInv(Prob As Double, Mean As Double) As Variant
  ' shg 2011, 2012-0519, 2014-0727, 2015-0414

  ' For a Poisson process with mean Mean, returns a three-element array:
  '   o The smallest integer N such that POISSON(N, Mean, True) >= Prob
  '   o The CDF for N-1 (which is < Prob)
  '   o The CDF for N (which is >= Prob)

  ' E.g., POISSON(5, 10, TRUE) returns 0.067085962
  ' PoissonInv(0.067085962, 10) returns 5

  ' Returns a descriptive error if
  '   Prob <= 0 or Prob >= 1
  '   Exp(-Mean) = 0 (i.e., Mean > ~745)
  '   The factor dK overflows

  Dim N             As Long     ' number of events
  Dim CDF           As Double   ' cumulative distribution function at N
  Dim CDFOld        As Double   ' CDF at N-1

  ' These two variables are used to simplify the
  ' probability mass function summation:
  '    CDF = CDF + Exp(-Mean) * Mean ^ N / N!
  Dim dExpMean      As Double   ' =Exp(-Mean)
  Dim dK            As Double   ' incremental power & factorial

  If Prob <= 0 Or Prob >= 1 Then
    PoissonInv = "Prob ]0,1["

  ElseIf Mean > 745 Then
    PoissonInv = "Mean ]0,745]"

  Else
    dExpMean = Exp(-Mean)

    If dExpMean > Prob Then
      ' Prob < POISSON(0, Mean, TRUE)
      PoissonInv = Array(0#, 0#, dExpMean)

    Else
      dK = 1#
      On Error GoTo Oops

      Do
        CDFOld = CDF
        CDF = CDF + dExpMean * dK
        N = N + 1
        dK = dK * Mean / N
      Loop While CDF < Prob - 0.000000000000001

      ' we return N-1 because N has already been incremented
      PoissonInv = Array(N - 1, CDFOld, CDF)
    End If
    Exit Function

Oops:
    PoissonInv = "dK overflow!"
  End If
End Function

 

[shg]

A further mod to accommodate large means:

Function PoissonInv(Prob As Double, Mean As Double) As Variant
  ' shg 2011, 2012, 2014, 2015-0415

  ' For a Poisson process with mean Mean, returns a three-element array:
  '   o The smallest integer N such that POISSON(N, Mean, True) >= Prob
  '   o The CDF for N-1 (which is < Prob)
  '   o The CDF for N (which is >= Prob)

  ' E.g., POISSON(5, 10, TRUE) returns 0.067085962
  ' PoissonInv(0.067085962, 10) returns 5

  ' Returns a descriptive error if Prob <= 0 or Prob >= 1
  ' If Mean >= 100, then uses a normal approximation

  Dim N             As Long     ' number of events
  Dim CDF           As Double   ' cumulative distribution function at N
  Dim CDFOld        As Double   ' CDF at N-1

  ' These two variables are used to simplify the probability mass
  ' function summation CDF = CDF + Exp(-Mean) * Mean ^ N / N!
  
  Dim dExpMean      As Double   ' =Exp(-Mean)
  Dim dK            As Double   ' incremental power & factorial

  If Prob <= 0 Or Prob >= 1 Then
    PoissonInv = "Prob ]0,1["

  ElseIf Mean < 100 Then
    dExpMean = Exp(-Mean)
    dK = 1#
    CDF = dExpMean

    Do While CDF < Prob - 0.000000000000001
      CDFOld = CDF
      N = N + 1
      dK = dK * Mean / N
      CDF = CDF + dExpMean * dK
    Loop

    PoissonInv = Array(N, CDFOld, CDF)

  Else
    ' Plan B, for large means; approximate the Poisson as a normal distribution
    ' http://en.wikipedia.org/wiki/Continuity_correction#Poisson
    Dim iInv        As Long

    With WorksheetFunction
      iInv = .Ceiling(.Norm_Inv(Prob, Mean, Sqr(Mean)) - 0.5, 1)
      PoissonInv = Array(iInv, _
                         .Norm_Dist(iInv - 0.5, Mean, Sqr(Mean), True), _
                         .Norm_Dist(iInv + 0.5, Mean, Sqr(Mean), True))
    End With
  End If
End Function

 

[willwillpower]

I found a formula-based solution to this that would work for most situations where the mean is not huge. The Poisson is really just the limit of a Binomial Distribution with very large number of trials.

So you get the same result as this User-Defined formula with the in-built formula "=BINOM.INV(100000,<mean>/100000,<alpha>)"

 

[willwillpower]

Sorry, that didn't come out right. The formula should read:
"=BINOM.INV(100000,mean<mean>/100000,alpha<alpha>)"</alpha></mean>

I found a formula-based solution to this that would work for most situations where the mean is not huge. The Poisson is really just the limit of a Binomial Distribution with very large number of trials.

So you get the same result as this User-Defined formula with the in-built formula "=BINOM.INV(100000,<mean>/100000,<alpha>)"

</alpha></mean>