Код для моделирования методом Монте-Карло: генерировать образцы заданного размера в R

Question

Я начал путем генерации образца из 500 равномерно распределенных случайных чисел между 0 и 1, используя этот код:

set.seed(1234) 
X<-runif(500, min=0, max=1) 

Теперь, мне нужно напишите psuedocode, который генерирует 10000 выборок N = 500 для моделирования MC, вычисляет среднее значение моего вновь созданного X и сохраняет номер итерации и среднее значение в объекте результата. Я никогда не пытался это, и до сих пор у меня есть это:

n.iter <-(10000*500) 
results <- matrix (0, n.iter, 4) 

Наконец, как только это будет сделано, я запустить его, а затем получить медиану, среднее, и мин/макс начисленной выборочных средних и сохраните их в кадре данных, который называется MC.table. (Также обратите внимание, выше, я понятия не имею, почему в матричном коде есть «4» --- я отработал предыдущие примеры). Любые советы или помощь будут очень признательны.

EDIT: У меня есть пример, который может работать, но я не понимаю, что происходит с ним, поэтому, пожалуйста, более подробной информацией о его действиях для этого:

Ni <- 10000 
n <- 500 
c <- 0 

for (i in n){ 
for (j in 1:Ni){ 
c <- c+ 1 
d <- data.frame (x= , y=) 
results [c,1] <- c 
results [c,2] <- j 
results [c,3] <- i 
results [c,4] <- something(d$x, d$y) 
rm (d) } } 

Если вы могли бы даже взять время чтобы объяснить, что это значит, это поможет мне! Благодаря!

Answer 1

Вы можете попробовать использовать data.table, пакет, который может быть установлен с использованием install.packages("data.table"). С установленным, вы будете запускать что-то вроде ...

> require(data.table) 
> dt <- data.table(x=runif(500*10000),iter=rep(1:500,each=10000)) 
        # x iter 
     # 1: 0.48293196 1 
     # 2: 0.61935416 1 
     # 3: 0.99831614 1 
     # 4: 0.26944687 1 
     # 5: 0.38027524 1 
    # ---     
# 4999996: 0.11314160 500 
# 4999997: 0.07958396 500 
# 4999998: 0.97690312 500 
# 4999999: 0.81670765 500 
# 5000000: 0.62934609 500 
> summaries <- dt[,list(mean=mean(x),median=median(x)),by=iter] 
    # iter  mean median 
    # 1: 1 0.5005310 0.5026592 
    # 2: 2 0.4971551 0.4950034 
    # 3: 3 0.4977677 0.4985360 
    # 4: 4 0.5034727 0.5052344 
    # 5: 5 0.4999848 0.4971214 
# ---       
# 496: 496 0.5013314 0.5048186 
# 497: 497 0.4955447 0.4941715 
# 498: 498 0.4983971 0.4910115 
# 499: 499 0.5000382 0.4997024 
# 500: 500 0.5009614 0.4988237 
> min_o_means <- min(summaries$mean) 
# [1] 0.4914826 

Я думаю, что синтаксис довольно прост. Вы можете посмотреть некоторые из функций, используя ? (например, ?rep). Строки, начинающиеся с #, просто отображают сгенерированные объекты. В data.tables число слева от : - это только номер строки, а --- - строки, пропущенные на дисплее.

Answer 2

Я предполагаю, что ответ, который я дал бы, будет зависеть от того, хотите ли вы научиться псевдокоде или если вы хотите изучить способ «R» ish сделать это. Этот ответ - это то, что я бы рекомендовал для кого-то, кто хотел бы научиться работать с R.

Сначала я сделал бы матрицу с N столбцами и 10000 строк. R ценит это, когда мы делаем пространство раньше времени, чтобы цифры вошли.

X=matrix(NA,nrow=10000,ncol=500)

Вы знаете, как генерировать 500 случайных величин, для одной строки.

runif(500,0,1)

Теперь вы должны выяснить, как получить, что произойдет 10000 раз и назначить каждый X. Возможно, цикл будет работать.

for(i in 1:10000) X[i,]=runif(500,0,1)

Затем вы должны выяснить, как получить резюме каждой строки. Одна из функций, которая может помочь, - rowMeans().Посмотрите на своей странице справки, а затем попытаться получить среднее значение каждой строки табличного X

, чтобы получить средства каждой итерации

rowMeans(X)

затем, чтобы получить представление о том, что эти цифры как Я мог бы быть склонен к

plot(rowMeans(X))

Answer 3

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

test<-function(
    seed=1234, 
    sample.size=500, 
    resample.number=1000, 
    alpha=0.05 
    ) 
    { 

     #initialize original sample 
     original.sample<-runif(sample.size, min=0, max=1) 

     #initialize data.frame 
     resample.results<-data.frame("Run.Number"=NULL,"mean"=NULL) 
     for(counter in 1:resample.number){ 
      temp<-sample(original.sample, size=length(original.sample), replace = TRUE) 
      temp.mean<-mean(temp) 
      temp.table.row<-data.frame("Run.Number"=counter,"mean"=temp.mean) 
      resample.results<-rbind(resample.results,temp.table.row) 
     } 
     resample.results<-resample.results[with(resample.results, order(mean)), ] 

     #for the mean information 
     lowerCI.row<-resample.number*alpha/2 
     upplerCI.row<-resample.number*(1-(alpha/2)) 
     median.row<-resample.number/2 

     #for the mean information 
     median<-resample.results$mean[median.row] 
     lowerCI<-resample.results$mean[lowerCI.row] 
     upperCI<-resample.results$mean[upplerCI.row] 

     #for the position information 
     median.run<-resample.results$Run.Number[median.row] 
     lowerCI.run<-resample.results$Run.Number[lowerCI.row] 
     upperCI.run<-resample.results$Run.Number[upplerCI.row] 

     mc.table<-data.frame("median"=NULL,"lowerCI"=NULL,"upperCI"=NULL) 
     values<-data.frame(median,lowerCI,upperCI) 
     #as.numeric because R doesn't like to mix data types 
     runs<-as.numeric(data.frame(median.run,lowerCI.run,upperCI.run)) 
     mc.table<-rbind(mc.table,values) 
     mc.table<-rbind(mc.table,runs) 

     print(mc.table) 
    } 

После передискретизации данных вы найдете среднее значение. Затем вы заказываете все свои повторно выбранные средства. Середина этого списка является медианной. И, например, с 10000 перевыборами, 250-й заказный средний результат будет вашим нижним 95% -ным CI. Хотя я не сделал этого здесь, значение min будет только в позиции 1, а максимальное значение будет в позиции 10000. Будьте осторожны, когда вы уменьшаете номер повторной выборки: способ, которым я рассчитал позиции, может стать десятичным значением, которое будет путать R.

Кстати, я положил это в форму функции. Если вам нравится делать строки по строкам, просто убедитесь, что все строки связаны между функциями() и следующими основными {}