r - R中的for循环用于创建向量的两组

Question

使用以下数据集

u_data
       rSLn rwave       rexpd     y_ij rwave2     u_ij
    1     1     1  199.929886 5.302956      1 5.302956
    2     1     2   27.738826 3.358249      4 3.358249
    3     1     3  144.000000 4.976734      9 4.976734
    4     1     4   72.000000 4.290459     16 4.290459
    5     1     5    0.000000 0.000000     25 0.000000
    6     2     1  392.606361 5.975351      1 5.975351
    7     2     2  749.524990 6.620773      4 6.620773
    8     2     3 3120.000000 8.045909      9 8.045909
    9     2     4 1600.000000 7.378384     16 7.378384
    10    2     5 1000.000000 6.908755     25 6.908755
    11    2     6 5840.000000 8.672657     36 8.672657
    12    2     7 3960.000000 8.284252     49 8.284252
    13    2     8 4700.000000 8.455531     64 8.455531
    14    2     9 1660.000000 7.415175     81 7.415175
    15    2    10 5620.000000 8.634265    100 8.634265
    16    3     1 1566.117441 7.356993      1 7.356993
    17    3     2  739.702016 6.607598      4 6.607598
    18    3     3    0.000000 0.000000      9 0.000000
    19    3     4    0.000000 0.000000     16 0.000000
    20    3     5    0.000000 0.000000     25 0.000000
    21    3     6    0.000000 0.000000     36 0.000000
    22    3     7    0.000000 0.000000     49 0.000000
    23    3     8    0.000000 0.000000     64 0.000000
    24    3     9  600.000000 6.398595     81 6.398595
    25    3    10  720.000000 6.580639    100 6.580639
    26    4     1  249.912358 5.525104      1 5.525104
    27    4     2    9.246275 2.326914      4 2.326914
    28    4     3  848.000000 6.744059      9 6.744059
    29    4     4  820.000000 6.710523     16 6.710523
    30    4     5  968.000000 6.876265     25 6.876265
    31    4     6 4800.000000 8.476580     36 8.476580
    32    4     7 1572.000000 7.360740     49 7.360740
    33    4     8 1960.000000 7.581210     64 7.581210
    34    4     9 1800.000000 7.496097     81 7.496097
    35    4    10 1700.000000 7.438972    100 7.438972
    36    5     1    0.000000 0.000000      1 0.000000
    37    5     2 6768.273444 8.820149      4 8.820149
    38    5     3  520.000000 6.255750      9 6.255750
    39    5     4 1020.000000 6.928538     16 6.928538
    40    5     5 1520.000000 7.327123     25 7.327123
    41    5     6 2075.000000 7.638198     36 7.638198
    42    5     7 1760.000000 7.473637     49 7.473637
    43    5     8 1270.000000 7.147559     64 7.147559
    44    5     9 5400.000000 8.594339     81 8.594339
    45    5    10 6550.000000 8.787373    100 8.787373

并具有以下值

ux_data=as.matrix(u_data[,c(2,5)])
ux_data=cbind(1, ux_data)
class=rbinom(length(unique(u_data$rSLn)),1,0.48)+1
thet.value=c(4.25,5.85,1.26,9.78,6.86)
n_g_i=numeric()
for ( d in unique(u_data$rSLn)){
  n_g_i[d]=length(u_data$rwave[u_data$rSLn==d])
}
sigma2=0.7849
SIGMA=matrix(c(100,0,0,
               0,1,0,
               0,0,1/100), nrow = 3, ncol = 3, byrow = T)

我想执行以下代码，这些代码运行良好。

 u_ij_C1=(u_data$u_ij[rep(class,times=n_g_i)==1] #u_ij_new belongs to cluster-1
             -rep(thet.value[class==1], n_g_i[class==1]))


    m_beta_C1=(solve((t(ux_data[rep(class,times=n_g_i)==1,])%*%ux_data[rep(class,times=n_g_i)==1,]/
                         (sigma2))+solve(SIGMA)) %*%(t(ux_data[rep(class,times=n_g_i)==1,])%*%u_ij_c1/sigma2))

    sig2_beta_C1=(solve((t(ux_data[rep(class,times=n_g_i)==1,])
                        %*%ux_data[rep(class,times=n_g_i)==1,]/(sigma2))+solve(SIGMA)))

    u_ij_C2=(u_data$u_ij[rep(class,times=n_g_i)==2] #u_ij_new belongs to cluster-2
             -rep(thet.value[class==2], n_g_i[class==2]))
    m_beta_C2=(solve((t(ux_data[rep(class,times=n_g_i)==2,])%*%ux_data[rep(class,times=n_g_i)==2,]/
                        (sigma2))+solve(SIGMA)) %*%(t(ux_data[rep(class,times=n_g_i)==2,])%*%u_ij_C2/sigma2))

    sig2_beta_C2=(solve((t(ux_data[rep(class,times=n_g_i)==2,])
                         %*%ux_data[rep(class,times=n_g_i)==2,]/(sigma2))+solve(SIGMA)))

每个m_beta都是大小为 3 的向量，并且是 我尝试使用 for 循环执行sig2_beta的顺序矩阵，不幸的是，它不起作用3x3

ngrp=2
u_ij_New_C12=numeric()
mu_beta_C12=numeric()
sig_beta_C12=array()
for ( k in 1:ngrp){
  u_ij_New_C12[k]=(u_data$u_ij[rep(class,times=n_g_i)==k] #u_ij_new belongs to cluster-k
                -rep(theta_i[class==k], n_g_i[class==k])) #repeting thetas belongs to cluster-k
  sig_beta_C12[k]=(solve((t(ux_data[rep(class,times=n_g_i)==k,])
                          %*%ux_data[rep(class,times=n_g_i)==k,]/
                            (sigma2))+solve(SIGMA)))

  mu_beta_C12[k]=(sig_beta_C12[k] %*%(t(ux_data[rep(class,times=n_g_i)==k,])%*%u_ij_New_C12[k]/sigma2))


} 

对于 k=1，我期望 和 相同的cluster-1结果cluster-2。例如mu_beta_C12[1]andsig_beta_C12[1]应该与m_beta_C1and完全相似sig2_beta_C1。任何帮助表示赞赏。

标签： rfor-loopmatrixvectorgroup-by