Calculate Lambda given

Question

I am using R and trying to calculate an unknown lambda. Known is an area in Poisson distribution.

x <- cbind("q" = rpois(n = 4, lambda = 3),  "ppois" = runif(n = 4), "lambda_unknown" = rep(NA, times = 4))
x
#      q      ppois lambda_unknown
# [1,] 4 0.05207818             NA
# [2,] 5 0.61127960             NA
# [3,] 3 0.83317758             NA
# [4,] 4 0.94495935             NA

I would like a function that helps me calculate my unknown lambda so "ppois" =ppois(q = q, "lambda_unknown")

desired (roughly)output:

x
#      q  ppois  lambda_unknown
# [1,] 4   0.05           ~9.15
# [2,] 5   0.61           ~5.02
# [3,] 3   0.83            ~1.4
# [4,] 4   0.94           ~2.05

Final row as an example:

ppois(4, 2.05)
# [1] 0.9427231

Matthew · Accepted Answer · 2016-01-20T22:55:56.603

Unfortunately, I don't think there is a built in function for that, but we can easily program a function to use a numeric approach.

For example, using the bisection method, we can do

approx<-function(q,p,epsilon) {
  lower<-0
  upper<-1
  while(ppois(q,upper)>p) {upper<-upper*2}
  while (upper-lower>epsilon) {
    middle<-(upper+lower)/2
    if (ppois(q,middle)<p) {upper<-middle}
    else {lower<-middle}
  }
  return ((upper+lower)/2)
}

This function will find the approximate value of lambda which results in a probability of p with a desired q within some desired epsilon (actually within epsilon/2). In order to use this we do have to use the fact that the ppois function is monotonically decreasing in lambda on the interval [0,infinity). It would still work, with modification, if the function were monotonically increasing, but monotonicity is required near our solution for the bisection method to work.

Using this

approx(4,0.05,0.01) # 9.152344
approx(5,0.61,0.01) # 5.035156
approx(3,0.83,0.01) # 2.144531
approx(4,0.94,0.01) # 2.082031

See here for more information on the bisection method. Other numeric methods are faster, but much harder to code.

In order to replace the lambda_unknown column with the desired values, we can use the apply function like so:

x[,"lambda_unknown"]<-apply(x,1,function(z){approx(z["q"],z["ppois"],0.01)})

This will apply the inline function to each row of the matrix x (the 1 indicates apply by rows, a 2 means apply by columns). The inline function takes the given row and computes the approx function feeding in the correct parameters from the row.

Works perfect, thanks a lot! One question, i probably miss something but I would like to use function like this. x$lambda <- approx(q = x$q, p = x$ppois, epsilon = 0.01). This command gives me only lambda for first row and 14 error messages that follows. Do I have to create a loop for this to work? — Tomas Ericsson, Jan 20 '16 at 20:48
You will have to use the apply function for that to work. For example, as your x is a matrix, we can do `x[,"lambda_unknown"]<-apply(x,1,function(x){approx(x[1],x[2],0.01)})`. — Matthew, Jan 20 '16 at 22:50

Calculate Lambda given

1 Answers1