With the following data set
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
And with following values
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)
I would like to execute the following code, which are working perfectly.
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)))
Each m_beta is a vector of size 3 and sig2_beta is a matrix of order 3x3
I am trying to do it with for loop, Unfortunately, it is not working
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))
}
For k=1, I am expecting the same result for the cluster-1 and the same for cluster-2. For example mu_beta_C12[1] and sig_beta_C12[1] should be exactly similar to m_beta_C1 and sig2_beta_C1 respectively.
Any help is appreciated.
