-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfigure_2.R
188 lines (171 loc) · 7.23 KB
/
figure_2.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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
library(scatterplot3d)
dire <- "figures/"
if(!dir.exists(dire)) dir.create(dire)
fileprefix <- paste0(dire, "figure_2")
set.seed(123)
n <- 3
r <- 100000
## HS vanilla
zeta <- abs(rt(r*3, 1))
theta <- matrix(abs(rt(r*3*n,1)), r*3,n)
phi_hs <- matrix(rnorm(r*3*n,0,theta*zeta),r*3,n)
png(filename = paste0(fileprefix, "a.png"),
width = 5, height = 5, units = "in", res = 300)
# jpeg(filename = paste0(fileprefix, "a.jpg"),
# width = 5, height = 5, units = "in", res = 300, quality = 100)
# pdf(paste0(fileprefix, "a.pdf"),
# width = 5, height = 5)
scatterplot3d::scatterplot3d(phi_hs, pch=".",
xlab = "",ylab = "",
zlab = "", tick.marks = FALSE,
color = scales::alpha(1, alpha = .5),
cex.main=.5,
box = FALSE,y.margin.add = 0 ,mar = c(.5,.5,-.1, 0) + 0.1,
xlim = c(-100,100),
ylim = c(-100,100),
zlim = c(-100,100))
mtext(expression(phi[1]), side = 1, adj = 0.3, line = -.5)
dims <- par("usr")
text(dims[1]+ 0.77*diff(dims[1:2]),dims[3]+ 0.15*diff(dims[3:4]),expression(phi[2]),srt=45)
text(dims[1]+ 0.015*diff(dims[1:2]), dims[3]+ 0.35*diff(dims[3:4]), expression(phi[3]))
dev.off()
## HS structured
zeta2 <- abs(rt(r*3, 1))
phi_hs2 <- phi_hs
phi_hs2[,3] <- rnorm(r*3,0,theta[,3]*zeta2)
png(filename = paste0(fileprefix, "d.png"),
width = 5, height = 5, units = "in", res = 300)
# jpeg(filename = paste0(fileprefix, "d.jpg"),
# width = 5, height = 5, units = "in", res = 300, quality = 100)
# pdf(paste0(fileprefix, "d.pdf"),
# width = 5, height = 5)
scatterplot3d::scatterplot3d(phi_hs2, pch=".",
xlab = "",ylab = "",
zlab = "", tick.marks = FALSE,
color = scales::alpha(1, alpha = .5),
cex.main=.5,
box = FALSE,y.margin.add = 0 ,mar = c(.5,.5,-.1, 0) + 0.1,
xlim = c(-100,100),
ylim = c(-100,100),
zlim = c(-100,100))
mtext(expression(phi[1]), side = 1, adj = 0.3, line = -.5)
dims <- par("usr")
text(dims[1]+ 0.77*diff(dims[1:2]),dims[3]+ 0.15*diff(dims[3:4]),expression(phi[2]),srt=45)
text(dims[1]+ 0.015*diff(dims[1:2]), dims[3]+ 0.35*diff(dims[3:4]), expression(phi[3]))
dev.off()
## SSVS vanilla
p <- rbeta(r, .1,.1)
gammas <- matrix(
rbinom(r*n, 1, p),
nrow = r
)
tau <- gammas
tau[which(gammas==0)] <- 0.01
tau[which(gammas==1)] <- 100
phi_ssvs <- matrix(
rnorm(r*n, 0, as.vector(tau)),
nrow = r
)
png(filename = paste0(fileprefix, "b.png"),
width = 5, height = 5, units = "in", res = 300)
# jpeg(filename = paste0(fileprefix, "b.jpg"),
# width = 5, height = 5, units = "in", res = 300, quality = 100)
# pdf(paste0(fileprefix, "b.pdf"),
# width = 5, height = 5)
scatterplot3d::scatterplot3d(phi_ssvs, pch=".",
xlab = "",ylab = "",
zlab = "", tick.marks = FALSE,
color = scales::alpha(1, alpha = 0.5),
cex.main=.5,
box = FALSE,y.margin.add = 0 ,
mar = c(.5,.5,-.1, 0) + 0.1,
xlim = c(-200,200),
ylim = c(-200,200),
zlim = c(-200,200))
mtext(expression(phi[1]), side = 1, adj = 0.3, line = -.5)
dims <- par("usr")
text(dims[1]+ 0.75*diff(dims[1:2]),dims[3]+ 0.15*diff(dims[3:4]),expression(phi[2]),srt=45)
text(dims[1]+ 0.015*diff(dims[1:2]), dims[3]+ 0.35*diff(dims[3:4]), expression(phi[3]))
dev.off()
## SSVS structured
p2 <- rbeta(r, .1,.1)
phi_ssvs2 <- phi_ssvs
gammas2 <- rbinom(r,1,p2)
tau2 <- gammas2*100
tau2[gammas2==0] <- .01
phi_ssvs2[,3] <- rnorm(r, 0, tau2)
png(filename = paste0(fileprefix, "e.png"),
width = 5, height = 5, units = "in", res = 300)
# jpeg(filename = paste0(fileprefix, "e.jpg"),
# width = 5, height = 5, units = "in", res = 300, quality = 100)
# pdf(paste0(fileprefix, "e.pdf"),
# width = 5, height = 5)
scatterplot3d::scatterplot3d(phi_ssvs2, pch=".",
xlab = "",ylab = "",
zlab = "", tick.marks = FALSE,
color = scales::alpha(1, alpha = 0.5),
cex.main=.5,
box = FALSE,y.margin.add = 0 ,
mar = c(.5,.5,-.1, 0) + 0.1,
xlim = c(-200,200),
ylim = c(-200,200),
zlim = c(-200,200))
mtext(expression(phi[1]), side = 1, adj = 0.3, line = -.5)
dims <- par("usr")
text(dims[1]+ 0.75*diff(dims[1:2]),dims[3]+ 0.15*diff(dims[3:4]),expression(phi[2]),srt=45)
text(dims[1]+ 0.015*diff(dims[1:2]), dims[3]+ 0.35*diff(dims[3:4]), expression(phi[3]))
dev.off()
## SHM within
lambda1 <- rgamma(r, .01,.01)
lambda2 <- rgamma(r, .01,.01)
phi1 <- matrix(
rnorm(n*r,0,sqrt(lambda1)),
r,n
)
png(filename = paste0(fileprefix, "c.png"),
width = 5, height = 5, units = "in", res = 300)
# jpeg(filename = paste0(fileprefix, "c.jpg"),
# width = 5, height = 5, units = "in", res = 300, quality = 100)
# pdf(paste0(fileprefix, "c.pdf"),
# width = 5, height = 5)
scatterplot3d::scatterplot3d(phi1, pch=".",
xlab = "",ylab = "",
zlab = "", tick.marks = FALSE,
color = scales::alpha(1, alpha = 0.5),
cex.main=.5,
box = FALSE,y.margin.add = 0 ,
mar = c(.5,.5,-.1, 0) + 0.1,
xlim = c(-2,2),
ylim = c(-2,2),
zlim = c(-2,2))
mtext(expression(phi[1]), side = 1, adj = 0.3, line = -.5)
dims <- par("usr")
text(dims[1]+ 0.77*diff(dims[1:2]),dims[3]+ 0.15*diff(dims[3:4]),expression(phi[2]),srt=45)
text(dims[1]+ 0.015*diff(dims[1:2]), dims[3]+ 0.35*diff(dims[3:4]), expression(phi[3]))
dev.off()
## SHM cross- vs ownlag
phi2 <- matrix(
rnorm(n*r,0,sqrt(lambda2)),
r,n
)
png(filename = paste0(fileprefix, "f.png"),
width = 5, height = 5, units = "in", res = 300)
# jpeg(filename = paste0(fileprefix, "f.jpg"),
# width = 5, height = 5, units = "in", res = 300, quality = 100)
# pdf(paste0(fileprefix, "f.pdf"),
# width = 5, height = 5)
scatterplot3d::scatterplot3d(phi1[,1],phi1[,2],phi2[,1], pch=".",
xlab = "",ylab = "",
zlab = "", tick.marks = FALSE,
color = scales::alpha(1, alpha = 0.5),#.2
cex.main=.5,
box = FALSE,y.margin.add = 0,
mar = c(.5,.5,-.1, 0) + 0.1,
xlim = c(-2,2),
ylim = c(-2,2),
zlim = c(-2,2))
mtext(expression(phi[1]), side = 1, adj = 0.3, line = -.5)
dims <- par("usr")
text(dims[1]+ 0.77*diff(dims[1:2]),dims[3]+ 0.15*diff(dims[3:4]),expression(phi[2]),srt=45)
text(dims[1]+ 0.015*diff(dims[1:2]), dims[3]+ 0.35*diff(dims[3:4]), expression(phi[3]))
dev.off()