CS计算机代考程序代写 chain Bayesian # Regression example, taken from “A First Course in Bayesian Statistical

# Regression example, taken from “A First Course in Bayesian Statistical
# Methods” by Peter D. Hoff

# Arizona diabetes data

rm(list=ls()); set.seed(20210325);

# Read in data

file <- "http://www2.stat.duke.edu/~pdh10/FCBS/Exercises/azdiabetes.dat" Data <- read.table(file=file, header=T); rm(file); dim(Data); names(Data); # Data pre-processing y <- Data$glu; X <- as.matrix(Data[,c(1,3,4,5,6,7)]); rm(Data); names <- colnames(X); m <- apply(X, 2, mean); s <- apply(X, 2, sd); X <- sweep(X, 2, m, "-"); rm(m); X <- sweep(X, 2, s, "/"); rm(s); apply(X, 2, mean); apply(X, 2, sd); X <- cbind(1, X); colnames(X) <- c("ones", names); n <- dim(X)[1]; k <- dim(X)[2]; # We assume the standard noninformative prior on the variance, # p(beta, sigma2) = const*(sigma2)^(-1) # This gives the marginal posterior # sigma2 | y ~ Inv-Chisq(n-k, s2) # where s2 = 1/(n-k) * (y - X %*% beta.hat)' (y - X %*% beta.hat) # and # beta | sigma2, y ~ Normal with mean vector and covariance matrix: # (X'X)^(-1) %*% X' y & sigma2 * (X'X)^(-1) library(rstan) stan_model <- " data{ int n;
int k;
vector [n] y;
matrix[n,k] X;
# vector[k] beta0;
# matrix[k,k] Omega;
}

parameters{
vector[k] beta;
# real sigma2;
real log_sigma;
}

transformed parameters{
vector[n] mu;
real sigma;
mu = X * beta;
sigma <- exp(log_sigma); # sqrt(sigma2); } model{ # beta ~ multi_normal(beta0, sigma2*Omega) ; # sigma2 ~ inv_gamma(.001, .001); y ~ normal(mu, sigma) ; } " beta0 <- rep(0, k) data <- list(n=n, k=k, y=y, X=X) # , beta0=beta0, Omega=Omega) fit <- stan(model_code=stan_model, data=data, iter=2000, chains=5) Results <- summary(fit)$summary[c(1:7, 541), ] Results round(Results[, 4:8], 2) sims <- extract(fit)