# Load packages
library(bayesrules)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(rstan)
## Loading required package: StanHeaders
##
## rstan version 2.32.6 (Stan version 2.32.2)
##
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
##
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
##
## Attaching package: 'rstan'
##
## The following object is masked from 'package:tidyr':
##
## extract
library(rstanarm)
## Loading required package: Rcpp
## This is rstanarm version 2.32.1
## - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
## - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
## - For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores())
##
## Attaching package: 'rstanarm'
##
## The following object is masked from 'package:rstan':
##
## loo
library(bayesplot)
## This is bayesplot version 1.11.1
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
## * Does _not_ affect other ggplot2 plots
## * See ?bayesplot_theme_set for details on theme setting
library(tidybayes)
library(janitor)
##
## Attaching package: 'janitor'
##
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
library(broom.mixed)
library(ggplot2)
library(patchwork)
# Load and plot data
data(bikes)
head(bikes)
## date season year month day_of_week weekend holiday temp_actual
## 1 2011-01-01 winter 2011 Jan Sat TRUE no 57.39952
## 2 2011-01-03 winter 2011 Jan Mon FALSE no 46.49166
## 3 2011-01-04 winter 2011 Jan Tue FALSE no 46.76000
## 4 2011-01-05 winter 2011 Jan Wed FALSE no 48.74943
## 5 2011-01-07 winter 2011 Jan Fri FALSE no 46.50332
## 6 2011-01-08 winter 2011 Jan Sat TRUE no 44.17700
## temp_feel humidity windspeed weather_cat rides
## 1 64.72625 80.5833 10.74988 categ2 654
## 2 49.04645 43.7273 16.63670 categ1 1229
## 3 51.09098 59.0435 10.73983 categ1 1454
## 4 52.63430 43.6957 12.52230 categ1 1518
## 5 50.79551 49.8696 11.30464 categ2 1362
## 6 46.60286 53.5833 17.87587 categ2 891
# Excercise 9.9
#a)
plot1 <- plot_normal(mean = 5000, sd = 1400)+
labs(x = "beta_0c", y = "pdf")
plot2 <- plot_normal(mean =-10, sd = 3)+
labs(x = "beta_1", y = "pdf")
plot3 <- plot_gamma(shape = 1, rate = 0.0005)+
labs(x = "sigma", y = "pdf")
plot1+plot2+plot3
#b)
# 9.3.1 Simulation via rstanarm
bike_prior_model <- stan_glm(rides ~ humidity, data = bikes,
family = gaussian,
prior_intercept = normal(5000, 1400),
prior = normal(-10,3),
prior_aux = exponential(0.0005),
chains = 5, iter = 4000*2, seed = 8473,prior_PD = TRUE)
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000837 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 8.37 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 8000 [ 0%] (Warmup)
## Chain 1: Iteration: 800 / 8000 [ 10%] (Warmup)
## Chain 1: Iteration: 1600 / 8000 [ 20%] (Warmup)
## Chain 1: Iteration: 2400 / 8000 [ 30%] (Warmup)
## Chain 1: Iteration: 3200 / 8000 [ 40%] (Warmup)
## Chain 1: Iteration: 4000 / 8000 [ 50%] (Warmup)
## Chain 1: Iteration: 4001 / 8000 [ 50%] (Sampling)
## Chain 1: Iteration: 4800 / 8000 [ 60%] (Sampling)
## Chain 1: Iteration: 5600 / 8000 [ 70%] (Sampling)
## Chain 1: Iteration: 6400 / 8000 [ 80%] (Sampling)
## Chain 1: Iteration: 7200 / 8000 [ 90%] (Sampling)
## Chain 1: Iteration: 8000 / 8000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.45 seconds (Warm-up)
## Chain 1: 0.113 seconds (Sampling)
## Chain 1: 0.563 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 8000 [ 0%] (Warmup)
## Chain 2: Iteration: 800 / 8000 [ 10%] (Warmup)
## Chain 2: Iteration: 1600 / 8000 [ 20%] (Warmup)
## Chain 2: Iteration: 2400 / 8000 [ 30%] (Warmup)
## Chain 2: Iteration: 3200 / 8000 [ 40%] (Warmup)
## Chain 2: Iteration: 4000 / 8000 [ 50%] (Warmup)
## Chain 2: Iteration: 4001 / 8000 [ 50%] (Sampling)
## Chain 2: Iteration: 4800 / 8000 [ 60%] (Sampling)
## Chain 2: Iteration: 5600 / 8000 [ 70%] (Sampling)
## Chain 2: Iteration: 6400 / 8000 [ 80%] (Sampling)
## Chain 2: Iteration: 7200 / 8000 [ 90%] (Sampling)
## Chain 2: Iteration: 8000 / 8000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.43 seconds (Warm-up)
## Chain 2: 0.15 seconds (Sampling)
## Chain 2: 0.58 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.7e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 8000 [ 0%] (Warmup)
## Chain 3: Iteration: 800 / 8000 [ 10%] (Warmup)
## Chain 3: Iteration: 1600 / 8000 [ 20%] (Warmup)
## Chain 3: Iteration: 2400 / 8000 [ 30%] (Warmup)
## Chain 3: Iteration: 3200 / 8000 [ 40%] (Warmup)
## Chain 3: Iteration: 4000 / 8000 [ 50%] (Warmup)
## Chain 3: Iteration: 4001 / 8000 [ 50%] (Sampling)
## Chain 3: Iteration: 4800 / 8000 [ 60%] (Sampling)
## Chain 3: Iteration: 5600 / 8000 [ 70%] (Sampling)
## Chain 3: Iteration: 6400 / 8000 [ 80%] (Sampling)
## Chain 3: Iteration: 7200 / 8000 [ 90%] (Sampling)
## Chain 3: Iteration: 8000 / 8000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.404 seconds (Warm-up)
## Chain 3: 0.105 seconds (Sampling)
## Chain 3: 0.509 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 8e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 8000 [ 0%] (Warmup)
## Chain 4: Iteration: 800 / 8000 [ 10%] (Warmup)
## Chain 4: Iteration: 1600 / 8000 [ 20%] (Warmup)
## Chain 4: Iteration: 2400 / 8000 [ 30%] (Warmup)
## Chain 4: Iteration: 3200 / 8000 [ 40%] (Warmup)
## Chain 4: Iteration: 4000 / 8000 [ 50%] (Warmup)
## Chain 4: Iteration: 4001 / 8000 [ 50%] (Sampling)
## Chain 4: Iteration: 4800 / 8000 [ 60%] (Sampling)
## Chain 4: Iteration: 5600 / 8000 [ 70%] (Sampling)
## Chain 4: Iteration: 6400 / 8000 [ 80%] (Sampling)
## Chain 4: Iteration: 7200 / 8000 [ 90%] (Sampling)
## Chain 4: Iteration: 8000 / 8000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.441 seconds (Warm-up)
## Chain 4: 0.164 seconds (Sampling)
## Chain 4: 0.605 seconds (Total)
## Chain 4:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 5).
## Chain 5:
## Chain 5: Gradient evaluation took 1e-05 seconds
## Chain 5: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 5: Adjust your expectations accordingly!
## Chain 5:
## Chain 5:
## Chain 5: Iteration: 1 / 8000 [ 0%] (Warmup)
## Chain 5: Iteration: 800 / 8000 [ 10%] (Warmup)
## Chain 5: Iteration: 1600 / 8000 [ 20%] (Warmup)
## Chain 5: Iteration: 2400 / 8000 [ 30%] (Warmup)
## Chain 5: Iteration: 3200 / 8000 [ 40%] (Warmup)
## Chain 5: Iteration: 4000 / 8000 [ 50%] (Warmup)
## Chain 5: Iteration: 4001 / 8000 [ 50%] (Sampling)
## Chain 5: Iteration: 4800 / 8000 [ 60%] (Sampling)
## Chain 5: Iteration: 5600 / 8000 [ 70%] (Sampling)
## Chain 5: Iteration: 6400 / 8000 [ 80%] (Sampling)
## Chain 5: Iteration: 7200 / 8000 [ 90%] (Sampling)
## Chain 5: Iteration: 8000 / 8000 [100%] (Sampling)
## Chain 5:
## Chain 5: Elapsed Time: 0.324 seconds (Warm-up)
## Chain 5: 0.105 seconds (Sampling)
## Chain 5: 0.429 seconds (Total)
## Chain 5:
# Effective sample size ratio and Rhat
neff_ratio(bike_prior_model)
## (Intercept) humidity sigma
## 0.78470 0.71245 0.93525
rhat(bike_prior_model)
## (Intercept) humidity sigma
## 1.0002557 0.9999558 0.9999925
#Autocorrelation
mcmc_acf(bike_prior_model)
# Trace plots of parallel chains
mcmc_trace(bike_prior_model, size = 0.1)
# Density plots of parallel chains
mcmc_dens_overlay(bike_prior_model)
#c)
# 100 simulated model lines (prior model)
set.seed(8473)
bikes %>%
add_epred_draws(bike_prior_model, ndraws = 100) %>%
ggplot(aes(x = humidity, y = rides)) +
geom_line(aes(y = .epred, group = .draw), alpha = 0.15) +
geom_point(data = bikes, size = 0.05)
# Simulate four sets of data
set.seed(8473)
bikes %>%
add_predicted_draws(bike_prior_model, ndraws = 4) %>%
ggplot(aes(x = humidity, y = rides)) +
geom_point(aes(y = .prediction, group = .draw), size = 0.2) +
facet_wrap(~ .draw)
#Exercise 9.10
ggplot(bikes, aes(x = humidity, y = rides)) +
geom_point(size = 0.5) +
geom_smooth(method = "lm", se = FALSE)
## `geom_smooth()` using formula = 'y ~ x'
#Exercise 9.11
#a)
bike_posterior_model <- update(bike_prior_model, prior_PD = FALSE)
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.00107 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.7 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 8000 [ 0%] (Warmup)
## Chain 1: Iteration: 800 / 8000 [ 10%] (Warmup)
## Chain 1: Iteration: 1600 / 8000 [ 20%] (Warmup)
## Chain 1: Iteration: 2400 / 8000 [ 30%] (Warmup)
## Chain 1: Iteration: 3200 / 8000 [ 40%] (Warmup)
## Chain 1: Iteration: 4000 / 8000 [ 50%] (Warmup)
## Chain 1: Iteration: 4001 / 8000 [ 50%] (Sampling)
## Chain 1: Iteration: 4800 / 8000 [ 60%] (Sampling)
## Chain 1: Iteration: 5600 / 8000 [ 70%] (Sampling)
## Chain 1: Iteration: 6400 / 8000 [ 80%] (Sampling)
## Chain 1: Iteration: 7200 / 8000 [ 90%] (Sampling)
## Chain 1: Iteration: 8000 / 8000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.262 seconds (Warm-up)
## Chain 1: 0.285 seconds (Sampling)
## Chain 1: 0.547 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1.3e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 8000 [ 0%] (Warmup)
## Chain 2: Iteration: 800 / 8000 [ 10%] (Warmup)
## Chain 2: Iteration: 1600 / 8000 [ 20%] (Warmup)
## Chain 2: Iteration: 2400 / 8000 [ 30%] (Warmup)
## Chain 2: Iteration: 3200 / 8000 [ 40%] (Warmup)
## Chain 2: Iteration: 4000 / 8000 [ 50%] (Warmup)
## Chain 2: Iteration: 4001 / 8000 [ 50%] (Sampling)
## Chain 2: Iteration: 4800 / 8000 [ 60%] (Sampling)
## Chain 2: Iteration: 5600 / 8000 [ 70%] (Sampling)
## Chain 2: Iteration: 6400 / 8000 [ 80%] (Sampling)
## Chain 2: Iteration: 7200 / 8000 [ 90%] (Sampling)
## Chain 2: Iteration: 8000 / 8000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.627 seconds (Warm-up)
## Chain 2: 0.271 seconds (Sampling)
## Chain 2: 0.898 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.1e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 8000 [ 0%] (Warmup)
## Chain 3: Iteration: 800 / 8000 [ 10%] (Warmup)
## Chain 3: Iteration: 1600 / 8000 [ 20%] (Warmup)
## Chain 3: Iteration: 2400 / 8000 [ 30%] (Warmup)
## Chain 3: Iteration: 3200 / 8000 [ 40%] (Warmup)
## Chain 3: Iteration: 4000 / 8000 [ 50%] (Warmup)
## Chain 3: Iteration: 4001 / 8000 [ 50%] (Sampling)
## Chain 3: Iteration: 4800 / 8000 [ 60%] (Sampling)
## Chain 3: Iteration: 5600 / 8000 [ 70%] (Sampling)
## Chain 3: Iteration: 6400 / 8000 [ 80%] (Sampling)
## Chain 3: Iteration: 7200 / 8000 [ 90%] (Sampling)
## Chain 3: Iteration: 8000 / 8000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.8 seconds (Warm-up)
## Chain 3: 0.279 seconds (Sampling)
## Chain 3: 1.079 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.3e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 8000 [ 0%] (Warmup)
## Chain 4: Iteration: 800 / 8000 [ 10%] (Warmup)
## Chain 4: Iteration: 1600 / 8000 [ 20%] (Warmup)
## Chain 4: Iteration: 2400 / 8000 [ 30%] (Warmup)
## Chain 4: Iteration: 3200 / 8000 [ 40%] (Warmup)
## Chain 4: Iteration: 4000 / 8000 [ 50%] (Warmup)
## Chain 4: Iteration: 4001 / 8000 [ 50%] (Sampling)
## Chain 4: Iteration: 4800 / 8000 [ 60%] (Sampling)
## Chain 4: Iteration: 5600 / 8000 [ 70%] (Sampling)
## Chain 4: Iteration: 6400 / 8000 [ 80%] (Sampling)
## Chain 4: Iteration: 7200 / 8000 [ 90%] (Sampling)
## Chain 4: Iteration: 8000 / 8000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.567 seconds (Warm-up)
## Chain 4: 0.267 seconds (Sampling)
## Chain 4: 0.834 seconds (Total)
## Chain 4:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 5).
## Chain 5:
## Chain 5: Gradient evaluation took 1e-05 seconds
## Chain 5: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 5: Adjust your expectations accordingly!
## Chain 5:
## Chain 5:
## Chain 5: Iteration: 1 / 8000 [ 0%] (Warmup)
## Chain 5: Iteration: 800 / 8000 [ 10%] (Warmup)
## Chain 5: Iteration: 1600 / 8000 [ 20%] (Warmup)
## Chain 5: Iteration: 2400 / 8000 [ 30%] (Warmup)
## Chain 5: Iteration: 3200 / 8000 [ 40%] (Warmup)
## Chain 5: Iteration: 4000 / 8000 [ 50%] (Warmup)
## Chain 5: Iteration: 4001 / 8000 [ 50%] (Sampling)
## Chain 5: Iteration: 4800 / 8000 [ 60%] (Sampling)
## Chain 5: Iteration: 5600 / 8000 [ 70%] (Sampling)
## Chain 5: Iteration: 6400 / 8000 [ 80%] (Sampling)
## Chain 5: Iteration: 7200 / 8000 [ 90%] (Sampling)
## Chain 5: Iteration: 8000 / 8000 [100%] (Sampling)
## Chain 5:
## Chain 5: Elapsed Time: 0.495 seconds (Warm-up)
## Chain 5: 0.258 seconds (Sampling)
## Chain 5: 0.753 seconds (Total)
## Chain 5:
#b)
# Effective sample size ratio and Rhat
neff_ratio(bike_posterior_model)
## (Intercept) humidity sigma
## 1.00565 1.04070 0.98125
rhat(bike_posterior_model)
## (Intercept) humidity sigma
## 1.000023 1.000030 1.000152
#Autocorrelation
mcmc_acf(bike_posterior_model)
# Trace plots of parallel chains
mcmc_trace(bike_posterior_model, size = 0.1)
# Density plots of parallel chains
mcmc_dens_overlay(bike_posterior_model)
#c)
# 100 simulated model lines (prior model)
set.seed(8473)
prior100 <- bikes %>%
add_epred_draws(bike_prior_model, ndraws = 100) %>%
ggplot(aes(x = humidity, y = rides)) +
geom_line(aes(y = .epred, group = .draw), alpha = 0.15) +
geom_point(data = bikes, size = 0.05)+
labs(title = "Simulación 100 rectas del modelo a priori")
# 100 simulated model lines (posterior model)
set.seed(8473)
post100 <- bikes %>%
add_epred_draws(bike_posterior_model, ndraws = 100) %>%
ggplot(aes(x = humidity, y = rides)) +
geom_line(aes(y = .epred, group = .draw), alpha = 0.15) +
geom_point(data = bikes, size = 0.05)+
labs(title = "Simulación 100 rectas del modelo a posteriori")
prior100 + post100
#Exercise 9.12
#a) b) y c)
# Posterior summary statistics
tidy(bike_posterior_model, effects = c("fixed", "aux"),
conf.int = TRUE, conf.level = 0.95)
## # A tibble: 4 × 5
## term estimate std.error conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 4067. 176. 3732. 4410.
## 2 humidity -9.14 2.52 -14.1 -4.30
## 3 sigma 1573. 49.7 1481. 1677.
## 4 mean_PPD 3486. 99.3 3295. 3682.
#d)
# Store the 4 chains for each parameter in 1 data frame
bike_posterior_model_df <- as.data.frame(bike_posterior_model)
head(bike_posterior_model_df)
## (Intercept) humidity sigma
## 1 4055.334 -9.104262 1618.898
## 2 4340.841 -13.985910 1659.860
## 3 4246.110 -12.485416 1652.535
## 4 4109.626 -10.774143 1616.585
## 5 3836.199 -6.050123 1570.407
## 6 3876.629 -6.112626 1519.148
# Tabulate the beta_1 values that exceed 0
bike_posterior_model_df %>%
mutate(exceeds_0 = humidity < 0) %>%
tabyl(exceeds_0)
## exceeds_0 n percent
## FALSE 2 0.0001
## TRUE 19998 0.9999
# Exercise 9.13
#a)
# Predict rides for each parameter set in the chain
set.seed(84735)
predict_90 <- bike_posterior_model_df %>%
mutate(mu = `(Intercept)` + humidity*90,
y_new = rnorm(20000, mean = mu, sd = sigma))
head(predict_90, 3)
## (Intercept) humidity sigma mu y_new
## 1 4055.334 -9.104262 1618.898 3235.951 4316.132
## 2 4340.841 -13.985910 1659.860 3082.109 2876.551
## 3 4246.110 -12.485416 1652.535 3122.423 4412.834
#b)
# Plot the posterior model of the typical ridership on 90% humidity days
d_mu <- ggplot(predict_90, aes(x = mu)) +
geom_density() +
xlim(-3000, 9000) +
ylim (0, 0.004)
# Plot the posterior predictive model of tomorrow's ridership
d_ynew <- ggplot(predict_90, aes(x = y_new)) +
geom_density() +
xlim(-3000, 9000) +
ylim (0, 0.004)
d_mu + d_ynew
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_density()`).
#c)
# 80% posterior prediction interval or the number of riders tomorrow
predict_90 %>%
summarize(lower_new = quantile(y_new, 0.1),
upper_new = quantile(y_new, 0.9))
## lower_new upper_new
## 1 1235.484 5258.443
#d)
# Simulate a set of predictions
set.seed(84735)
shortcut_prediction <-
posterior_predict(bike_posterior_model, newdata = data.frame(humidity = 90))
# Construct a 80% posterior credible interval
posterior_interval(shortcut_prediction, prob = 0.8)
## 10% 90%
## 1 1235.484 5258.443
# Plot the approximate predictive model
mcmc_dens(shortcut_prediction)+
xlab("predicted ridership on a 90% humidity day")
