आर

Question

में मैट्रिक्स संचयी मानक विचलन की कुशल गणना मैंने हाल ही में आर-सहायता मेलिंग सूची पर इस प्रश्न को पोस्ट किया लेकिन मुझे कोई जवाब नहीं मिला, इसलिए मैंने सोचा कि मैं इसे यहां भी पोस्ट करूंगा और देख सकता हूं कि कोई सुझाव था या नहीं।

मैं एक मैट्रिक्स के संचयी मानक विचलन की गणना करने की कोशिश कर रहा हूं। मुझे एक ऐसा फ़ंक्शन चाहिए जो मैट्रिक्स स्वीकार करता है और उसी आकार का मैट्रिक्स देता है जहां आउटपुट सेल (i, j) पंक्तियों 1 और i के बीच इनपुट कॉलम जे के मानक विचलन पर सेट होता है। एनएएस को अनदेखा किया जाना चाहिए, जब तक कि इनपुट मैट्रिक्स का सेल (i, j) स्वयं NA न हो, आउटपुट मैट्रिक्स के मामले में सेल (i, j) भी NA होना चाहिए।

मुझे एक अंतर्निहित फ़ंक्शन नहीं मिला, इसलिए मैंने निम्नलिखित कोड लागू किया। दुर्भाग्यवश, यह एक लूप का उपयोग करता है जो बड़े matrices के लिए कुछ हद तक धीमा हो जाता है। क्या कोई तेज़ अंतर्निहित कार्य है या कोई बेहतर दृष्टिकोण सुझा सकता है?

cumsd <- function(mat) 
{ 
    retval <- mat*NA 
    for (i in 2:nrow(mat)) retval[i,] <- sd(mat[1:i,], na.rm=T) 
    retval[is.na(mat)] <- NA 
    retval 
} 

धन्यवाद।

Answer 1

आप cumsum इस्तेमाल कर सकते हैं मैट्रिक्स पर vectorized संचालन के विचरण/एसडी के लिए प्रत्यक्ष सूत्रों से आवश्यक रकम की गणना करने के:

cumsd_mod <- function(mat) { 
    cum_var <- function(x) { 
     ind_na <- !is.na(x) 
     nn <- cumsum(ind_na) 
     x[!ind_na] <- 0 
     cumsum(x^2)/(nn-1) - (cumsum(x))^2/(nn-1)/nn 
    } 
    v <- sqrt(apply(mat,2,cum_var)) 
    v[is.na(mat) | is.infinite(v)] <- NA 
    v 
} 

सिर्फ तुलना के लिए:

set.seed(2765374) 
X <- matrix(rnorm(1000),100,10) 
X[cbind(1:10,1:10)] <- NA # to have some NA's 

all.equal(cumsd(X),cumsd_mod(X)) 
# [1] TRUE 

और समय के बारे में:

X <- matrix(rnorm(100000),1000,100) 
system.time(cumsd(X)) 
# user system elapsed 
# 7.94 0.00 7.97 
system.time(cumsd_mod(X)) 
# user system elapsed 
# 0.03 0.00 0.03 

Answer 2

एक और प्रयास (मरेक तेज है)

cumsd2 <- function(y) { 
n <- nrow(y) 
apply(y,2,function(i) { 
    Xmeans <- lapply(1:n,function(z) rep(sum(i[1:z])/z,z)) 
    Xs <- sapply(1:n, function(z) i[1:z]) 
    sapply(2:n,function(z) sqrt(sum((Xs[[z]]-Xmeans[[z]])^2,na.rm = T)/(z-1))) 
}) 
}