Создание функции 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
...

где h колеблется от 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