Следующая функция использует vapply и массивы высокой размерности для вычисления некоторых значений p независимости для заданных полей в наборе данных. Проблема в том, что это так медленно ... Я даю вам dput некоторых примеров данных и саму функцию. Если вы найдете что-нибудь, что я мог бы сделать, чтобы сделать это быстрее, пожалуйста, скажите мне;)
exemple_data = list(z = structure(c(0.627450980392157, 0.372549019607843, 0.637254901960784,
0.843137254901961, 0.450980392156863, 0.686274509803922, 0.0980392156862745,
0.901960784313726, 0.686274509803922, 0.666666666666667, 0.823529411764706,
0.137254901960784, 0.92156862745098, 0.784313725490196, 0.647058823529412,
0.235294117647059, 0.647058823529412, 0.431372549019608, 0.0980392156862745,
0.0196078431372549, 0.803921568627451, 0.666666666666667, 0.0980392156862745,
0.509803921568627, 0.764705882352941, 0.372549019607843, 0.431372549019608,
0.598039215686274, 0.941176470588235, 0.333333333333333, 0.0196078431372549,
0.549019607843137, 0.392156862745098, 0.490196078431373, 0.392156862745098,
0.666666666666667, 0.764705882352941, 0.0784313725490196, 0.274509803921569,
0.862745098039216, 0.196078431372549, 0.92156862745098, 0.215686274509804,
0.235294117647059, 0.529411764705882, 0.549019607843137, 0.980392156862745,
0.274509803921569, 0.392156862745098, 0.431372549019608, 0.941176470588235,
0.156862745098039, 0.862745098039216, 0.901960784313726, 0.686274509803922,
0.117647058823529, 0.882352941176471, 0.294117647058824, 0.254901960784314,
0.303921568627451, 0.588235294117647, 0.313725490196078, 0.529411764705882,
0.705882352941177, 0.705882352941177, 0.745098039215686, 0.215686274509804,
0.980392156862745, 0.803921568627451, 0.745098039215686, 0.725490196078431,
0.0784313725490196, 0.745098039215686, 0.882352941176471, 0.568627450980392,
0.519607843137255, 0.235294117647059, 0.686274509803922, 0.588235294117647,
0.843137254901961, 0.0980392156862745, 0.823529411764706, 0.117647058823529,
0.294117647058824, 0.137254901960784, 0.313725490196078, 0.96078431372549,
0.0392156862745098, 0.156862745098039, 0.549019607843137, 0.0392156862745098,
0.529411764705882, 0.666666666666667, 0.764705882352941, 0.0980392156862745,
0.411764705882353, 0.607843137254902, 0.176470588235294, 0.0196078431372549,
0.156862745098039, 0.607843137254902, 0.411764705882353, 0.882352941176471,
0.607843137254902, 0.490196078431373, 0.882352941176471, 0.176470588235294,
0.784313725490196, 0.647058823529412, 0.588235294117647, 0.980392156862745,
0.274509803921569, 0.470588235294118, 0.627450980392157, 0.941176470588235,
0.137254901960784, 0.627450980392157, 0.156862745098039, 0.117647058823529,
0.823529411764706, 0.490196078431373, 0.0392156862745098, 0.833333333333333,
0.862745098039216, 0.176470588235294, 0.92156862745098, 0.490196078431373,
0.549019607843137, 0.450980392156863, 0.92156862745098, 0.470588235294118,
0.254901960784314, 0.803921568627451, 0.823529411764706, 0.627450980392157,
0.901960784313726, 0.196078431372549, 0.725490196078431, 0.725490196078431,
0.901960784313726, 0.519607843137255, 0.509803921568627, 0.705882352941177,
0.666666666666667, 0.196078431372549, 0.274509803921569, 0.784313725490196,
0.343137254901961, 0.313725490196078, 0.372549019607843, 0.156862745098039,
0.705882352941177, 0.313725490196078, 0.411764705882353, 0.607843137254902,
0.0588235294117647, 0.588235294117647, 0.196078431372549, 0.137254901960784,
0.0784313725490196, 0.764705882352941, 0.745098039215686, 0.372549019607843,
0.372549019607843, 0.0588235294117647, 0.784313725490196, 0.862745098039216,
0.254901960784314, 0.0784313725490196, 0.245098039215686, 0.705882352941177,
0.352941176470588, 0.598039215686274, 0.431372549019608, 0.882352941176471,
0.568627450980392, 0.470588235294118, 0.509803921568627, 0.470588235294118,
0.303921568627451, 0.843137254901961, 0.450980392156863, 0.411764705882353,
0.215686274509804, 0.245098039215686, 0.647058823529412, 0.294117647058824,
0.637254901960784, 0.549019607843137, 0.725490196078431, 0.254901960784314,
0.0196078431372549, 0.96078431372549, 0.96078431372549, 0.509803921568627,
0.862745098039216, 0.117647058823529, 0.833333333333333, 0.92156862745098,
0.411764705882353, 0.215686274509804, 0.686274509803922, 0.235294117647059,
0.352941176470588, 0.470588235294118, 0.0784313725490196, 0.843137254901961,
0.343137254901961, 0.196078431372549, 0.117647058823529, 0.352941176470588,
0.0588235294117647, 0.941176470588235, 0.745098039215686, 0.274509803921569,
0.294117647058824, 0.392156862745098, 0.764705882352941, 0.980392156862745,
0.352941176470588, 0.431372549019608, 0.901960784313726, 0.137254901960784,
0.568627450980392, 0.0392156862745098, 0.96078431372549, 0.803921568627451,
0.0196078431372549, 0.0588235294117647, 0.803921568627451, 0.333333333333333,
0.568627450980392, 0.450980392156863, 0.333333333333333, 0.96078431372549,
0.450980392156863, 0.333333333333333, 0.568627450980392, 0.529411764705882,
0.215686274509804, 0.392156862745098, 0.725490196078431, 0.490196078431373,
0.0392156862745098, 0.980392156862745, 0.176470588235294, 0.941176470588235,
0.0588235294117647, 0.176470588235294, 0.784313725490196), .Dim = c(5L,
50L), .Dimnames = list(c("sr", "pop15", "pop75", "dpi", "ddpi"
), c("Australia", "Austria", "Belgium", "Bolivia", "Brazil",
"Canada", "Chile", "China", "Colombia", "Costa Rica", "Denmark",
"Ecuador", "Finland", "France", "Germany", "Greece", "Guatamala",
"Honduras", "Iceland", "India", "Ireland", "Italy", "Japan",
"Korea", "Luxembourg", "Malta", "Norway", "Netherlands", "New Zealand",
"Nicaragua", "Panama", "Paraguay", "Peru", "Philippines", "Portugal",
"South Africa", "South Rhodesia", "Spain", "Sweden", "Switzerland",
"Turkey", "Tunisia", "United Kingdom", "United States", "Venezuela",
"Zambia", "Jamaica", "Uruguay", "Libya", "Malaysia"))), bp = c(0.359829289416585,
0.565182697194852, 0.451504788384914, 0.411767914767176, 0.607858899990815
), binary_repr = structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0,
0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1,
0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0,
0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0,
1, 0), .Dim = c(32L, 5L)), n = 50L, d = 5L)
compute_independence_p_values <- function(z,bp,binary_repr,n,d,N = 13){
# prerequisites :
d = nrow(z)
n = ncol(z)
min = bp*t(binary_repr)
max = bp^t(1-binary_repr)
lambda_l = apply(max-min,2,prod) # 2^d
lambda_k = vapply(1:d,function(d_rem){lambda_l/((max-min)[d_rem,])},lambda_l) # 2^d * d
# first, we comput the empirical value of the statistic :
core = vapply(1:2^d,function(.x){((min[,.x]<=z)*(max[,.x]>z))==1},FUN.VALUE = z) # dims : d, n, 2^d
f_l = apply(apply(core,2:3,all),2,mean) # 2^d
f_k = vapply(1:d,function(d_rem){ apply(apply(core[-d_rem,,],2:3,all),2,mean)},f_l) # 2^d, d
statistic <- sum(f_l^2/lambda_l) -2 * colSums(f_k * f_l / lambda_k) # d
# then we bootstrap it :
z_rep = vapply(1:N,function(i){z},z) # d, n, N
z_repeats = vapply(1:d,function(i){
z = z_rep
z[i,,] = runif(n*N)
z
},z_rep) # d, n, N, D=d
cores = vapply(1:d,function(d_rem){
vapply(1:2^d,function(.x){
((min[,.x]<=z_repeats[,,,d_rem])*(max[,.x]>z_repeats[,,,d_rem]))==1
}, FUN.VALUE = z_repeats[,,,d_rem])
}, FUN.VALUE = array(0,c(d,n,N,2^d))) # d, n, N, 2^d, d
f_l = vapply(1:d,function(d_rem){apply(apply(cores[ ,,,,d_rem], 2:4, all), 2:3, mean)},array(0.,c(N,2^d))) #(N,2^d,d)
f_k = vapply(1:d,function(d_rem){apply(apply(cores[-d_rem,,,,d_rem], 2:4, all), 2:3, mean)},array(0.,c(N,2^d))) #(N,2^d,d)
samples = apply(aperm(f_l^2,c(2,3,1))/lambda_l - 2 * aperm(f_k*f_l,c(2,3,1))/vapply(1:N,function(i){lambda_k},lambda_k),c(2,3),sum)
p_val = rowMeans(statistic < samples)
return(p_val)
}
result = compute_independence_p_values(z = exemple_data$z,
bp = exemple_data$bp,
binary_repr = exemple_data$binary_repr,
n = exemple_data$n,
d = exemple_data$d,
N=13) # <<<< The number of bootstrap resamples. I need this to be like 10^3 or 10^4
Код широко использует функцию vapply для выполнения простых вычислений над массивами. Но это так медленно. Можете ли вы придумать лучший способ?