首页 > 解决方案 > 从每一行产生一个矩阵(tidyverse)

问题描述

我正在尝试跨后验矩阵进行分析。我开始的是k ^2 列的小标题,其中k是矩阵的维度。第i行形成第i次迭代的矩阵。

因此,例如对于 3x3 矩阵,这是:

set.seed(12)
n <- 1000
z1z1 <- rnorm(n, 5, 1)
z2z2 <- rnorm(n, 5, 1)
z3z3 <- rnorm(n, 5, 1)
z1z2 <- rnorm(n, 0, 1)
z1z3 <- rnorm(n, 0, 1)
z2z3 <- rnorm(n, 0, 1)
post3 <- as_tibble(matrix(c(z1z1, z1z2, z1z3,
                            z1z2, z2z2, z2z3,
                            z1z3, z2z3, z3z3), 
                          ncol = 9))
post3

给予:

# A tibble: 1,000 x 9
      V1     V2      V3     V4    V5      V6      V7      V8    V9
   <dbl>  <dbl>   <dbl>  <dbl> <dbl>   <dbl>   <dbl>   <dbl> <dbl>
 1  3.52 -0.618  2.96   -0.618  2.48 -0.634   2.96   -0.634   5.98
 2  6.58 -0.827  0.0909 -0.827  5.52 -1.84    0.0909 -1.84    6.20
 3  4.04  1.48  -1.66    1.48   6.58  0.166  -1.66    0.166   5.58
 4  4.08 -1.01   0.809  -1.01   5.49  0.607   0.809   0.607   6.55
 5  3.00  0.582 -0.485   0.582  6.20  0.0765 -0.485   0.0765  6.38
 6  4.73  0.718  1.97    0.718  4.00 -0.147   1.97   -0.147   4.35
 7  4.68 -0.372  0.572  -0.372  4.65 -1.68    0.572  -1.68    3.83
 8  4.37 -0.809  0.883  -0.809  3.96  0.985   0.883   0.985   4.97
 9  4.89  0.405  0.686   0.405  6.02  0.252   0.686   0.252   6.29
10  5.43  0.124  0.199   0.124  5.75  0.354   0.199   0.354   4.20
# ... with 990 more rows

这是第一次迭代中的矩阵:

k <- sqrt(length(post3))
matrix(post3[1,], nrow = k)
     [,1]      [,2]       [,3]      
[1,] 3.519432  -0.618137  2.962622  
[2,] -0.618137 2.479522   -0.6338298
[3,] 2.962622  -0.6338298 5.977552

然后我沿着这个后验计算第一个特征向量的优势:

post3 %>%
  rowwise %>%
  mutate(
    pre_eig = list(eigen(matrix(c(V1, V2, V3, V4, V5, V6, V7, V8, V9), nrow = k))),
    dom = pre_eig[[1]][1] / sum(pre_eig[[1]][1:k])) %>%
  select('dom')

给予:

# A tibble: 1,000 x 1
     dom
   <dbl>
 1 0.676
 2 0.437
 3 0.462
 4 0.427
 5 0.414
 6 0.504
 7 0.474
 8 0.429
 9 0.394
10 0.383
# ... with 990 more rows

我想做的是使这个脚本具有通用性,以便它可以对任何 k 值进行后验。我遇到的问题是如何定义矩阵而不必手写所有列名- 将其应用于 2000x2000 矩阵时,我不想写出V1, V2, V3... V4000000!我尝试了一些东西(包括...... eigen(matrix(c(paste0('V', 1:(k^2))), nrow = k))),我认为它不起作用,因为它想要V1, V2...而不是"V1", "V2"...)并且我都没有想法。如何让它自动从后小标题中获取列名?

然后我就可以使用完全相同的脚本,例如 post3 <- as_tibble(matrix(c(z1z1, z1z2, z1z2, z2z2), ncol = 4))...

标签: rdplyr

解决方案


如果将每一行的值收集到键值对中,则可以避免显式命名所有列:

library(tidyr)

post3 %>%
  # add row ID (so that results can be sorted back into original order)
  mutate(row.id = seq(1, n())) %>%

  # convert each row to long format, with values sorted from 1st to k^2th column
  gather(position, value, -row.id) %>%
  mutate(position = as.numeric(gsub("^V", "", position))) %>%
  arrange(row.id, position) %>% 
  select(-position) %>%

  # group by row ID & calculate
  group_by(row.id) %>%
  summarise(pre_eig = list(eigen(matrix(value, nrow = k))[["values"]]),
            dom = pre_eig[[1]][1] / sum(pre_eig[[1]][1:k])) %>%
  ungroup() %>%

  # sort results in original order
  arrange(row.id) %>%
  select(dom)

结果应该和以前一样:

# A tibble: 1,000 x 1
     dom
   <dbl>
 1 0.676
 2 0.437
 3 0.462
 4 0.427
 5 0.414
 6 0.504
 7 0.474
 8 0.429
 9 0.394
10 0.383
# ... with 990 more rows

推荐阅读