использование списков для моделирования - PullRequest
2 голосов
/ 10 июня 2011

Я поставил перед собой небольшую проблему на пути к обучению R. Вопрос был, учитывая выборку из 500 чисел в нормальном распределении со средним значением 20, сколько чисел до 20 я получу для стандартных отклонений от 6 до 10. Просто Чтобы узнать больше, я решил получить 4 образца для каждого SD. Итак, к концу у меня должно быть:

sd6samp1: ...

sd6samp2: ...

....

sd10samp4: ...

Мой первый подход, который сработал, был:

 ddss<-c(6:10) # sd's
 sam<-c(1:4) # 4 samples for each
 k=0  # counter in 0
 for (i in ddss) {   # for each sd
   for (j in sam) {  # for each sample
     nam <- paste("sam",i,".",j, sep="") # building a name
     n <- assign(nam,rnorm(500, 20, i))  # the great assign function
     k <- k+sum(n<=0)
   }
   print(assign(paste("ds",i,sep=""), k)) # ohh assign you're great
   k=0 # reset counter
 }

Ища, как создавать имена переменных с помощью цикла 'i', обнаружил, что 'assign' выполняет свою работу, но он также сказал:

Обратите внимание, что если вы планируете некоторые симуляции, многие guRus скажут, что вы должны использовать список.

Так что я подумал, что было бы хорошо выучить списки ...

В то же время я также обнаружил отличный вариант ... ддсс <- с (6:10) </p>

for (i in ddss) {
   print(paste('prob. x<=0), with sd=',i))
   print(pnorm(0,mean=20,sd=i)*500)
}

Это сработало, чтобы ответить на вопрос, но списки еще предстояло составить ... и многое из R еще предстоит изучить. Основная идея состояла не в том, чтобы узнать саму вероятность или количество негативов ... а в том, чтобы узнать R и, в частности, некоторые циклы.

Итак, я пытался использовать упомянутые списки

Мой самый близкий подход был:

ddss<-c(6:10) # sd's to be calculated.
sam<-c(1:4) # 4 samples for each sd
liss<-list()  # initializing the list
for (i in ddss) {   # for each sd
   liss[[i]] <- list()
   for (j in sam) {  # for each sample
      liss[[i]][[j]] <- rnorm(500, 20, i)
      print(paste('ds',i,'samp',j,'=',sum(liss[[i]][[j]]<0)))
   }
}

С этим я получаю информацию, но меня интересуют два вопроса (1 и 2) и некоторые другие вопросы (3 и 4):

  1. Я получаю список из 10 элементов, 6 пустых и затем 4 с подсписками. Кажется, я не могу понять, как работать с элементами 1: 4 списка (sd) с именами 6: 9 (сами sd).

  2. Даже при том, что я пытался, я не мог назвать элементы списков через циклы 'for'. Любое понимание этих вопросов было бы замечательно.

  3. Так как в этом контексте симуляции. Как вы думаете, что лучше: вложенные списки (списки с подсписками) или простые (более длинные) списки?

  4. Мне было интересно, будут ли здесь полезны функции apply, я пытался что-то сделать, например:

vbv<-matrix(c(6,6,6,6,7,7,7,7,8,8,8,8,9,9,9,9))
lsl<-apply(vbv, 2, function(x) rnorm(500,20,x))

Но, похоже, я даже близко не подхожу ...

Спасибо за ваше время, если вы прочитали это далеко!

Вы можете также взять еще немного, чтобы ответить; -).

Ответы [ 3 ]

4 голосов
/ 10 июня 2011

Проблема в ваших индексах: вы работаете над индексатором i из ddss, который работает с 6 до 10. Таким образом, в первом рабочем цикле во внешнем цикле ваше первое утверждение действительно говорит: liss[[6]]<-list(), подразумевая, чтопервые 5 из них равны NULL.

Так что, если вы настаиваете на работе с циклами, вам следует сделать следующее (отметьте ?seq_along):

ddss<-c(6:10) # sd's to be calculated.
sam<-c(1:4) # 4 samples for each sd
liss<-list()  # initializing the list
for (i in seq_along(ddss)) {   # now, i runs from 1 to 5
   liss[[i]] <- list()
   for (j in sam) {  # for each sample
      liss[[i]][[j]] <- rnorm(500, 20, i)
      print(paste('ds',ddss[i],'samp',j,'=',sum(liss[[i]][[j]]<0)))
   }
   names(liss[[i]])<-as.character(sam)#this should solve your naming issue (1/2)
}
names(liss)<-as.character(ddss)#this should solve your naming issue (2/2)

Обратите внимание, что, как всегда,было бы неплохо назвать ваши переменные чем-то более полезным, чем i или j: если бы вы назвали его curd, может быть, вы бы не использовали его немедленно в качестве индексатора в списке?

