Массовая функция биномиальной вероятности с доверительным интервалом - PullRequest
Новая
       73

Массовая функция биномиальной вероятности с доверительным интервалом

0 голосов
/ 19 декабря 2018

Следующий вопрос нам нужно решить.

Рассмотрим следующее binomial probability mas function (pmf):

f(x;m,p) = (m¦x) p^x * (1-p)^(m-x), для x = 0, 1, 2,.....,m, и в остальном равно 0.Пусть X_1, X_2,....,Xn - это независимые и одинаково распределенные случайные выборки из f(x;m = 20; p = 0:45).

1). Предположим, что n = 15, и вычислим 95-процентный доверительный интервал для p, используя p-hat = Σ_(i=1)^n X_i/mn (оценка p).Смоделируйте эти доверительные интервалы 10000 раз и посчитайте, как часто значение параметра p лежит в этих 10000 доверительных интервалах. 

m <- 20
p <- 0.45
n <- 15
x <- m
nsim <- 10000
counter <- 0

for (i in 1:nsim) {
  bpmf <- rbinom(x,m,p)
  esti_p <- bpmf/(m*n)
  var_bpmf <- var(bpmf) 
  CI_lower <- esti_p - qnorm(0.975)*sqrt(var_bpmf/n) 
  CI_upper <- esti_p + qnorm(0.975)*sqrt(var_bpmf/n) 
  if ((CI_lower<p) & (CI_upper>p)) counter <- counter + 1
}

Это не работает должным образом, и я не вижу, что делаюнеправильно.Кто-нибудь может мне помочь с этим?

Когда я запускаю свой код, я верю, что ответ сейчас правильный, но он дает следующее предложение: «Было 50 или более предупреждений (используйте предупреждения(), чтобы увидеть первые 50) «Когда я запускаю это, это даст: 

"1: In if ((CI_lower < p) & (CI_upper > p)) counter <- counter +  ... :
the condition has length > 1 and only the first element will be used".

Также я не знаю точно, если; 

 CI_lower <- esti_p - qnorm(0.975)*sqrt(var_bpmf/n) 
 CI_upper <- esti_p + qnorm(0.975)*sqrt(var_bpmf/n)

это правильная формуларассчитать доверительный интервал.

1 Ответ

0 голосов
/ 23 декабря 2018 
m <- 20
p <- 0.45
nsim <- 10000

  bpmf <- rbinom(size=m,prob=p,n=nsim)
  esti_p <- bpmf/m
  var_bpmf <- esti_p*(1-esti_p)/m 
  CI_lower <- esti_p - qnorm(0.975)*sqrt(var_bpmf) 
  CI_upper <- esti_p + qnorm(0.975)*sqrt(var_bpmf) 
  counter <-((CI_lower<p) & (CI_upper>p))
table(counter)
...