# Define the log-likelihood function for linear regression
<- function(beta, X, y, sigma_sq) {
log_likelihood <- X %*% beta
y_hat <- y - y_hat
residuals <- -0.5 * length(y) * log(2 * pi * sigma_sq) - 0.5 * sum(residuals^2) / sigma_sq
log_likelihood return(log_likelihood)
}
# Define the prior distribution for beta
<- function(beta) {
prior_beta <- rep(0, length(beta))
prior_mean <- rep(10, length(beta))
prior_sd <- sum(dnorm(beta, mean = prior_mean, sd = prior_sd, log = TRUE))
log_prior return(log_prior)
}
# Define the prior distribution for sigma
<- function(sigma_sq) {
prior_sigma <- 2
alpha <- 2
beta # log_prior <- dgamma(1/sigma_sq, shape = alpha, rate = beta, log = TRUE)
<- extraDistr::dinvgamma(sigma_sq, alpha = alpha, beta = beta, log = TRUE)
log_prior
return(log_prior)
}
# Define the proposal distribution for beta
<- function(beta, scale) {
proposal_beta <- rnorm(length(beta), mean = beta, sd = scale)
beta_proposal return(beta_proposal)
}
# Define the proposal distribution for sigma
<- function(sigma_sq, scale) {
proposal_sigma # sigma_proposal <- rgamma(1, shape = sigma_sq / scale, rate = scale)
<- extraDistr::rinvgamma(1, alpha = sigma_sq / scale, beta = scale)
sigma_proposal return(sigma_proposal)
}
# Set up the data
# set.seed(123)
# n <- 100
# X <- cbind(1, rnorm(n), rnorm(n), rnorm(n))
# beta_true <- c(1, 2, 3, 4)/4
# sigma_true <- 1
# y <- X %*% beta_true + rnorm(n, sd = sigma_true)
# Set up the Metropolis-Hastings algorithm
# n_iter <- 10000
# Run the Metropolis-Hastings algorithm
= function(
mh
X,
y,beta = rep(0, ncol(X)),
sigma_sq = .5,
scale_beta = 0.1,
scale_sigma = 1,
chains = 2,
warmup = 1000,
n_iter = 2000,
seed = 123
) {set.seed(seed)
<- list()
result <- beta
beta_start <- sigma_sq
sigma_sq_start
for (c in 1:chains){
<- 0
acceptance_beta <- 0
acceptance_sigma <- matrix(0, n_iter, ncol(X))
beta_samples <- rep(0, n_iter)
sigma_sq_samples
if (c > 1) {
<- beta_start
beta <- sigma_sq_start
sigma_sq
}
for (i in 1:n_iter) {
# Update beta
<- proposal_beta(beta, scale_beta)
beta_proposal <- log_likelihood(beta_proposal, X, y, sigma_sq) + prior_beta(beta_proposal) -
log_ratio_beta log_likelihood(beta, X, y, sigma_sq) - prior_beta(beta)
if (log(runif(1)) < log_ratio_beta) {
<- beta_proposal
beta <- acceptance_beta + 1
acceptance_beta
}<- beta
beta_samples[i, ]
# Update sigma_sq
<- proposal_sigma(sigma_sq, scale_sigma)
sigma_sq_proposal <- log_likelihood(beta, X, y, sigma_sq_proposal) + prior_sigma(sigma_sq_proposal) -
log_ratio_sigma log_likelihood(beta, X, y, sigma_sq) - prior_sigma(sigma_sq)
if (log(runif(1)) < log_ratio_sigma) {
<- sigma_sq_proposal
sigma_sq <- acceptance_sigma + 1
acceptance_sigma
}<- sigma_sq
sigma_sq_samples[i]
}
message("Acceptance rate for beta:", acceptance_beta / n_iter, "\n")
message("Acceptance rate for sigma:", acceptance_sigma / n_iter, "\n")
= list(
result[[c]] beta = beta_samples[-(1:warmup), ],
sigma_sq = sigma_sq_samples[-(1:warmup)],
# y_rep = X %*% t(beta_samples[-(1:warmup), ])
# +rnorm(n_iter - warmup, sd = sqrt(sigma_sq_samples[-(1:warmup)]))
y_rep = t(X %*% t(beta_samples[-(1:warmup), ]) + rnorm(n_iter - warmup, sd = sqrt(sigma_sq_samples[-(1:warmup)])))
)
}
result
}
= df_happiness |>
X_train select(life_exp, gdp_pc, corrupt) |>
as.matrix()
= mh(
our_result X = cbind(1, X_train),
y = df_happiness$happiness,
beta = c(mean(df_happiness$happiness), rep(0, ncol(X_train))),
sigma_sq = var(df_happiness$happiness),
scale_sigma = .5,
warmup = 1000,
n_iter = 2000
)
str(our_result)
Appendix E — Other to come
E.1 Simulation
E.2 Bayesian Demonstration
Metropolis-Hastings demo
reference nice shiny app https://github.com/tomicapretto/shiny-hmc
E.3 Linear Programming
E.4 Boosting
import boosting_demo.R