Исправление функции на симуляции Бернулли - PullRequest
Новая
       25

Исправление функции на симуляции Бернулли

0 голосов
/ 28 января 2019

Я хотел бы создать функцию с тремя параметрами: размер выборки, количество выборок, истинное значение p успеха в исследовании Бернулли.

Эта функция выдаст следующие результаты: средняя оценкаp (т. е. среднее значение p оценок для каждой выборки), средняя оценка истинного стандартного отклонения (т. е. среднее значение sd оценок для каждой выборки) и, наконец, доля выборок, в которых 95% -й ДИ содержалtrue p.

Я придумываю следующий код: 

   draw_function <- function(si=1, proba= 0.5){
   sample_draw <- rbinom(n= si, 1, prob= proba)
   mean_estimate <-  mean(sample_draw)
   s_estimate <- sqrt(mean_estimate*(1-mean_estimate))
  confidence_interval <- list(c(mean_estimate - 
  1.96*s_estimate/sqrt(si), mean_estimate + 1.96*s_estimate/sqrt(si)))
    list(mean_sample = mean_estimate, s_sample = s_estimate, c_sample = 
  confidence_interval)
   }



   draw_sample <- function(s=1, r=1, p=0.5) {
   for (i in 1:r) {
     inter <-  list(mean_p=c(), mean_s = c(), mean_c = c() )
     samp <- draw_function(si=s, proba= p)
    inter$mean_p = c(inter$mean_p, samp$mean_sample)
  inter$mean_s = c(inter$mean_s, samp$s_sample)
  inter$mean_c <- c(inter$mean_c, samp$c_sample)
    if ( inter[[3]][1] > p & inter[[3]][2] < p ) {
   mean_c <- (mean_c + 1)/i
   } 
  else {
  meanc <- (mean_c + 0)/i
  }
    }

  return(list(mean(inter$mean_p), mean(inter$mean_s), inter$mean_c))

Однако этот код не работает даже после внесения изменений в ошибки, которые были любезно доведены до моего сведения.

Я получаю эту ошибку:

Ошибка в draw_sample (s = 30, r = 1000, p = 0,05):
(список) объект не может быть приведен к типу'double'

Поэтому я бы хотел попросить вас помочь найти проблему и создать такую ​​функцию!Спасибо!

1 Ответ

0 голосов
/ 28 января 2019

Ошибка возникает в следующем разделе: 

inter[[3]][1] > p & inter[[3]][2] < p

Используя браузер () или построчно запуская код, вы заметите, что 

confidence_interval <- list(c(mean_estimate - 1.96*s_estimate/sqrt(si), 
  mean_estimate + 1.96*s_estimate/sqrt(si)))

возвращает сам список.Так что, если вы хотите взять части сами, это должны быть inter[[3]][[i]][1] и inter[[3]][[i]][2] соответственно.

Возможны несколько оптимизаций, я, скорее всего, отредактирую этот комментарий в последнее время с несколькими предложениями.:-)

::: Edit ::: Несколько более глубокий ответ.Функции r[dens] в R связаны с некоторыми издержками.В качестве такого простого способа ускорения кода можно запустить все симуляции один раз и выполнить вычисления для подмножеств умным способом.Ниже показан один из способов (если имитация находится в пределах памяти), вставка всего в матрицу и последующее использование rowMeans для быстрого вычисления необходимой информации. 

draw_sample <- function(s = 1, r = 1, p = .5){
  samples <- matrix(rbinom(n = r * s, size = 1, p = p), ncol = s, nrow = r, byrow = TRUE)
  mean_p <- rowMeans(samples)
  mean_s <- mean_p * (1 - mean_p)
  bound = 1.96 * mean_s / sqrt(s)
  mean_c <- mean(ifelse(mean_p - bound < p & mean_p + bound > p, 1, 0))
  list(mean_p = mean(mean_p), mean_s = mean(mean_s), mean_c = mean_c)
}
...