Re: [R] Coverage probability for a Poisson parameter
Hi,
In your function cover, lambda1 and lambda2 are used but not in the
argument of the function. I suppose that you need to have lambda1 and
lambda2 in the argument of the function cover, like function(lambda1,
lambda2, n, significance.level).
Give it a try.
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}
}
--
View this message in context:
http://r.789695.n4.nabble.com/Coverage-probability-for-a-Poisson-parameter-tp4702535p4703230.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Coverage probability for a Poisson parameter
Hi,
Given the function cover, it's very likely that you will get 0 for both s1
and s1 with small value of lambda1 and lambda2. In that case the sum s will
be 0. With s being 0, you will have issue with the expression in pi -
s2/s and root - ((s2/s)*(1-s2/s)+k/(4*s))^(1/2). You need to take care of
the case that s is 0 before proceeding calculating pi and root.
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}
}
--
View this message in context:
http://r.789695.n4.nabble.com/Coverage-probability-for-a-Poisson-parameter-tp4702535p4703238.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Coverage probability for a Poisson parameter
Hi,
Some suggestion about the arguments of the function defined below. Since
theta is calculated with the value of lambda1 and lambda2, there is no need
to include theta in the argument. Or, your function can be defined as
function(lambda1, lambda2, significance.level)
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}
}
--
View this message in context:
http://r.789695.n4.nabble.com/Coverage-probability-for-a-Poisson-parameter-tp4702535p4703248.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Coverage probability for a Poisson parameter
Hi,
After some thought, I found the treatment of sample mean equal 0 was not
appropriate. I modified the function likelihood.ratio.test.Poisson.
resulting.matrix now has 0.0512 as the average of type I error.
function(lambda, sample.size, significance.level)
{
reject - 0
sample.mean - mean(rpois(sample.size, lambda))
if (sample.mean == 0)
{
test.statistics - 2 * sample.size * lambda
}
else
{
test.statistics - 2 * sample.size * (lambda - sample.mean + sample.mean
* log(sample.mean / lambda))
}
if (test.statistics = qchisq(1 - significance.level, 1)) {reject - 1}
else {reject - 0}
return(reject)
}
for (i in 1:500){
+ resulting.matrix[i,1] - 0.01 * i
+ resulting.matrix[i,2] - mean(sapply(1:100,function(x)
likelihood.ratio.test.Poisson(0.01*i,10,0.05)))
+ }
mean(resulting.matrix[,2])
[1] 0.05102
--
View this message in context:
http://r.789695.n4.nabble.com/Coverage-probability-for-a-Poisson-parameter-tp4702535p4702909.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Coverage probability for a Poisson parameter
Hi,
Roughly reading the code, I find this statement phat - x / m is
probably incorrect since this will give you the set of 100 observed x
values /100. I redefine the function cover with three inputs: lambda
for the parameter of the poisson distribution, sample.size and
significance.level. The output is 1 or 0, depending on whether lambda is
inside the confidence interval or not. With 5% level of significance I
expect to get 95% of the time the parameter will be included in the
confidence interval. The two runs below shows 96% and 94.8%, pretty close.
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:100, function(x)cover(2.5,100,0.05)))
[1] 0.96
mean(sapply(1:1000, function(x)cover(2.5,10,0.05)))
[1] 0.948
--
View this message in context:
http://r.789695.n4.nabble.com/Coverage-probability-for-a-Poisson-parameter-tp4702535p4702537.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Coverage probability for a Poisson parameter
Hi,
My first question is what is your n when you say fixed n. I assume the
lambda is the mean of the poisson distribution that you want to take sample
from.
Another question is about the sample size. It does not make too much
sense to make a sample of size 1.
Let's assume that you want to fix the sample size to be 100 and change
lambda from 0.1 to 5 with an increment of 0.1. For each lambda, plan to
run, say, 1000 times. Then the following will be my approach. Recall that
the function cover returns 1 when lambda is in the confidence interval and 0
otherwise. resulting_matrix is created with size 50 x 2 with 0 populated.
The matrix is to store lambda and the proportion of samples with lambda
inside the confidence interval calculated from samples. With the resulting
matrix, one can see that lambdas are in the first column with values of 0.1
to 5 with increment of 0.1. The corresponding proportions are in the second
column. All of the proportions are from 0.917 to 0.969 as the last line
shows.
Hope this helps.
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}
+ }
resulting.matrix - matrix(0, nrow=50,ncol=2)
for (i in 1:50)
+ {
+ resulting.matrix[i,1] - 0.1 * i
+ resulting.matrix[i,2] - mean(sapply(1:1000,function(x)
cover(0.1*i,100,0.05)))
+ }
resulting.matrix
[,1] [,2]
[1,] 0.1 0.917
[2,] 0.2 0.949
[3,] 0.3 0.928
[4,] 0.4 0.939
[5,] 0.5 0.943
[6,] 0.6 0.949
[7,] 0.7 0.942
[8,] 0.8 0.939
[9,] 0.9 0.945
[10,] 1.0 0.943
[11,] 1.1 0.962
[12,] 1.2 0.933
[13,] 1.3 0.947
[14,] 1.4 0.951
[15,] 1.5 0.946
[16,] 1.6 0.939
[17,] 1.7 0.946
[18,] 1.8 0.953
[19,] 1.9 0.964
[20,] 2.0 0.943
[21,] 2.1 0.937
[22,] 2.2 0.944
[23,] 2.3 0.945
[24,] 2.4 0.950
[25,] 2.5 0.954
[26,] 2.6 0.946
[27,] 2.7 0.945
[28,] 2.8 0.949
[29,] 2.9 0.956
[30,] 3.0 0.953
[31,] 3.1 0.941
[32,] 3.2 0.949
[33,] 3.3 0.943
[34,] 3.4 0.956
[35,] 3.5 0.950
[36,] 3.6 0.944
[37,] 3.7 0.952
[38,] 3.8 0.958
[39,] 3.9 0.938
[40,] 4.0 0.944
[41,] 4.1 0.950
[42,] 4.2 0.945
[43,] 4.3 0.948
[44,] 4.4 0.962
[45,] 4.5 0.969
[46,] 4.6 0.956
[47,] 4.7 0.950
[48,] 4.8 0.955
[49,] 4.9 0.946
[50,] 5.0 0.945
range(resulting.matrix[,2])
[1] 0.917 0.969
--
View this message in context:
http://r.789695.n4.nabble.com/Coverage-probability-for-a-Poisson-parameter-tp4702535p4702551.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
