-
Notifications
You must be signed in to change notification settings - Fork 1
/
Analysis.R
123 lines (98 loc) · 3.07 KB
/
Analysis.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
library(sandwich)
data(PublicSchools); ps <- na.omit(PublicSchools)
ps$Income <- ps$Income * 0.01
#Robust Regression
library(MASS); library(invgamma); library(MCMCpack)
K <- 5000
#Store distribution
LAMBDA <- BETA <- NULL
S2 <- rep(NA, K)
#Prior Values
v <- 4 #df
#Jeffrey's prior
a <- c(1, 3/2)
#Quadratic
#Data
ps[, 'Income2'] <- ps[,'Income']^2; ps[, 'Intercept'] <- 1
ps <- ps[c("Expenditure", "Intercept", "Income", "Income2")]
X <- as.matrix.data.frame(ps[-c(1)])
Y <- as.matrix.data.frame(ps[c(1)])
N <- dim(ps)[1]
# Starting value of BETA and S2
m1 <- lm(Expenditure ~ Income + I(Income^2), data = ps)
beta <- as.matrix(m1$coefficients)
sigma2 <- summary(m1)$sigma
#Current value
lambda <- rep(NA, N)
RSS <- rep(NA, N)
#Robust Regression
set.seed(143)
for (i in 1:K) {
#lambdas
for (n in 1:N){
rss.i <- (Y[n] - t(beta) %*% X[n, ])^2
RSS[n] <- rss.i
lambda[n] <- rinvgamma(1, (v +1)/2, (v + ((sigma2)^(-1))*rss.i)/2)
}
lambda.matrix <- diag(lambda, nrow = N, ncol = N)
LAMBDA <- rbind(lambda, LAMBDA)
#sigma2
wRSS <- sum(RSS/lambda)
sigma2 <- rinvgamma(1, (N + 2*(a[1] - 1))/2, (1/2)*wRSS)
S2[i] <- sigma2
#Beta
beta.sigma <- solve(t(X)%*% solve(sigma2*lambda.matrix) %*% X )
beta.mu <- beta.sigma %*% (t(X) %*% solve(sigma2*lambda.matrix) %*%Y )
beta <- as.matrix(mvrnorm(1, beta.mu, beta.sigma))
BETA <- rbind(t(beta),BETA)
}
#Standard Regression
#Quadratic
ols.beta <- as.array(m1$coefficients)
library(mvtnorm)
df <- N + a[1] - 3 - 1
s2 <- sum((m1$residuals)^2)
BETA.st <- rmvt(K, sigma = s2*solve((t(X) %*% X)), df =df, delta = ols.beta )
S2.st <- rinvgamma(K, df/2, s2/2)
#With no quardratic term
X <- as.matrix.data.frame(ps[-c(1, 4)])
Y <- as.matrix.data.frame(ps[c(1)])
N <- dim(ps)[1]
# Starting value of BETA and S2
m2 <- lm(Expenditure ~ Income, data = ps[-c(4)])
beta <- as.matrix(m2$coefficients)
sigma2 <- summary(m2)$sigma
#Store distribution
l.LAMBDA <- l.BETA <- NULL
l.S2 <- rep(NA, K)
set.seed(143)
for (i in 1:K) {
#lambdas
for (n in 1:N){
rss.i <- (Y[n] - t(beta) %*% X[n, ])^2
RSS[n] <- rss.i
lambda[n] <- rinvgamma(1, (v +1)/2, (v + ((sigma2)^(-1))*rss.i)/2)
}
lambda.matrix <- diag(lambda, nrow = N, ncol = N)
l.LAMBDA <- rbind(lambda, l.LAMBDA)
#sigma2
wRSS <- sum(RSS/lambda)
sigma2 <- rinvgamma(1, (N + 2*(a[1] - 1))/2, (1/2)*wRSS)
l.S2[i] <- sigma2
#Beta
beta.sigma <- solve(t(X)%*% solve(sigma2*lambda.matrix) %*% X )
beta.mu <- beta.sigma %*% (t(X) %*% solve(sigma2*lambda.matrix) %*%Y )
beta <- as.matrix(mvrnorm(1, beta.mu, beta.sigma))
l.BETA <- rbind(t(beta),l.BETA)
}
library("coda")
effectiveSize(BETA)
# Standard Regression
#Linear
ols.l.beta <- as.array(m2$coefficients)
df <- N + a[1] - 2 - 1
s2 <- sum((m2$residuals)^2)
l.BETA.st <- rmvt(K, sigma = s2*solve((t(X) %*% X)), df =df, delta = ols.l.beta )
l.S2.st <- rinvgamma(K, df/2, s2/2)
save.image(file="/Users/MacUser/Desktop/MA578/Final Project/analysis.RData")
save(BETA, LAMBDA, S2, BETA.st, l.BETA, l.BETA.st, S2.st, l.S2, l.LAMBDA, l.S2, l.S2.st, file="/Users/MacUser/Desktop/MA578/Final Project/AnalysisMain.RData")