-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathBayes_Normal.R
94 lines (60 loc) · 1.83 KB
/
Bayes_Normal.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
require(runjags)
set.seed(432104)
n <- 1000
x <- rnorm(n, 0, 5)
## Usando o modelo já implemetado:
model.Gauss <-"
model {
for (i in 1:n){
x[i] ~ dnorm(mu, phi)
}
mu ~ dnorm(0,.0001)
phi <- pow(sigma, -2)
sigma ~ dunif(0,100)
}
"
dados <- list(x=x, n=n)
inits.gen <- list(mu=1, sigma=100)
param<- c("mu", "sigma", "phi")
runjags.options(method = "rjags") ## sets it back to run everything on just one core
### Run JAGS
### --------------------
jagsfit <- run.jags(model = model.Gauss,
monitor = param,
data = dados, inits = inits.gen,
adapt = 1, n.chains = 1, thin = 1,
burnin = 1000, sample = 10000)
plot(jagsfit)
require(runjags)
set.seed(432104)
n <- 1000
x <- rnorm(n, 0, 5)
# Implementando a função de verossimilhança:
model.Gauss.imp <-"
model {
for (i in 1:n){
loglik[i] <- -0.5*log(2*pi)-0.5*log(pow(sigma,2))-0.5*pow(sigma,-2)*pow(x[i]-mu,2)
dummy[i] ~ dpois( -loglik[i] + const )
}
mu ~ dnorm(0,.0001)
phi <- pow(sigma, -2)
sigma ~ dunif(0,100)
}
"
dados <- list(x=x, n=n, pi=pi, const=10000, dummy=rep(0, n))
inits.gen <- list(mu=1, sigma=100)
param<- c("mu", "sigma", "phi")
runjags.options(method = "runjags") ## sets it back to run everything on just one core
### Run JAGS
### --------------------
jagsfit <- run.jags(model = model.Gauss.imp,
monitor = param,
data = dados, inits = inits.gen,
adapt = 1, n.chains = 1, thin = 1,
burnin = 1, sample = 1000)
plot(jagsfit)
jagsfit$summary
par(mfrow=c(3,1))
plot(1:1000, unlist(jagsfit$mcmc[,"mu"]), type="l", ylab="mu")
plot(1:1000, unlist(jagsfit$mcmc[,"sigma"]), type="l", ylab="sigma")
plot(1:1000, unlist(jagsfit$mcmc[,"tau"]), type="l", ylab="phi")