2013-07-19 5 views
3

Мне нужно поместить код в кусочную функцию и сохранить сгенерированные значения в кадре данных. Правила таковы:Как закодировать кусочную функцию в R для небольшого моделирования и сохранения значений в кадре данных

  • У меня есть объект X, который генерируется Бернулли (1/3).
  • Если X = 0, другой объект, Y, генерируется E = Экспоненциальный (1).
  • Если X = 1, Y генерируется E, если E < = P и (P + EL), если E> P, где P - константа (1, например) и EL = экспоненциальная (Lambda) независимая от E.
  • Я хочу создать кадр данных с 100 образцами X и Y, полученными с использованием этого метода, и дополнительно выполнить этот процесс 10000 раз (или, другими словами, сгенерировать 10000 кадров данных).

Я попытался сделать что-то вроде этого, но поскольку я полный новичок с помощью R, я не мог понять, как правильно определить каждый элемент, и, очевидно, как сохранить результаты в кадре данных ,

Вот «код» я сделал:

test <- function(lambda,p) { 
for (i in 1:10000) { # Number of simulations that I want to do. 
for(j in 1:100) { # Sample size of each simulation/data frame. 
    x <- rbinom(1,1,1/3) 
    e <- rexp(1,1) 
    if (x==0) {y <- e} 
    else { 
    if (e<=p) {y <- e} 
    else {y <- p + rexp(1, lambda)} 
    } 
    } 
}  

Но еще до тестирования я знаю, что это невозможно, что она работает должным образом. Кадр данных, который я хочу сделать, может содержать только значения для X и Y.

Я знаю, что это может быть очень простой вопрос, поэтому большое спасибо за ваши ответы.

ответ

3

Я думаю:

simfun <- function(n=100,lambda=1,P=1,ret.df=TRUE) { 
    X <- rbinom(n,prob=1/3,size=1) 
    E <- rexp(n,1) 
    EL <- rexp(n,lambda) 
    Y <- ifelse(X==0,E, 
       ifelse(E<=P,E,P+E*EL)) 
    ## return data frame or matrix 
    if (ret.df) data.frame(X,Y) else cbind(X,Y) 
} 
system.time(res <- replicate(1e4,simfun(),simplify=FALSE)) 
## 6 seconds 

## or: 
library(plyr) 
system.time(res2 <- raply(1e4,simfun(ret.df=FALSE),.progress="text") 
## returns a 1e4 x 100 x 2 array, maybe more convenient: use 
## rlply instead of raply (and ret.df=TRUE) to match the previous 
## results 
+1

+1! Я тоже думаю :) 'ifelse' +' replicate'! – agstudy

+0

+1. Большое спасибо Вам! Но как насчет того, как создается каждый кадр данных? Я имею в виду ... Если мне нужно применить тест для каждого фрейма данных, каково будет имя или индикатор для каждого из них? И еще. Если я хочу сделать симуляции Z, мне нужно только поставить эту функцию между '{}' в аргументе 'for (i in 1: Z)'? Как насчет «оценки» каждого генерируемого кадра данных? – anxoestevez

+1

вы можете 'lapply()' применять тест для каждого фрейма данных или индексировать их как 'res [[i]]' где 'i' идет от 1 до 1e4 ... как для симуляции Z - вам нужно еще * другой * уровень петли ??? В приведенном выше коде уже генерируется 10 000 фреймов данных ... –