STAT-431-EX-18_merged

pdf

School

University of Illinois, Urbana Champaign *

*We aren’t endorsed by this school

Course

431

Subject

Statistics

Date

Apr 3, 2024

Type

pdf

Pages

Uploaded by rohaanraheja11

STAT 431 EX 18 2023-12-07 Question 1 a) library (rjags) ## Loading required package: coda ## Linked to JAGS 4.3.2 ## Loaded modules: basemod,bugs # Data and initial values data <- list ( x = c ( 1997 , 1998 , 1999 , 2000 , 2001 , 2002 , 2003 , 2004 , 2005 , 2006 , 2007 , 2008 , 2009 , 2010 ), y = c ( 2.2952 , 2.3435 , 2.5512 , 2.5531 , 2.3918 , 2.1546 , 2.3596 , 2.2431 , 2.1725 , 2.3162 , 2.3504 , 2.1926 , N = 14 ) inits <- list ( list ( beta0 = 0 , beta1 = 0 , tausq = 1 ), list ( beta0 = 1 , beta1 = - 1 , tausq = 0.5 ), list ( beta0 = - 1 , beta1 = 1 , tausq = 2 ) ) # Write the model in a string model_string <- " model { for (i in 1:N) { mu[i] <- beta0 + beta1 * x[i] y[i] ~ dnorm(mu[i], tausq) } beta0 ~ dnorm(0, 1E-6) beta1 ~ dnorm(0, 1E-6) tausq ~ dgamma(0.001, 0.001) sigma <- 1 / sqrt(tausq) } " # Save the model to a file writeLines (model_string, "regression_model.txt" ) # Run the model jags_model <- jags.model ( "regression_model.txt" , data = data, inits = inits, n.chains = 3 , n.adapt = 100 ## Compiling model graph ## Resolving undeclared variables 1

## Allocating nodes ## Graph information: ## Observed stochastic nodes: 14 ## Unobserved stochastic nodes: 3 ## Total graph size: 66 ## ## Initializing model jags_samples <- coda.samples (jags_model, c ( "beta0" , "beta1" , "tausq" ), n.iter = 5000 ) # Print summary and diagnostic plots summary (jags_samples) ## ## Iterations = 1:5000 ## Thinning interval = 1 ## Number of chains = 3 ## Sample size per chain = 5000 ## ## 1. Empirical mean and standard deviation for each variable, ## plus standard error of the mean: ## ## Mean SD Naive SE Time-series SE ## beta0 16.708529 1614.4090 13.181594 18.712826 ## beta1 -0.007201 0.8058 0.006579 0.009318 ## tausq 14.051884 21.9573 0.179281 0.082953 ## ## 2. Quantiles for each variable: ## ## 2.5% 25% 50% 75% 97.5% ## beta0 -2.029e+03 -1.907e+03 4.942637 1948.1254 2078.425 ## beta1 -1.036e+00 -9.712e-01 -0.001322 0.9529 1.014 ## tausq 2.403e-02 4.842e-02 0.072412 29.9570 68.257 plot (jags_samples) 2

0 1000 2000 3000 4000 5000 -2000 Iterations Trace of beta0 -3000 -1000 0 1000 2000 3000 0e+00 Density of beta0 N = 5000 Bandwidth = 250.1 0 1000 2000 3000 4000 5000 -1.0 Iterations Trace of beta1 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 0.0 1.0 Density of beta1 N = 5000 Bandwidth = 0.1248 0 1000 2000 3000 4000 5000 0 80 Iterations Trace of tausq 0 20 40 60 80 100 120 0.00 Density of tausq N = 5000 Bandwidth = 3.401 a) continued # Assuming ' jags_model ' is your JAGS model object # Load the coda package library (coda) # Extract the coda object from the JAGS model coda_samples <- coda.samples (jags_model, c ( "beta0" , "beta1" , "tausq" ), n.iter = 5000 ) # Check Gelman-Rubin diagnostic gelman.diag (coda_samples) ## Potential scale reduction factors: ## ## Point est. Upper C.I. ## beta0 78.64 198.5 ## beta1 78.65 198.6 ## tausq 3.96 23.5 ## ## Multivariate psrf ## ## 60.9 # Plot trace plots for each parameter traceplot (coda_samples) 3

Your preview ends here

Eager to read complete document? Join bartleby learn and gain access to the full version

Access to all documents
Unlimited textbook solutions
24/7 expert homework help

5000 6000 7000 8000 9000 10000 -2000 -1000 0 1000 2000 Iterations Trace of beta0 5000 6000 -1.0 -0.5 0.0 0.5 1.0 5000 6000 7000 8000 9000 10000 0 20 40 60 80 100 120 Iterations Trace of tausq It converges. b) # Frequentist (lm) model lm_model <- lm (y ~ x, data = data) summary (lm_model) 4

