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
|
set.seed(12345)
x = 1:60
beta0_actual = 20
beta1_actual = 5
sigma_actual = 15
y = x*5 + 20 + rnorm(60, 0, 15)
plot(y ~ x)
linmod <- lm(y ~ x)
summary(linmod)
LogLikelihood <- function(beta0, beta1, sigma){
R = y - x*beta1 - beta0
sum(dnorm(R, mean = 0, sigma, log = TRUE))
}
prior <- function(beta0, beta1, sigma){
beta0prior = dunif(beta0, min=-100, max=100, log = TRUE)
beta1prior = dunif(beta1, min=-100, max=100, log = TRUE)
sigmaprior = dnorm(sigma, mean = 0, sd = 10, log = TRUE)
return(beta0prior + beta1prior + sigmaprior)
}
proposalfunction <- function(beta0, beta1, sigma){
return(rnorm(3, mean = c(beta0, beta1, sigma), sd= c(0.1, 0.5, 0.3)))
}
posterior <- function(beta0, beta1, sigma){
return (LogLikelihood(beta0, beta1, sigma) + prior(beta0, beta1, sigma))
}
MCMC <- function(intialvalue, generations){
chain = matrix(rep(NA, ((generations+1) * 3)), ncol = 3)
colnames(chain) <- c("beta0", "beta1", "sigma")
chain[1,] = intialvalue
for (i in 1:generations){
proposal = proposalfunction(chain[i,1], chain[i,2], chain[i,3])
probab = exp(posterior(proposal[1], proposal[2], proposal[3]) - posterior(chain[i,1], chain[i,2], chain[i,3]))
if (runif(1) < probab){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
}
return(chain)
}
res.chain1 <- MCMC(intialvalue = c(4,20,1), generations = 100000)
res.chain2 <- MCMC(intialvalue = c(4,20,1), generations = 100000)
par(mfrow = c(1, 3))
plot(res.chain1[,1], res.chain1[,2], type = "l")
plot(res.chain1[,1], res.chain1[,3], type = "l")
plot(res.chain1[,2], res.chain1[,3], type = "l")
图1 估计的参数的Burning In及马尔可夫链收敛后的变化
library(coda)
mcmc.res <- mcmc(res.chain1)
summary(mcmc.res)
plot(mcmc.res)
mcmcchains = mcmc.list(mcmc(res.chain1), mcmc(res.chain2))
plot(mcmcchains)
图2 trace plot以及概率密度分布
gelman.diag(mcmcchains)
gelman.plot(mcmcchains)
geweke.plot(mcmcchains)
geweke.diag(mcmcchains)
图3 用gelman.plot检验马尔可夫链是否收敛
autocorr.plot(mcmc(res.chain1))
densityplot(mcmcchains)
par(mfrow = c(1, 3))
densplot(mcmcchains)
par(mfrow = c(1, 3))
traceplot(mcmc(res.chain1))
xyplot(mcmc(res.chain1))
beta0_actual = 20
beta1_actual = 5
sigma_actual = 15
burnIn = 20000
acceptance = 1 - mean(duplicated(res.chain1[-(1:burnIn),]))
par(mfrow = c(2,3))
hist(res.chain1[-(1:burnIn),1],nclass=50, , main="Posterior of beta0", xlab="True value = red line" )
abline(v = mean(res.chain1[-(1:burnIn),1]), col = "blue")
abline(v = linmod$coefficients[1], col="red" )
hist(res.chain1[-(1:burnIn),2],nclass=50, main="Posterior of beta1", xlab="True value = red line")
abline(v = mean(res.chain1[-(1:burnIn),2]), col = "blue")
abline(v = linmod$coefficients[2], col="red" )
hist(res.chain1[-(1:burnIn),3],nclass=50, main="Posterior of sigma", xlab="True value = red line")
abline(v = mean(res.chain1[-(1:burnIn),3]) , col = "blue")
abline(v = sigma_actual, col="red" )
plot(res.chain1[-(1:burnIn),1], type = "l", xlab="Generations" , ylab = "Beta0", main = "Chain values of beta0", )
abline(h = linmod$coefficients[1], col="red" )
plot(res.chain1[-(1:burnIn),2], type = "l", xlab="Generations" , ylab = "Beta1", main = "Chain values of beta1", )
abline(h = linmod$coefficients[2], col="red" )
plot(res.chain1[-(1:burnIn),3], type = "l", xlab="Generations" , ylab = "Sigma", main = "Chain values of sd", )
abline(h = sigma_actual, col="red" )
|