創建的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