## ## Call: ## lm(formula = y ~ x, data = data) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.20180 -0.08265 0.02306 0.07616 0.18042 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 48.38982 17.03075 2.841 0.0149 * ## x -0.02301 0.00850 -2.707 0.0191 * ## --- ## Signif. codes: 0 ' *** ' 0.001 ' ** ' 0.01 ' * ' 0.05 ' . ' 0.1 ' ' 1 ## ## Residual standard error: 0.1282 on 12 degrees of freedom ## Multiple R-squared: 0.3791, Adjusted R-squared: 0.3273 ## F-statistic: 7.326 on 1 and 12 DF, p-value: 0.01907 # Bayesian (JAGS) model summary (jags_samples) ## ## Iterations = 1:5000 ## Thinning interval = 1 ## Number of chains = 3 ## Sample size per chain = 5000 ## ## 1. Empirical mean and standard deviation for each variable, ## plus standard error of the mean: ## ## Mean SD Naive SE Time-series SE ## beta0 16.708529 1614.4090 13.181594 18.712826 ## beta1 -0.007201 0.8058 0.006579 0.009318 ## tausq 14.051884 21.9573 0.179281 0.082953 ## ## 2. Quantiles for each variable: ## ## 2.5% 25% 50% 75% 97.5% ## beta0 -2.029e+03 -1.907e+03 4.942637 1948.1254 2078.425 ## beta1 -1.036e+00 -9.712e-01 -0.001322 0.9529 1.014 ## tausq 2.403e-02 4.842e-02 0.072412 29.9570 68.257 jags_samples <- coda.samples (jags_model, c ( "beta0" , "beta1" , "tausq" ), n.iter = 5000 ) # Assuming ' jags_model ' is the object containing the JAGS model # Extract the relevant samples for beta1 jags_samples <- coda.samples (jags_model, c ( "beta1" ), n.iter = 5000 ) beta1_samples <- as.matrix (jags_samples[[ 1 ]]) # Calculate the posterior probability that beta1 > 0 posterior_prob_beta1_positive <- mean (beta1_samples > 0 , na.rm = TRUE ) # Print the result cat ( "Posterior probability that beta1 > 0:" , posterior_prob_beta1_positive, " \n " ) 5

## Posterior probability that beta1 > 0: 0.063 2. library (rjags) # Load the data data <- list ( x = c ( 32 , 33 , 34 , 35 , 36 , 37 , 38 , 39 , 40 , 41 ), y = c ( 2.4 , 2.7 , 3.2 , 2.7 , 3.0 , 3.6 , 3.4 , 4.0 , 3.7 , 3.9 ), N = 10 , prior_intercept_mean = 2 , prior_intercept_precision = 1 / ( 0.5 ˆ 2 ), # Variance = 0.5ˆ2 prior_slope_mean = 0 , prior_slope_precision = 1 / ( 0.5 ˆ 2 ) # Variance = 0.5ˆ2 ) # Specify the model in JAGS language model_string <- " model { for (i in 1:N) { mu[i] <- beta0 + beta1 * x[i] y[i] ~ dnorm(mu[i], tau) } beta0 ~ dnorm(prior_intercept_mean, prior_intercept_precision) beta1 ~ dnorm(prior_slope_mean, prior_slope_precision) tau ~ dgamma(0.001, 0.001) sigma <- 1 / sqrt(tau) } " # Save the model to a file writeLines (model_string, "informative_priors_model.txt" ) # Set initial values for each chain inits <- list ( list ( beta0 = 2 , beta1 = 0 , tau = 1 ), list ( beta0 = 1 , beta1 = - 1 , tau = 0.5 ), list ( beta0 = - 1 , beta1 = 1 , tau = 2 ) ) # Run the model jags_model_informative_priors <- jags.model ( "informative_priors_model.txt" , data = data, inits = inits, ## Compiling model graph ## Resolving undeclared variables ## Allocating nodes ## Graph information: ## Observed stochastic nodes: 10 ## Unobserved stochastic nodes: 3 ## Total graph size: 52 ## ## Initializing model 6

Your preview ends here

Eager to read complete document? Join bartleby learn and gain access to the full version

Access to all documents
Unlimited textbook solutions
24/7 expert homework help

jags_samples_informative_priors <- coda.samples (jags_model_informative_priors, c ( "beta0" , "beta1" , "tau" # Print summary and diagnostic plots summary (jags_samples_informative_priors) ## ## Iterations = 1:5000 ## Thinning interval = 1 ## Number of chains = 3 ## Sample size per chain = 5000 ## ## 1. Empirical mean and standard deviation for each variable, ## plus standard error of the mean: ## ## Mean SD Naive SE Time-series SE ## beta0 1.64340 0.58729 0.0047952 0.0203298 ## beta1 0.04501 0.01655 0.0001351 0.0005759 ## tau 5.44212 2.77889 0.0226895 0.0391238 ## ## 2. Quantiles for each variable: ## ## 2.5% 25% 50% 75% 97.5% ## beta0 0.65811 1.30762 1.64319 1.98811 2.62834 ## beta1 0.01698 0.03528 0.04512 0.05463 0.07274 ## tau 1.55397 3.46147 4.94423 6.87909 12.12741 plot (jags_samples_informative_priors) 0 1000 2000 3000 4000 5000 -20 Iterations Trace of beta0 -20 -10 0 10 20 0.0 0.8 Density of beta0 N = 5000 Bandwidth = 0.07867 0 1000 2000 3000 4000 5000 -0.4 Iterations Trace of beta1 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 0 20 Density of beta1 N = 5000 Bandwidth = 0.002237 0 1000 2000 3000 4000 5000 0 30 Iterations Trace of tau 0 10 20 30 40 0.00 Density of tau N = 5000 Bandwidth = 0.3951 7

\ ' o..) . In R ~ --- - --- ,__ ) _ __e 'a _ }:} _~ o tk 'r0 __2_ AA ~ - v-o ,r , OU S ~ ~ l CVt' } __ ) ___ '(O~ _L ____ - ) ~ b e ~ Cc:) r- s. ·; s. ~ t " +e--r rn S of c;L· v..Q ~ o " q " ~ s i 5\ rv: .C,- C..O.t'\c L f 0f (Oc-t.f:i · cA~"6 · O ,- ---- j ~ - ~ - ~ ------- ~._~_ --- - ~ ~- ~ 1' --- ... \1 -- -. - "'W ----- --- ---- 1 1,