forked from davezes/Stats_141SL
-
Notifications
You must be signed in to change notification settings - Fork 0
/
_03_fixed_X_01.R
executable file
·74 lines (37 loc) · 1.04 KB
/
_03_fixed_X_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
options(stringsAsFactors=FALSE)
xbool_save_file <- FALSE
##################### barchart for proportions
nn <- 2 * 10^(5)
n <- 30
set.seed(780)
x <- rcauchy(n)
hist(x)
xsig <- 3
y <- x + rnorm(n, 0, xsig)
if(xbool_save_file) {
png(file.path("~", "Desktop", "highLeverageSim_01.png"), width=1000, height=1000, pointsize=24)
}
plot(x,y, cex=1.4, lwd=2)
if(xbool_save_file) { dev.off() }
xlm <- lm(y~x-1)
summary(xlm)
xsig_hat <- sqrt( sum( (y - xlm$fitted)^2 ) / (n-1) ) ; xsig_hat
xSEest <- xsig_hat / sqrt(sum(x^2)) ; xSEest
xSEtrue <- xsig / sqrt(sum(x^2)) ; xSEtrue
xSEs_vec <- numeric(nn)
kk <- 1
while(kk <= nn) {
y <- x + rnorm(n, 0, xsig)
xlm <- lm(y~x-1)
xSEs_vec[kk] <- summary(xlm)$coef[1,2]
kk <- kk + 1
if(kk %% 10000 == 0) { cat(kk, "\n") }
}
if(xbool_save_file) {
png(file.path("~", "Desktop", "distributionSimBhat_01.png"), width=1000, height=1000, pointsize=24)
}
hist(xSEs_vec, lwd=1.4)
abline(v=xSEtrue, lwd=3)
if(xbool_save_file) { dev.off() }
mean(xSEs_vec)
xSEtrue