Создание функции R, чтобы использовать mclapply из пакета многоядерной

Question

мне нужно проанализировать некоторые смоделированные данные со следующей структурой:

h c x1    y1    x1c10 
1 0 37.607056431 104.83097593 5 
1 1 27.615251557 140.85532974 10 
1 0 34.68915314  114.59312842 2 
1 1 30.090387454 131.60485642 9 
1 1 39.274429397 106.76042522 10 
1 0 33.839385007 122.73681319 2 
... 

где Н находится в диапазоне от 1 до 2500, а индексы образца методом Монте-Карло, каждый из образцов с 1000 наблюдений. Я анализировать эти данные с помощью следующего кода, который дает мне два объекта (fnN1, fdQB101):

mc<-2500 ##create loop index 
fdN1<-matrix(0,mc,1000) 
fnQB101 <- matrix(0,mc,1000) ##create 2500x1000 storage matrices, elements zero 

for(j in 1:mc){ 

fdN1[j,] <- dnorm(residuals(lm(x1 ~ c,data=s[s$h==j,])), 
        mean(residuals(lm(x1 ~ c,data=s[s$h==j,]))), 
          sd(residuals(lm(x1 ~ c,data=s[s$h==j,])))) 

x1c10<-as.matrix(subset(s,s$h==j,select=x1c10)) 

fdQB100 <- as.matrix(predict(polr(as.factor(x1c10) ~ c , 
            method="logistic", data=s[s$h==j,]), 
             type="probs")) 

indx10<- as.matrix(cbind(as.vector(seq(1:nrow(fdQB100))),x1c10)) 

fdQB101[j,] <- fdQB100[indx10] 

} 

Объекты fdN1 и fdQB101 являются 2500x1000 матрицы с предсказанными вероятностями в качестве элементов. Мне нужно создать функцию из этого цикла, которую я могу вызвать с помощью lapply() или mclapply(). Когда я обернуть это в следующей команде функции:

ndMC <- function(mc){ 

for(j in 1:mc){ 
... 
} 
return(list(fdN1,fdQB101)) 

} 
lapply(mc,ndMC) 

объекты fdN1 и fdQB101 каждый возвращается в виде матриц 2500x1000 нулей вместо предсказанных вероятностей. Что я делаю не так?

Answer 1

Вы должны быть в состоянии сделать это с помощью пакета data.table. Вот пример:

library(data.table) 
dt<-data.table(h=rep(1L,6), c=c(0L,1L,0L,1L,1L,0L), 
      X1=c(37.607056431,27.615251557,34.68915314,30.090387454,39.274429397,33.839385007), 
      y1=c(104.83097593,140.85532974,114.59312842,131.60485642,106.76042522,122.73681319), 
      x1c10=c(5L,10L,2L,9L,10L,2L)) 

## Create a linear model for every grouping of variable h: 
fdN1.partial<-dt[,list(lm=list(lm(X1~c))),by="h"] 

## Retrieve the linear model for h==1: 
fdN1.partial[h==1,lm] 
## [[1]] 
## 
## Call: 
## lm(formula = X1 ~ c) 
## 
## Coefficients: 
## (Intercept)   c 
##  35.379  -3.052 

Вы также могли бы написать функцию, чтобы обобщить это решение:

f.dnorm<-function(y,x) { 
    f<-lm(y ~ x) 
    out<-list(dnorm(residuals(f), mean(residuals(f)), sd(residuals(f)))) 
    return(out) 
} 

## Generate two dnorm lists for every grouping of variable h: 
dt.lm<-dt[,list(dnormX11=list(f.dnorm(X1,rep(1,length(X1)))), dnormX1c=list(f.dnorm(X1,c))),by="h"] 

## Retrieve one of the dnorm lists for h==1: 
unlist(dt.lm[h==1,dnormX11]) 
##   1   2   3   4   5   6 
## 0.06296194 0.03327407 0.08884549 0.06286739 0.04248756 0.09045784