Анализ нескольких повторений

Question

У меня есть набор данных с 1000 повторений.Каждое повторение содержит 50 записей.Мне нужно запустить и вывести некоторую статистику для каждого повторения.Кажется, цикл for - это правильный путь, но я не могу заставить его работать.

Пример набора данных (усеченный до 3 выборок на повторение) и код, который я хотел бы запустить на каждом наборе повторений

'''R
####DATASET
Rep Y_Star_T m
1 3 3 
1 8 7
1 9 6
2 13 2 
2 5 5
2 19 16
3 12 7 
3  7
3 9 6
####Global variables needed outside the loop
a <- .25
L <- 75
W <- 20
big_N <- (L*W)/a

####Begin appended calculations needed for the loop
w <- netdata$y_star_T/netdata$m
mu_hat <- (1/n *sum(w))/a
tau_hat <- (big_N*a)*mu_hat
var_mu_hat <- (1/(n*(n-1))*sum((w-mu_hat)^2))/a^2
var_tau_hat <- (big_N*a)^2*var_mu_hat

Я бы хотел, чтобы вектор длины равнялся числуповторы, которые содержат значения, рассчитанные выше (например, mu_hat).Для 1000 повторений длина вектора была бы 1000 Один вектор на mu_hat, tau_hat и т. Д. - это хорошо, если это просто, или список списков также будет работать.

Answer 1

Вы должны определить вспомогательную функцию для вычисления необходимой статистики.Тогда *apply, которые функционируют по группам, определенным столбцом Rep.

funStats <- function(DF, a = 0.25, L = 75, W = 20, big_N = L*W/a){
    n <- nrow(DF)
    w <- DF[["Y_Star_T"]]/DF["m"]
    mu_hat <- (1/n *sum(w))/a
    tau_hat <- (big_N*a)*mu_hat
    var_mu_hat <- (1/(n*(n-1))*sum((w-mu_hat)^2))/a^2
    var_tau_hat <- (big_N*a)^2*var_mu_hat
    res <- c(mu_hat = mu_hat,
      tau_hat = tau_hat,
      var_mu_hat = var_mu_hat,
      var_tau_hat = var_tau_hat
    )
    res
}

sapply(split(netdata, netdata$Rep), funStats)
#                       1            2            3
#mu_hat      4.857143e+00 1.158333e+01 6.428571e+00
#tau_hat     7.285714e+03 1.737500e+04 9.642857e+03
#var_mu_hat  1.065170e+02 6.557882e+02 3.721224e+02
#var_tau_hat 2.396633e+08 1.475523e+09 8.372755e+08


result <- lapply(split(netdata, netdata$Rep), funStats)
result <- do.call(rbind, result)

result
#     mu_hat   tau_hat var_mu_hat var_tau_hat
#1  4.857143  7285.714   106.5170   239663265
#2 11.583333 17375.000   655.7882  1475523437
#3  6.428571  9642.857   372.1224   837275510

Данные.

Обратите внимание, что строка 8 закомментирована.

netdata <- read.table(text = "
Rep Y_Star_T m
1 3 3 
1 8 7
1 9 6
2 13 2 
2 5 5
2 19 16
3 12 7 
#3  7
3 9 6
", header = TRUE)