创建的R函数从多核包

Question

我需要分析一些模拟数据具有以下结构mclapply使用：

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