Coverage probability for a Poisson parameter

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
13 messages Options
Reply | Threaded
Open this post in threaded view
|

Coverage probability for a Poisson parameter

Angle
This post has NOT been accepted by the mailing list yet.
Hi, I am stuck while trying to compute the coverage probability for a variable Poisson parameter.

I want to plot the coverage probability of the Wald Interval as a function of np, in order to show that it performs badly. Since the sum of poisson distributed random variables is a poisson with a new parameter corresponding to the sum of the parameters of the random variables, I generated data x from a poisson with parameter np.

I just get a graph where the line of the coverage probability is stable to 0, while it should fluctuate around 0.95.

My aim is to express the coverage probability as the frequency of times that the true parameter is inside the Wald Interval over m simulations. But I am not really good at creating loops and I don't know if I need one loop in order to achieve what I am aiming to.

m <- 1000000
conf.level <- 0.95
np <- 2:50
x <- rpois(m, np)

cover <- function(np) {
  stopifnot(is.numeric(np))
  stopifnot(length(np) == 1)
  stopifnot(2 <= np & np <= 50)
  fpx <- dpois(x, np)
  phat <- x / m
  low <- phat - qnorm((1 + conf.level) / 2) * sqrt(phat / m)
  hig <- phat + qnorm((1 + conf.level) / 2) * sqrt(phat / m)
  inies <- as.numeric(low <= np & np <= hig)
  sum(inies * fpx)
}

np <- seq(2, 50, 1)
plot(np, vapply(np, cover, 0.5), type = "l", ylab = "coverage probability")
abline(h = conf.level, lty = 2)
Reply | Threaded
Open this post in threaded view
|

Re: Coverage probability for a Poisson parameter

JS Huang
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Coverage probability for a Poisson parameter

Angle
This post has NOT been accepted by the mailing list yet.
This post was updated on .
Thank you, I am starting to understand now.

Your code sets the function cover that every time generates observations from a poisson distribution and then compute the confidence intervals and the coverage, right?

And with sapply does it repeat the simulation 10000 times and then get the mean of the coverage probabilities?

My final aim is to realize two graphs where I can plot coverage probabilities: the first graph for fixed n and the second for fixed lambda.

In order to do that, I think I have to initialise a for loop, where I make the value of lambda varying from a to b and the same in another for loop.

Say, I want to hold fixed lambda at 0.5 and n from 1 to 100.

significance.level <- 0.05
lambda <- 0.5
for (sample.size in 1:100){
     cover <- function(lambda, sample.size, significance.level)  {
     x <- rpois(sample.size,lambda)
     estimate <- mean(x)
     lower <- estimate - qnorm(1 - significance.level/2) * sqrt(estimate/sample.size)
     upper <- estimate + qnorm(1 - significance.level/2) * sqrt(estimate/sample.size)
     if (lambda > lower & lambda < upper){1}else{0}
     mean(sapply(1:100000, function(x)cover(lambda,sample.size,significance.level)))
     }
}


I would like to obtain the mean of the coverage probabilities for (n=1, ..., n=100) and then plot the coverage probabilities, but this loop doesn't seem to work properly. I am a beginner with R and I still need to learn how it works.

While for the graph, I figured out, but first I would like to compute correctly the coverage probabilities obviously.

sample.size <- seq(1, 100, 1)
plot(sample.size, vapply(sample.size, cover, 0.5), type = "l", ylab = "coverage probability")
abline(h = conf.level, lty = 2)
Reply | Threaded
Open this post in threaded view
|

Re: Coverage probability for a Poisson parameter

JS Huang
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Coverage probability for a Poisson parameter

Angle
This post has NOT been accepted by the mailing list yet.
This post was updated on .
EDIT: I solved the problem, it was just a typo.


Thank you, that solved many of my problems! With n I meant the sample size.

With your function I computed the coverage probability of many confidence intervals.

Now I am a bit stuck with the confidence interval derived from the likelihood ratio test.

I tried this way but the graph plotted does not behave as I would expect. The bound I use should be right, but it doesn't display the expected results and I don't get where I am wrong.

#lrt Poisson
#FIXED N AND LAMBDA FROM 0.1 TO 5
n<-10
significance.level <- 0.05
cover <- function(lambda, n, significance.level)  {
  x <- rpois(n,lambda)
  estimate <- mean(x)
  ratio <- -2*log(((lambda^(n*estimate))*exp(-n*lambda)))/((estimate^(n*estimate))*exp(-n*estimate))
  if (ratio <= (qnorm(significance.level/2))^2){1}else{0}
}

resulting.matrix <- matrix(0, nrow=500,ncol=2)
for (i in 1:500) {
  resulting.matrix[i,1] <- 0.01 * i
  resulting.matrix[i,2] <- mean(sapply(1:100,function(x) cover(0.01*i,10,0.05)))
}

resulting.matrix
lambda <- seq (0.01, 5, 0.01)
nlambda <- n*lambda
plot(nlambda, resulting.matrix[1:500,2], type ="l", ylab = "coverage probability")
abline(h = 0.95, lty = 2)
Reply | Threaded
Open this post in threaded view
|

