This report presents the analysis and comparison of two linear models, A and B, aimed at predicting the Actual_Weight of objects created by a 3D printer based on their CAD_Weight. These models are evaluated using the filament1 dataset, employing various statistical methods and metrics for a comprehensive comparison.
Model Estimation and Prediction Data Loading
Model Estimation Model A and Model B are estimated from the filament1 dataset. These models differ in their variance structure, intending to capture the relationship between CAD_Weight and Actual_Weight accurately.
Prediction and Interval Calculation Using the estimated models, predictive distributions and 95% prediction intervals are computed for new datasets.
Performance Evaluation Squared Error and Dawid-Sebastiani Score Functions to calculate the Squared Error (ES) and Dawid-Sebastiani Score (DS) are utilized to assess the model predictions’ accuracy. Leave-one-out cross-validation is performed for both models, providing insights into their predictive performance on unseen data.
head(loo_a) Monte Carlo P-Value Estimate A Monte Carlo estimate of the p-value tests the exchangeability between the predictions from Models A and B.
This report presents a statistical analysis of archaeological findings from a gravesite outside the walls of Visby, Gotland. In the 1920s, excavations uncovered 493 femurs (256 left and 237 right), raising questions about the number of individuals buried at the site. Through the application of Bayesian analysis and Monte Carlo integration methods, we aim to estimate the total number of people buried and the likelihood of finding femurs.
The analysis is based on the assumption that the number of left and right femurs represents two independent observations from a Binomial distribution with unknown parameters: the total number of individuals (N) and the probability of finding a femur (ϕ). We adopt a Bayesian approach, setting prior distributions for N and ϕ to reflect historical beliefs about the gravesite.
The posterior probabilities are estimated using Monte Carlo integration, leveraging the arch_loglike function for log-likelihood calculations and the estimate function to approximate posterior expectations of N and ϕ.
Estimation of Parameters The Monte Carlo integration method was applied with parameters y = c(237, 256), xi = 1/1001, a = 0.5, b = 0.5, and K = 10000 to estimate the posterior expectations of the total number of individuals buried and the probability of finding a femur.
#' Name, UUN
#' Add your own function definitions on this file.
# Load necessary library
#' neg_log_lik
#
#' @description Evaluate the negated log-likelihood for model A and B
#' @param beta A vector with the beta parameters
#' @param data A `data.frame` with the same variables as the `filament1` data set.
#' Must have columns `CAD_Weight` and `Actual_Weight`
#' @param model Either "A" for a log-linear variance model, or "B" for a proportional
#' scaling error model
neg_log_lik <- function(beta, data, model){
mu <- beta[1] + beta[2]*data[["CAD_Weight"]]
# distinguish between the two models to find the particular standard deviation for the betas
if(model == "A") {
sigma <- sqrt(exp(beta[3] + beta[4]*data[["CAD_Weight"]]))
}else{
sigma <- sqrt(exp(beta[3])+exp(beta[4]) * (data[["CAD_Weight"]]^2))
}
- sum(dnorm(data[["Actual_Weight"]],
mean = mu,
sd=sigma,
log = TRUE))
}
#' filament_estimate
#
#' @description Estimate filament models with different variance structure
#' @param data A `data.frame` with the same variables as the `filament1` data set.
#' Must have columns `CAD_Weight` and `Actual_Weight`
#' @param model Either "A" for a log-linear variance model, or "B" for a proportional
#' scaling error model
#' @return An estimation object suitable for use with [filament1_predict()]
filament1_estimate <- function(data, model) {
model <- match.arg(model, c("A", "B"))
if (model == "A") {
beta_start <- c(-0.1, 1.07, -2, 0.05)
} else {
beta_start <- c(-0.15, 1.07, -13.5, -6.5)
}
opt <- optim(beta_start,
neg_log_lik,
data = data,
model = model,
hessian = TRUE,
method = "Nelder-Mead",
control = list(maxit = 5000)
)
fit <- list(
model = model,
par = opt$par,
hessian = opt$hessian
)
class(fit) <- c("filament1_estimate", "list")
fit
}
#' filament1_aux_EV
#'
#' @description Evaluate the expectation and variance for model A and B
#' @param beta A vector with the beta parameters
#' @param data A `data.frame` containing the required predictors, including `CAD_Weight`
#' @param model Either "A" for a log-linear variance model, or "B" for a proportional
#' scaling error model
#' @param Sigma_beta : If not NULL, an estimate of the covariance matrix for
# the uncertainty of estimated betas
#' @return A list with four elements:
# E : E(y|beta,x)
# V : Var(y|beta,x)
# VE : Var(E(y|beta,x)|x) or NULL
# EV : E(Var(y|beta,x)|x) or NULL
filament1_aux_EV <- function(beta, data, model = c("A", "B"),
Sigma_beta = NULL) {
model <- match.arg(model)
if (model == "A") {
ZE.0 <- model.matrix( ~ 1 + CAD_Weight, data = data)
ZV.0 <- model.matrix( ~ 1 + CAD_Weight, data = data)
ZE = cbind(ZE.0, ZV.0 * 0)
ZV = cbind(ZE.0 * 0, ZV.0)
VE <- EV <- NULL
if (!is.null(Sigma_beta)) {
# E(Var(y|beta,x)|x)
EV <- exp(ZV %*% beta + rowSums(ZV * (ZV %*% Sigma_beta)) / 2)
# Var(E(y|beta,x)|x)
VE <- rowSums(ZE * (ZE %*% Sigma_beta))
}
out <- list(
E = ZE %*% beta,
V = exp(ZV %*% beta),
VE = VE,
EV = EV
)
} else {
ZE.0 <- model.matrix( ~ 1 + CAD_Weight, data = data)
ZV.0 <- model.matrix( ~ 1 + I(CAD_Weight^2), data = data)
ZE = cbind(ZE.0, ZV.0 * 0)
ZV = cbind(ZE.0 * 0, ZV.0)
VE <- EV <- NULL
if (!is.null(Sigma_beta)) {
# E(Var(y|beta,x)|x)
# (pmin: Ignore large Sigma_beta values)
EV <- ZV %*% exp(beta + pmin(0.5^2, diag(Sigma_beta)) / 2)
# Var(E(y|beta,x)|x)
VE <- rowSums(ZE * (ZE %*% Sigma_beta))
}
out <- list(
E = ZE %*% beta,
V = ZV %*% exp(beta),
VE = VE,
EV = EV
)
}
out
}
#' filament1_predict
#'
#' @description Computes predictive distributions and 95% prediction intervals for a new dataset
#' @param newdata A `data.frame` with the same variables as the `filament1` dataset
#' @param fit An object returned by `filament1_estimate()`
#' @return A data.frame with variables mean, sd, lwr, and upr
filament1_predict <- function(newdata, fit) {
ev <- filament1_aux_EV(fit$par, newdata, model = fit$model)
z <- qnorm(0.975) # 95% prediction interval
mean <- ev$E
sd <- sqrt(ev$V)
lwr <- mean - z * sd
upr <- mean + z * sd
return(data.frame(mean = mean, sd = sd, lwr = lwr, upr = upr))
}
#' Squared Error (ES)
#'
#' @description Calculates the squared error between actual and predicted values
#' @param actual The actual values
#' @param predicted The predicted values
#' @return The squared error
ES <- function(actual, predicted) {
sum((actual - predicted)^2)
}
#' Dawid-Sebastiani Score (DS)
#'
#' @description Calculates the Dawid-Sebastiani score
#' @param actual The actual values
#' @param predicted The predicted mean values
#' @param sd The standard deviation of the predicted values
#' @return The Dawid-Sebastiani score
DS <- function(actual, predicted, sd) {
sum(((actual - predicted)/sd)^2 - 2 * log(sd))
}
#' leave1out
#'
#' @description Performs leave-one-out cross-validation for the selected model
#' @param data A `data.frame` with the same variables as the `filament1` data set.
#' @param model Either "A" for a log-linear variance model, or "B" for a proportional scaling error model
#' @return A data.frame that extends the original data frame with four additional columns:
#' mean, sd, se, and ds for leave-one-out prediction means, standard deviations, and prediction scores.
leave1out <- function(data, model) {
n <- nrow(data)
predictions <- data.frame(mean=numeric(n), sd=numeric(n), se=numeric(n), ds=numeric(n))
for(i in 1:n) {
train_data <- data[-i, ]
test_data <- data[i, , drop = FALSE]
fit <- filament1_estimate(train_data, model)
pred <- filament1_predict(test_data, fit)
predictions$mean[i] <- pred$mean
predictions$sd[i] <- pred$sd
actual <- test_data$Actual_Weight
predicted <- pred$mean
predictions$se[i] <- ES(actual, predicted)
predictions$ds[i] <- DS(actual, predicted, pred$sd)
}
return(cbind(data, predictions))
}
#' Monte Carlo P-value Estimate
#'
#' @description Constructs a Monte Carlo estimate of the p-value to test the exchangeability
#' between model predictions from A and B against the alternative hypothesis that B is better than A.
#' @param data A `data.frame` with the same variables as the `filament1` data set.
#' @return A list with p-value and Monte Carlo standard errors
monte_carlo_pvalue <- function(data) {
set.seed(123) # Ensure reproducibility
n_sims <- 1000 # Reduced for quicker execution, adjust as needed
diff_scores <- numeric(n_sims)
for(i in 1:n_sims) {
shuffled_data <- data
shuffled_data$Actual_Weight <- sample(data$Actual_Weight)
loo_a <- leave1out(shuffled_data, "A")
loo_b <- leave1out(shuffled_data, "B")
score_a <- mean(loo_a$ds)
score_b <- mean(loo_b$ds)
diff_scores[i] <- score_b - score_a
if (i %% 100 == 0) { # Print progress every 100 iterations
cat("Iteration:", i, "Difference:", score_b - score_a, "\n")
}
}
p_value <- mean(diff_scores < 0)
mc_se <- sd(diff_scores < 0) / sqrt(n_sims)
return(list(p_value = p_value, mc_se = mc_se))
}
library(dplyr)
# Define the log-likelihood function
arch_loglike <- function(N, phi, y1 = 256, y2 = 237) {
l <- -lgamma(y1 + 1) - lgamma(y2 + 1) -
lgamma(N - y1 + 1) - lgamma(N - y2 + 1) +
2 * lgamma(N + 1) + (y1 + y2) * log(phi) + (2 * N - y1 - y2) * log(1 - phi)
return(l)
}
# Function to sample from the posterior distribution using Monte Carlo integration
monte_carlo_integration <- function(K = 10000) {
xi <- 1 / (1 + 1000)
a <- 2
b <- 2
n_samples <- rgeom(K, xi)
phi_samples <- rbeta(K, a, b)
log_likelihoods <- sapply(1:K, function(k) arch_loglike(n_samples[k], phi_samples[k]))
p_y_hat <- mean(exp(log_likelihoods))
E_N_hat <- sum(n_samples * exp(log_likelihoods)) / (K * p_y_hat)
E_phi_hat <- sum(phi_samples * exp(log_likelihoods)) / (K * p_y_hat)
return(list(E_N = E_N_hat, E_phi = E_phi_hat))
}
library(dplyr)
# Function to calculate the combined log-likelihood
arch_loglike <- function(N, phi, y) {
y1 <- y[1]
y2 <- y[2]
l <- -lgamma(y1 + 1) - lgamma(y2 + 1) -
lgamma(N - y1 + 1) - lgamma(N - y2 + 1) +
2 * lgamma(N + 1) + (y1 + y2) * log(phi) + (2 * N - y1 - y2) * log(1 - phi)
return(l)
}
# Monte Carlo integration function to estimate py(y), E(N|y), and E(ϕ|y)
estimate <- function(y, xi, a, b, K) {
n_samples <- rgeom(K, xi)
phi_samples <- rbeta(K, a, b)
log_likelihoods <- sapply(1:K, function(k) arch_loglike(n_samples[k], phi_samples[k], y))
p_y_hat <- mean(exp(log_likelihoods))
E_N_hat <- sum(n_samples * exp(log_likelihoods)) / (K * p_y_hat)
E_phi_hat <- sum(phi_samples * exp(log_likelihoods)) / (K * p_y_hat)
return(list(p_y = p_y_hat, E_N = E_N_hat, E_phi = E_phi_hat))
}