-
Notifications
You must be signed in to change notification settings - Fork 0
/
pROBLEM 2
94 lines (75 loc) · 3.04 KB
/
pROBLEM 2
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
Consider a single observation (y1, y2) = (0,0) from a bivariate normally distributed population
(θ1, θ2, 1, 1, 0.8). Assume the prior distribution of θ1 and θ2 to be uniform. Apply Gibb’s sampler to generate samples on θ1 and θ2 from their posterior
distributions. Compute their posterior means and variances and show their posterior distributions.
# Define the likelihood function
likelihood_function <- function(y, theta1, theta2) {
dnorm(y, mean = c(theta1, theta2), sd = c(1, 1), log = TRUE)
}
# Define the prior distributions
prior_distribution_theta1 <- function(theta1) {
if (theta1 > 0 && theta1 < 1) {
1
} else {
0
}
}
prior_distribution_theta2 <- function(theta2) {
if (theta2 > 0 && theta2 < 1) {
1
} else {
0
}
}
# Define the conditional distributions
conditional_distribution_theta1 <- function(theta1, theta2, y) {
likelihood_function(y, theta1, theta2) * prior_distribution_theta1(theta1)
}
conditional_distribution_theta2 <- function(theta1, theta2, y) {
likelihood_function(y, theta1, theta2) * prior_distribution_theta2(theta2)
}
# Implement Gibbs sampling
gibbs_sampler <- function(y, num_samples) {
theta1_samples <- rep(NA, num_samples)
theta2_samples <- rep(NA, num_samples)
# Initialize the starting values
theta1 <- runif(1)
theta2 <- runif(1)
for (i in 2:num_samples) {
# Sample theta1 from the conditional distribution
theta1_proposal <- rnorm(1, mean = theta1, sd = 1)
alpha <- conditional_distribution_theta1(theta1_proposal, theta2, y) / conditional_distribution_theta1(theta1, theta2, y)
if (runif(1) < alpha) {
theta1 <- theta1_proposal
}
# Sample theta2 from the conditional distribution
theta2_proposal <- rnorm(1, mean = theta2, sd = 1)
alpha <- conditional_distribution_theta2(theta1, theta2_proposal, y) / conditional_distribution_theta2(theta1, theta2, y)
if (runif(1) < alpha) {
theta2 <- theta2_proposal
}
# Store the samples
theta1_samples[i] <- theta1
theta2_samples[i] <- theta2
}
return(theta1_samples, theta2_samples)
}
# Set the parameters
y <- c(0, 0)
num_samples <- 10000
# Run the Gibbs sampler
theta1_samples <- gibbs_sampler(y, num_samples)[1, ]
theta2_samples <- gibbs_sampler(y, num_samples)[2, ]
# Compute the posterior means
posterior_mean_theta1 <- mean(theta1_samples)
posterior_mean_theta2 <- mean(theta2_samples)
# Compute the posterior variances
posterior_variance_theta1 <- var(theta1_samples)
posterior_variance_theta2 <- var(theta2_samples)
# Plot the posterior distributions
hist(theta1_samples, main = "Posterior Distribution of theta1", xlab = "theta1", ylab = "Density")
hist(theta2_samples, main = "Posterior Distribution of theta2", xlab = "theta2", ylab = "Density")
# Print the posterior means and variances
print(paste("Posterior mean of theta1:", posterior_mean_theta1))
print(paste("Posterior variance of theta1:", posterior_variance_theta1))
print(paste("Posterior mean of theta2:", posterior_mean_theta2))
print(paste("Posterior variance of theta2:", posterior_variance_theta2))