-
Notifications
You must be signed in to change notification settings - Fork 0
/
GC_Sim7.R
128 lines (92 loc) · 2.95 KB
/
GC_Sim7.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 Seven
#Sim 7
#
#setwd("M:/GCorr")
source("/uufs/chpc.utah.edu/common/home/u6011224/GenCorr/GenCorrFunct.R")
#source("M:/GCorr/GenCorr Functions.R")
#print("pig")
library(data.table)
library(R.utils)
library(MASS)
#arg = list("FIRST", 132,2)
arg = cmdArgs()
n = 100
p = 1000
q = 4
d = 5
score = matrix(-1, nrow = 100, ncol = 4) ###These are the scores for each method.
#### Best rank, worst rank. Perfect score is 1,2.
colnames(score) = c("GCorr_Best", "GCorr_Worst", "X_1X_2", "X_3X_4")
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, 2)
#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(2, rep(3, q), Sigma))
# ##scalarMat = matrix(rnorm(18, 0,2.5), nrow = q, ncol = 3)
#
# for(m in 1:q)
# {
#
# Y[,m] = cbind(X[,1]*X[,2], X[,3]*X[,4]) %*% scalarMat[m,]
# }###Should create a nice little linear model like above (around lines 28-34).
# #
Y[,1] = X[, 1]*X[,2] + 2*X[,4]*X[,3] + rnorm(n,0,1)### plus some small error!
Y[,2] = -X[, 1]*X[,2] + 1.5*X[,4]*X[,3] + rnorm(n,0,1)
Y[,3] = -2.5*X[, 1]*X[,2] - 2*X[,4]*X[,3] + rnorm(n,0,1)
Y[,4] = 2*X[, 1]*X[,2] +X[,4]*X[,3] + rnorm(n,0,1)
# # Y[,5] = -2*X[, 1] *-1*X[,3]# +X[,5]
# # Y[,6] = X[, 1]* -1.5*X[,3]# -2.5*X[,5]
# #
testData = cbind(Y, X)
#timestamp()
phi_num2 = CorrMat_MV_Interact(data = testData, q = q, normToUse = arg[[3]])
timestamp()
#head(phi_num2[1:20, 1:20])
sorted_phi_num2 = sort.int(phi_num2, decreasing = TRUE, index.return = TRUE)
###Using a temp cutoff of d=5
# selectedPairs = head(sorted_phi_num2$ix, d)
# for(k in 1:d)
# {
# j1 = (selectedPairs[k]-1) %% p + 1
# j2 = (selectedPairs[k]-1) %/% p + 1
#
# print(paste0("(", j1, ",", j2, ")"))
# print(selectedPairs[k])
#
# }
phi_num2MainScores = which(sorted_phi_num2$ix %in% c(1001,3003))
X1X2Score = which(sorted_phi_num2$ix == 1001)
X3X4Score = which(sorted_phi_num2$ix == 3003)
score[r, ] = c(phi_num2MainScores, X1X2Score, X3X4Score)
# # 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,4,5))
#
# #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,4,5))
if(r %% 25 == 0)
{
##timestamp()
print(paste0("We are at replicate ", r, " for Sim 7 run ", arg[[1]]))
print("The column scores are ")
print(colMeans(score[1:r,]))
}
}
fileName = paste0("/uufs/chpc.utah.edu/common/home/u6011224/GenCorr/", "Sim7Results_", arg[[1]],".csv")
fwrite(as.data.frame(score),
file = fileName,
showProgress = TRUE,
col.names=TRUE)