Теперь, есливы действительно стремитесь к улучшению (но хотите придерживаться списков), вы действительно хотите использовать функции стиля применения:

liss<-lapply(ddss, function(curds){ #apply the inline function to each ds and store results in a list
  return(lapply(sam, function(cursam){ #apply inline function to each sam and store results in a list
    rv<-rnorm(500, 20, curds)
    cat('ds',curds,'samp',cursam,'=',sum(rv<0), "\n") #maybe better for your purposes.
    return(rv)
  }))
}) 

Наконец, для вашего случая не так уж много причиниспользуйте списки (и при этом вам даже не нужно сохранять выборочные данные для каждого ds / sam): вы можете хранить все как трехмерный массив, но так как вы указываете это как учебное упражнение (эй, возможно, использование массива может быть вашим следующим упражнением:-)), я оставлю это на этом.

3 голосов
/ 10 июня 2011

Я собираюсь добавить другое решение, используя пакет plyr, который, я думаю, специально создан для таких упражнений.

library(plyr)

# generate a data frame of parameters, repeating some as required
parameters  = data.frame(mean = 20, sd = rep(6:10, each = 4))

# generate sample data for each combination of parameters
sample_data = mdply(df, rnorm, n = 500)

# generate answer by counting number of observations less than 20
answer = data.frame(
    parameters, 
    obs_less_20 = rowSums(sample_data[,-c(1, 2),] < 20)
)

head(answer)

mean sd obs_less_20
1   20  6         247
2   20  6         250
3   20  6         242
4   20  6         259
5   20  7         240
6   20  7         237
3 голосов
/ 10 июня 2011

lapply() полезен здесь, где мы можем просто применить к набору значений для SD.Это помогает написать пользовательскую оболочку вокруг функции rnorm(), чтобы мы могли передавать различные значения для различных аргументов rnorm() и обрабатывать реплики k ( k = 4в твоем примере) тоже неплохо.Эта оболочка foo() ниже:

foo <- function(sd, n, mean, reps = 1) {
    rands <- rnorm(n * reps, mean = mean, sd = sd)
    if(reps > 1)
        rands <- matrix(rands, ncol = reps)
    rands
}

Мы используем ее в вызове lapply() следующим образом:

sims <- lapply(6:10, FUN = foo, mean = 20, n = 500, reps = 4)

, который дает:

R> str(sims)
List of 5
 $ : num [1:500, 1:4] 30.3 22 15.6 20 19.4 ...
 $ : num [1:500, 1:4] 20.9 21.7 17.7 35 30 ...
 $ : num [1:500, 1:4] 17.88 26.48 5.19 19.25 15.59 ...
 $ : num [1:500, 1:4] 27.41 12.72 9.38 35.09 11.08 ...
 $ : num [1:500, 1:4] 16.2 11.6 20.5 35.4 27.3 ...

Затем мы можем вычислить количество наблюдений <20 на SD </p>

names(sims) <- paste("SD", 6:10, sep = "")
out <- lapply(sims, function(x) colSums(x < 20))

, что дает:

R> out
$SD6
[1] 218 251 253 227

$SD7
[1] 250 242 233 232

$SD8
[1] 258 241 246 274

$SD9
[1] 252 245 249 258

$SD10
[1] 253 259 241 242

@ Джорис предлагает мне показать, как получить доступ к элементам списка.Например, если вы хотите получить результаты моделирования для SD = 20, мы могли бы сделать out[[4]], потому что 20 было четвертым значением в векторе SD, к которому мы применили, или потому что я назвал элементы списка вывода out, мы можем в качестве результатов моделирования использовать out[["SD10"]].

Чтобы ответить на некоторые конкретные вопросы о ваших циклах и т. Д.,

  • , чтобы добавить имена киспользование списка names(), например names(mylist) <- c ("foo", "bar") <code>. You'd be better off in your loop calling names () `один раз за итерацию цикла для установки имен в одном кадре - вы, вероятно, не захотитеЯ хочу заполнить имена по мере продвижения, поскольку это было бы неэффективно.
  • Я не думаю, что это имеет слишком большое значение, если вы используете вложенный список или список, содержащий матрицу, как в моем примере.Чтобы изменить foo() для возврата списка, чтобы вывод lapply() представлял собой список списков, мы могли бы сделать:

Код:

bar <- function(sd, n, mean, reps = 1) {
    rands <- rnorm(n * reps, mean = mean, sd = sd)
    if(reps > 1)
        rands <- split(rands, rep(seq_len(reps), each = n))
    rands
}
sims2 <- lapply(6:10, FUN = bar, mean = 20, n = 500, reps = 4)
names(sims2) <- paste("SD", 6:10, sep = "")
out2 <- lapply(sims2, function(x) sapply(x, function(y) sum(y < 20)))

, который дает то же самоевывод как и раньше.

...