Re: Coverage probability for a Poisson parameter

JS Huang
This post has NOT been accepted by the mailing list yet.
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Coverage probability for a Poisson parameter

JS Huang
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Coverage probability for a Poisson parameter

Angle
This post has NOT been accepted by the mailing list yet.
Thank you very much, it was really useful and I could easily edit the methology to adapt it to many other cases!

I moved on and now I am trying to find confidence interval for the ratio of two Poisson means. Given S1 and S2 distributed as poisson with parameters lambda1 and lambda 2, I know that S2|(S1+S2) is distributed as binomial with parameters(S1+S2,lambda2/(lambda1+lambda2))

I am trying to found, as above, the coverage probability of the Wilson interval, for the parameter lambda2/(lambda1+lambda2)) and then tried to slighty modify the confidence interval in order to have the confidence interval for lambda2/lambda1, my actual parameter of interest.

This is my code but it gives me this error and I don't understand what I am doing wrong.

"""""Error in if (theta >= lower & theta <= upper) { :
  missing value where TRUE/FALSE needed
5 cover(0.001 * i, 20, 0.05)
4 FUN(1:100[[53L]], ...)
3 lapply(X = X, FUN = FUN, ...)
2 sapply(1:100, function(x) cover(0.001 * i, 20, 0.05))
1 mean(sapply(1:100, function(x) cover(0.001 * i, 20, 0.05))) """""""

#Wilson Poisson Ratio
n <- 20
significance.level <- 0.05
cover <- function(lambda, n, significance.level)  {
  s1 <- rpois(1,lambda1)
  s2 <- rpois(1,lambda2)
  theta <- lambda2/(lambda1+lambda2)
  s <- s1+s2
  z <- qnorm(1-0.05/2)
  k <- z^2
  pi <- s2/s
  lower <- (pi+(k/(2*s))-z*sqrt((pi*(1-pi)+(k/4*s))/s))/(1+k/s)
  upper <- (pi+(k/(2*s))+z*sqrt((pi*(1-pi)+(k/4*s))/s))/(1+k/s)
  if (theta >= lower & theta <= upper){1} else {0}
}
resulting.matrix <- matrix(0, nrow=500,ncol=2)
for (i in 1:500) {
  resulting.matrix[i,1] <- 0.001 * i
  resulting.matrix[i,2] <- mean(sapply(1:100,function(x) cover(0.001*i,20,0.05)))
}
resulting.matrix
theta <- seq(0.001, 0.5, 0.001)
lambda1 <- runif (1, 2, 5)
lambda2 <- runif (1, 2, 5)
plot(ylim=c(0.9,1),theta, resulting.matrix[1:500,2], type ="l", ylab = "coverage probability")
abline(h = 0.95, lty = 2)
Reply | Threaded
Open this post in threaded view
|

Re: Coverage probability for a Poisson parameter

JS Huang
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Coverage probability for a Poisson parameter

Angle
This post has NOT been accepted by the mailing list yet.
This post was updated on .
I slightly modified in this way but it seems to not work properly even now. It stops suddenly at theta=0.395 giving the same error as before. The error is the same as before, and depending on the times, it just stops working randomly, whatever the values of theta is. It just happens randomly and I don't understand why, since there shouldn't be any missing value.

cover <- function(theta, lambda1, lambda2, significance.level)  {
  s1 <- rpois(1,lambda1)
  s2 <- rpois(1,lambda2)
  theta <- lambda2/(lambda1+lambda2)
  s <- s1+s2
  z <- qnorm(1-0.05/2)
  k <- z^2
  pi <- s2/s
  root <- ((s2/s)*(1-s2/s)+k/(4*s))^(1/2)
  low <- (s2+k/2)/(s+k)-((z*sqrt(s))/(s+k))*root
  hig <- (s2+k/2)/(s+k)+((z*sqrt(s))/(s+k))*root
  if (theta >= low & theta <= hig){1} else {0}
}
resulting.matrix <- matrix(0, nrow=500,ncol=2)
for (i in 1:500) {
  resulting.matrix[i,1] <- 0.001 * i
  resulting.matrix[i,2] <- mean(sapply(1:100,function(x) cover(0.001*i,lambda1,lambda2,0.05)))
}
resulting.matrix
lambda1<-runif(1,0.5,5)
lambda2<-runif(1,0.5,5)
theta <- seq(0.001, 0.500, 0.001)
plot(ylim=c(0.9,1),theta, resulting.matrix[1:500,2], type ="l", ylab = "coverage probability")
abline(h = 0.95, lty = 2)
Reply | Threaded
Open this post in threaded view
|

Re: Coverage probability for a Poisson parameter

JS Huang
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Coverage probability for a Poisson parameter

Angle
This post has NOT been accepted by the mailing list yet.
You are right, such a stupid mistake. Thank you, you're helping me a lot with my initiation to R :)
Reply | Threaded
Open this post in threaded view
|

Re: Coverage probability for a Poisson parameter

JS Huang
In reply to this post by JS Huang
CONTENTS DELETED
The author has deleted this message.