-
Notifications
You must be signed in to change notification settings - Fork 0
/
GC_Sim1.R
128 lines (92 loc) · 2.88 KB
/
GC_Sim1.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
### SIMULATION ONE
#Sim 1
source("/uufs/chpc.utah.edu/common/home/u6011224/GenCorr/GenCorrFunct.R")
#source("GenCorr Functions.R")
#print("pig")
library(data.table)
library(R.utils)
library(MASS)
arg = cmdArgs()
#arg = list("FIRST", 890, 2)
# n = 60
# p = 3000
# q = 6
n = 120
p = 1500
q = 4
score = matrix(-1, nrow = 100, ncol = 3) ###These are the scores for each method.
#### Best rank, mid rank, worst rank. Perfect score is 1,2,3.
colnames(score) = c("GCorr_Best", "GCorr_Mid", "GCorr_Worst")
set.seed(arg[[2]])
for(r in 1:100)
{
X = matrix(NA, nrow = n, ncol = p)
for (i in 1:p)
{
#X[, i] = rnorm(n, 0, 5)
X[, i] = rnorm(n, 3, 1)
#X[, i] = rpois(n, 2)
}
Y = matrix(0, nrow = n, ncol = q)
# rho = 0.5
# Sigma = diag(q)
# Sigma = rho ^ abs(row(Sigma) - col(Sigma))
#
# scalarMat = t(mvrnorm(3, rep(0, q), Sigma))
# rho = 0.5
# Sigma = diag(q)
# Sigma = rho ^ abs(row(Sigma) - col(Sigma))
#
# scalarMat = t(mvrnorm(3, rep(3, q), Sigma))
sign = (-1)^(rbinom(n=3*q,size = 1, prob = 0.4))
scalarVals = sign*(4*log(n)/sqrt(n) + abs(rnorm(n = 3*q,0, 1)))
scalarMat = matrix(scalarVals, nrow = 3, ncol = q)
####This comes from DC-SIS. (I've adapted it to a MV setting)
for(m in 1:q)
{
Y[,m] = X[,c(1,2,3)] %*% scalarMat[,m]
}#
testData = cbind(Y, X)
#testData = matrix(rnorm(300, 0,1), ncol = 30, nrow = 10)
#phi = CovarProd_MV_main(data = testData, q = q)
phi_num2 = CorrMat_MV_main(data = testData, q = q, normToUse = arg[[3]])
omega = rep(2,3000)
#timestamp()
#omega = zllz(data = testData, q=q)
#timestamp()
dc = rep(1,3000)
#timestamp()
#dc = DCSIS_MV(data = testData, q=q)
#timestamp()
#
# head(phi)
# sorted_phi = sort.int(phi, decreasing = TRUE, index.return = TRUE)
# head(sorted_phi$ix, 10)
# phiMainScores = which(sorted_phi$ix %in% c(1,3,5))
#head(phi_num2)
sorted_phi_num2 = sort.int(phi_num2, decreasing = TRUE, index.return = TRUE)
#print(head(sorted_phi_num2$ix, 10))
phi_num2MainScores = which(sorted_phi_num2$ix %in% c(1,2,3))
# head(omega)
#sorted_omeg = sort.int(omega, decreasing = TRUE, index.return = TRUE)
# head(sorted_omeg$ix, 30)
#omegaMainScores = which(sorted_omeg$ix %in% c(1,2,3))
#head(dc)
#sorted_dc = sort.int(dc, decreasing = TRUE, index.return = TRUE)
#head(sorted_dc$ix, 30)
#dcMainScores = which(sorted_dc$ix %in% c(1,2,3))
score[r, 1:3] = phi_num2MainScores
#score[r, 4:6] = omegaMainScores
#score[r, 7:9] = dcMainScores
if(r %% 25 == 0)
{
print(paste0("We are at replicate ", r, " of Sim 1 run ", arg[[1]]))
print("The column scores are ")
print(colMeans(score[1:r,]))
}
}
fileName = paste0("/uufs/chpc.utah.edu/common/home/u6011224/GenCorr/", "Sim1Results_", arg[[1]],".csv")
fwrite(as.data.frame(score),
file = fileName,
showProgress = TRUE,
col.names=TRUE)