0

给定纵向数据，我如何计算一个矩阵，其中每列代表给定变量的加权分组平均值？

我开发了一种需要循环的方法，而且速度太慢。我认为这可能是矢量化的，但解决方案让我望而却步。

这是我目前的方法：

library(foreach)

# N is sample size
# g is the number of groups
# p is the number of variables
get_group_mean_matrix <- function(N, g, p){
  X <- matrix(rbinom(N*p, 10, .5), N)
  f <- sort((1:(N)) %% g + 1)
  w <- runif(N)
  dmmat <- foreach(i = unique(f), .combine = rbind) %do% {
    idx <- which(f == i)
    ws <- w[idx]/sum(w[idx])
    t((t(X[idx,]) %*% ws)) %x% rep(1, length(idx))
  }
  dmmat
}

> set.seed(666)
> get_group_mean_matrix(12, 3, 5)
          [,1]     [,2]     [,3]     [,4]     [,5]
 [1,] 5.261103 4.074266 5.828070 4.452703 5.990165
 [2,] 5.261103 4.074266 5.828070 4.452703 5.990165
 [3,] 5.261103 4.074266 5.828070 4.452703 5.990165
 [4,] 5.261103 4.074266 5.828070 4.452703 5.990165
 [5,] 5.560556 4.241942 3.698828 5.572523 4.212532
 [6,] 5.560556 4.241942 3.698828 5.572523 4.212532
 [7,] 5.560556 4.241942 3.698828 5.572523 4.212532
 [8,] 5.560556 4.241942 3.698828 5.572523 4.212532
 [9,] 4.289029 4.771115 5.150607 4.424339 6.346775
[10,] 4.289029 4.771115 5.150607 4.424339 6.346775
[11,] 4.289029 4.771115 5.150607 4.424339 6.346775
[12,] 4.289029 4.771115 5.150607 4.424339 6.346775
> library(microbenchmark)
> microbenchmark(get_group_mean_matrix(1200, 300, 50))
Unit: milliseconds
                                 expr      min       lq     mean   median       uq      max neval
 get_group_mean_matrix(1200, 300, 50) 76.33337 77.39607 80.76586 78.39808 84.46984 93.40047   100

最初，我尝试使用 执行此操作lfe::demeanlist，但它给了我错误的输出！

library(lfe)
get_group_mean_matrix_lfe <- function(N, g, p){
  X <- matrix(rbinom(N*p, 10, .5), N)
  f <- sort((1:(N)) %% g + 1)
  w <- runif(N)
  X - demeanlist(X, list(factor(f)), weights = w)
}
> set.seed(666)
> get_group_mean_matrix_lfe(12, 3, 5)
          [,1]     [,2]     [,3]     [,4]     [,5]
 [1,] 5.138068 4.001781 5.415467 4.722947 5.999827
 [2,] 5.138068 4.001781 5.415467 4.722947 5.999827
 [3,] 5.138068 4.001781 5.415467 4.722947 5.999827
 [4,] 5.138068 4.001781 5.415467 4.722947 5.999827
 [5,] 5.197308 4.067657 3.202478 5.866451 4.066385
 [6,] 5.197308 4.067657 3.202478 5.866451 4.066385
 [7,] 5.197308 4.067657 3.202478 5.866451 4.066385
 [8,] 5.197308 4.067657 3.202478 5.866451 4.066385
 [9,] 4.189951 4.887720 4.953305 4.501874 6.385846
[10,] 4.189951 4.887720 4.953305 4.501874 6.385846
[11,] 4.189951 4.887720 4.953305 4.501874 6.385846
[12,] 4.189951 4.887720 4.953305 4.501874 6.385846
> library(microbenchmark)
> microbenchmark(get_group_mean_matrix_lfe(1200, 300, 50))
Unit: milliseconds
                                     expr      min       lq     mean   median       uq      max neval
 get_group_mean_matrix_lfe(1200, 300, 50) 6.107421 6.202426 6.500411 6.293648 6.582943 8.350876   100

虽然速度快很多...

我会接受以下两种答案中的任何一种：

  1. 说明lfe::demeanlist在加权情况下所做的事情。当我从平均值中减去加权偏差时，我不应该得到加权平均值吗？知道这一点，我如何计算加权分组均值矩阵？
  2. 不涉及 demeanlist 计算加权分组均值矩阵的方法。

注意：用%*%矩阵乘法函数代替可以RcppEigen加快速度，但还不够。我认为问题在于循环。

这是一些示例输入： 

   f X1 X2 X3 X4 X5
1  1  6  5  7  3  6
2  1  6  4  5  5  6
3  1  5  6  3  6  6
4  1  3  5  4  3  5
5  2  5  4  7  7  7
6  2  4  1  4  2  6
7  2  5  6  6  6  5
8  2  6  7  2  5  4
9  3  5  3  4  6  9
10 3  6  6  5  5  6
11 3  5  7  4  6  8
12 3  5  3  7  8  6

哪里f是分组因子。

4

1 回答 1

1

Hurr durr 我所要做的就是将权重的平方根变成demeanlisthurr durr

library(foreach)
get_group_mean_matrix <- function(N, g, p){
  X <- matrix(rbinom(N*p, 10, .5), N)
  f <- sort((1:(N)) %% g + 1)
  w <- runif(N)
  dmmat <- foreach(i = unique(f), .combine = rbind) %do% {
    idx <- which(f == i)
    ws <- w[idx]/sum(w[idx])
    t((t(X[idx,]) %*% ws)) %x% rep(1, length(idx))
  }
  dmmat
}

set.seed(666)
A <- get_group_mean_matrix(12, 3, 5)

library(lfe)
get_group_mean_matrix_lfe <- function(N, g, p){
  X <- matrix(rbinom(N*p, 10, .5), N)
  f <- sort((1:(N)) %% g + 1)
  w <- runif(N)
  X - demeanlist(X, list(factor(f)), weights = w^.5)
}

set.seed(666)
B <- get_group_mean_matrix_lfe(12, 3, 5)

> all.equal(A, B)
[1] TRUE
于 2018-08-29T20:03:28.900 回答