forked from davezes/Stats_141SL
-
Notifications
You must be signed in to change notification settings - Fork 0
/
_04_fixed_XvsRndX_wishart_01.R
161 lines (91 loc) · 3.13 KB
/
_04_fixed_XvsRndX_wishart_01.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
options(stringsAsFactors=FALSE, width=350)
library(mvtnorm)
N <- 20
xbool_save_file <- TRUE
set.seed(777) ### try different seeds to see how different slope pred variance is between fixed and random inputs
b <- c(1, -1)
xSig <- diag(0.5, 2) + 0.5
###### look at Cov mx
xSig
xmu <- c(0, 0)
######### X is fixed -- not subject to sampling variation
X <- rmvnorm(N, xmu, xSig)
f.true <- X %*% b
sig <- 20
nn <- 300000
mxbhats <- matrix(NA, nn, length(b))
for(iinn in 1:nn) {
y <- f.true + rnorm(N, 0, sig)
b.hat <- solve(crossprod(X)) %*% (t(X) %*% y)
mxbhats[ iinn, ] <- b.hat
}
######## true SEs
sqrt( diag(solve( N * xSig )) ) * sig
xxx <- apply( t(mxbhats) - b, 1, function(x) { return( ( mean( x^2 ) ) ) } )
##### agrees with actual sim error
##### sim var bhat
cat(
"Under fixed X, this is our simulated pred error of bhat: ",
xxx,
"\n"
)
#### hist(mxbhats[ , 1])
####### var using fixed X
cat(
"Under fixed X, this is our theoretical (using true resid variance) pred error of bhat: ",
diag(solve( crossprod(X) )) * sig^2,
"\n"
)
## apply(mxbhats, 2, sd)
if(xbool_save_file) {
png(file.path("~", "Desktop", "distOfBhats_Xfixed_01.png"), width=1800, height=1000, pointsize=24)
}
par(mfrow=c(1, 2))
hist(mxbhats[ ,1], cex=1.4, lwd=2, xlim=c(-40, 40), main=paste0("Sim b1 with fixed X -- sim pred var = ", round(xxx[1], 2)))
abline(v=b[1], col="#00FF00", lwd=4)
hist(mxbhats[ ,2], cex=1.4, lwd=2, xlim=c(-40, 40), main=paste0("Sim b2 with fixed X -- sim pred var = ", round(xxx[2], 2)))
abline(v=b[2], col="#00FF00", lwd=4)
if(xbool_save_file) { dev.off() }
cat("\n")
##############################
N <- 20
sig <- 20
nn <- 300000
mxbhats <- matrix(NA, nn, length(b))
arinvXX <- array(NA, c(length(b), length(b), nn))
for(iinn in 1:nn) {
X <- rmvnorm(N, xmu, xSig) #### X is NOT fixed -- but rather subject to sampling variation
f.true <- X %*% b
y <- f.true + rnorm(N, 0, sig)
invXX <- solve(crossprod(X))
arinvXX[ , , iinn] <- invXX
b.hat <- invXX %*% (t(X) %*% y)
mxbhats[ iinn, ] <- b.hat
}
xxx2 <- apply( t(mxbhats) - b, 1, function(x) { return( ( mean( x^2 ) ) ) } )
apply(arinvXX, c(1,2), mean) * sig^2
##### var bhat
##### sim var bhat
cat(
"Under random X, this is our simulated pred error of bhat: ",
xxx2,
"\n"
)
######### the center of the inverse Wishart distribution:
xmeanInvWishart <- solve( xSig ) / (N - length(b) - 1)
cat(
"Under random X, this is our theoretical (using true resid variance) pred error of bhat: ",
diag(sig^2 * xmeanInvWishart),
"\n",
"note that we use the expectation of the inverse Wishart for this.",
"\n"
)
if(xbool_save_file) {
png(file.path("~", "Desktop", "distOfBhats_Xrnd_01.png"), width=1800, height=1000, pointsize=24)
}
par(mfrow=c(1, 2))
hist(mxbhats[ ,1], cex=1.4, lwd=2, xlim=c(-40, 40), main=paste0("Sim b1 with random X -- sim pred var = ", round(xxx2[1], 2)))
abline(v=b[1], col="#00FF00", lwd=4)
hist(mxbhats[ ,2], cex=1.4, lwd=2, xlim=c(-40, 40), main=paste0("Sim b2 with random X -- sim pred var = ", round(xxx2[2], 2)))
abline(v=b[2], col="#00FF00", lwd=4)
if(xbool_save_file) { dev.off() }