Installations
## INSTALL CMDSTANR
install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
## INSTALL CMDSTAN
cmdstanr::install_cmdstan()
install.packages("bayesplot")
install.packages("readr")
install.packages("dplyr")
Setup
library(cmdstanr)
## This is cmdstanr version 0.7.1
## - CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
## - CmdStan path: /Users/arvindsharma/.cmdstan/cmdstan-2.34.1
## - CmdStan version: 2.34.1
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(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ 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(ggplot2)
Data
data <- read.csv("golf_data.csv")
stan_data <- list(N = nrow(data), distance_feet = data$distance_feet,
tries = data$tries, successes = data$successes)
## plot data without axes
(ggplot(data, aes(x = distance_feet, y = successes / tries))
+ geom_point(size = 3)
+ xlim(0, 20) + ylim(0, 1)
+ theme_bw()
+ theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.margin = margin(t = 0, b = 0),
text = element_text(size = 20))
)
## plot data with axes
(ggplot(data, aes(x = distance_feet, y = successes / tries))
+ geom_point(size = 3)
+ geom_segment(aes(x = distance_feet, xend = distance_feet, y = successes / tries - 0.5, yend = successes / tries + 0.5))
+ xlim(0, 20) + ylim(0, 1)
+ theme_bw()
+ theme(plot.margin = margin(t = 0, b = 0),
text = element_text(size = 20))
)
## Warning: Removed 19 rows containing missing values (`geom_segment()`).
Check that cmdstanr
works
model <- cmdstanr::cmdstan_model("golf_setup.stan")
fit_mcmc <- model$sample(data = stan_data)
## Running MCMC with 4 sequential chains...
##
## Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2 finished in 0.0 seconds.
## Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3 finished in 0.0 seconds.
## Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.6 seconds.
fit_optim <- model$optimize(data = stan_data)
## Initial log joint probability = -0.0345396
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 2 -3.85186e-34 0.236546 2.77556e-17 1 1 5
## Optimization terminated normally:
## Convergence detected: gradient norm is below tolerance
## Finished in 0.2 seconds.
fit_advi <- model$variational(data = stan_data)
## ------------------------------------------------------------
## EXPERIMENTAL ALGORITHM:
## This procedure has not been thoroughly tested and may be unstable
## or buggy. The interface is subject to change.
## ------------------------------------------------------------
## Gradient evaluation took 1e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.01 seconds.
## Adjust your expectations accordingly!
## Begin eta adaptation.
## Iteration: 1 / 250 [ 0%] (Adaptation)
## Iteration: 50 / 250 [ 20%] (Adaptation)
## Iteration: 100 / 250 [ 40%] (Adaptation)
## Iteration: 150 / 250 [ 60%] (Adaptation)
## Iteration: 200 / 250 [ 80%] (Adaptation)
## Success! Found best value [eta = 1] earlier than expected.
## Begin stochastic gradient ascent.
## iter ELBO delta_ELBO_mean delta_ELBO_med notes
## 100 0.095 1.000 1.000
## 200 0.061 0.784 1.000
## 300 -0.004 5.709 1.000
## 400 0.055 4.550 1.076
## 500 -0.082 3.975 1.076
## 600 -0.182 3.405 1.076
## 700 -0.083 3.089 1.076
## 800 0.036 3.120 1.193
## 900 0.025 2.820 1.076
## 1000 -0.010 2.892 1.193
## 1100 -0.112 2.883 1.193 MAY BE DIVERGING... INSPECT ELBO
## 1200 0.111 3.027 1.676 MAY BE DIVERGING... INSPECT ELBO
## 1300 -0.040 1.850 1.676 MAY BE DIVERGING... INSPECT ELBO
## 1400 0.020 2.038 2.009 MAY BE DIVERGING... INSPECT ELBO
## 1500 0.033 1.908 2.009 MAY BE DIVERGING... INSPECT ELBO
## 1600 -0.004 2.709 2.955 MAY BE DIVERGING... INSPECT ELBO
## 1700 -0.035 2.678 2.955 MAY BE DIVERGING... INSPECT ELBO
## 1800 0.052 2.512 2.009 MAY BE DIVERGING... INSPECT ELBO
## 1900 0.015 2.722 2.520 MAY BE DIVERGING... INSPECT ELBO
## 2000 0.096 2.452 2.009 MAY BE DIVERGING... INSPECT ELBO
## 2100 0.128 2.386 2.009 MAY BE DIVERGING... INSPECT ELBO
## 2200 -0.041 2.599 2.520 MAY BE DIVERGING... INSPECT ELBO
## 2300 0.019 2.534 2.520 MAY BE DIVERGING... INSPECT ELBO
## 2400 -0.040 2.386 1.677 MAY BE DIVERGING... INSPECT ELBO
## 2500 0.027 2.600 2.514 MAY BE DIVERGING... INSPECT ELBO
## 2600 0.044 1.783 1.677 MAY BE DIVERGING... INSPECT ELBO
## 2700 -0.036 1.919 2.231 MAY BE DIVERGING... INSPECT ELBO
## 2800 -0.019 1.835 2.231 MAY BE DIVERGING... INSPECT ELBO
## 2900 0.070 1.710 1.477 MAY BE DIVERGING... INSPECT ELBO
## 3000 -0.021 2.053 2.231 MAY BE DIVERGING... INSPECT ELBO
## 3100 0.014 2.281 2.514 MAY BE DIVERGING... INSPECT ELBO
## 3200 -0.094 1.983 2.231 MAY BE DIVERGING... INSPECT ELBO
## 3300 0.147 1.833 1.640 MAY BE DIVERGING... INSPECT ELBO
## 3400 -0.004 5.887 2.231 MAY BE DIVERGING... INSPECT ELBO
## 3500 0.047 5.744 1.640 MAY BE DIVERGING... INSPECT ELBO
## 3600 0.003 7.349 2.231 MAY BE DIVERGING... INSPECT ELBO
## 3700 0.023 7.214 1.640 MAY BE DIVERGING... INSPECT ELBO
## 3800 0.062 7.193 1.640 MAY BE DIVERGING... INSPECT ELBO
## 3900 -0.055 7.278 2.122 MAY BE DIVERGING... INSPECT ELBO
## 4000 0.001 11.447 2.122 MAY BE DIVERGING... INSPECT ELBO
## 4100 -0.049 11.297 1.640 MAY BE DIVERGING... INSPECT ELBO
## 4200 -0.054 11.191 1.640 MAY BE DIVERGING... INSPECT ELBO
## 4300 -0.075 11.056 1.076 MAY BE DIVERGING... INSPECT ELBO
## 4400 0.054 7.094 1.076 MAY BE DIVERGING... INSPECT ELBO
## 4500 -0.012 7.544 2.122 MAY BE DIVERGING... INSPECT ELBO
## 4600 -0.045 5.973 1.025 MAY BE DIVERGING... INSPECT ELBO
## 4700 -0.024 5.973 1.025 MAY BE DIVERGING... INSPECT ELBO
## 4800 -0.085 5.981 1.025 MAY BE DIVERGING... INSPECT ELBO
## 4900 0.018 6.341 1.025 MAY BE DIVERGING... INSPECT ELBO
## 5000 -0.065 1.872 1.025 MAY BE DIVERGING... INSPECT ELBO
## 5100 -0.028 1.906 1.276 MAY BE DIVERGING... INSPECT ELBO
## 5200 0.047 2.055 1.358 MAY BE DIVERGING... INSPECT ELBO
## 5300 0.031 2.078 1.358 MAY BE DIVERGING... INSPECT ELBO
## 5400 0.015 1.938 1.276 MAY BE DIVERGING... INSPECT ELBO
## 5500 -0.080 1.500 1.193 MAY BE DIVERGING... INSPECT ELBO
## 5600 -0.048 1.493 1.193 MAY BE DIVERGING... INSPECT ELBO
## 5700 -0.038 1.432 1.193 MAY BE DIVERGING... INSPECT ELBO
## 5800 0.001 4.444 1.276 MAY BE DIVERGING... INSPECT ELBO
## 5900 0.133 3.970 1.193 MAY BE DIVERGING... INSPECT ELBO
## 6000 0.116 3.857 1.007 MAY BE DIVERGING... INSPECT ELBO
## 6100 -0.033 4.177 1.007 MAY BE DIVERGING... INSPECT ELBO
## 6200 -0.001 10.006 1.007 MAY BE DIVERGING... INSPECT ELBO
## 6300 -0.028 10.053 1.007 MAY BE DIVERGING... INSPECT ELBO
## 6400 0.052 10.106 1.193 MAY BE DIVERGING... INSPECT ELBO
## 6500 0.035 10.034 0.990 MAY BE DIVERGING... INSPECT ELBO
## 6600 0.076 10.019 0.990 MAY BE DIVERGING... INSPECT ELBO
## 6700 0.109 10.023 0.990 MAY BE DIVERGING... INSPECT ELBO
## 6800 -0.012 7.978 0.990 MAY BE DIVERGING... INSPECT ELBO
## 6900 0.128 7.988 1.090 MAY BE DIVERGING... INSPECT ELBO
## 7000 0.080 8.035 1.090 MAY BE DIVERGING... INSPECT ELBO
## 7100 -0.001 19.927 1.090 MAY BE DIVERGING... INSPECT ELBO
## 7200 -0.010 14.032 0.981 MAY BE DIVERGING... INSPECT ELBO
## 7300 -0.012 13.954 0.933 MAY BE DIVERGING... INSPECT ELBO
## 7400 0.101 13.912 0.933 MAY BE DIVERGING... INSPECT ELBO
## 7500 -0.093 14.074 1.090 MAY BE DIVERGING... INSPECT ELBO
## 7600 -0.087 14.028 1.090 MAY BE DIVERGING... INSPECT ELBO
## 7700 0.006 15.558 1.121 MAY BE DIVERGING... INSPECT ELBO
## 7800 -0.144 14.624 1.090 MAY BE DIVERGING... INSPECT ELBO
## 7900 -0.083 14.589 1.041 MAY BE DIVERGING... INSPECT ELBO
## 8000 -0.053 14.585 1.041 MAY BE DIVERGING... INSPECT ELBO
## 8100 -0.211 2.311 0.933 MAY BE DIVERGING... INSPECT ELBO
## 8200 -0.051 2.534 1.041 MAY BE DIVERGING... INSPECT ELBO
## 8300 -0.054 2.520 1.041 MAY BE DIVERGING... INSPECT ELBO
## 8400 0.027 2.709 1.041 MAY BE DIVERGING... INSPECT ELBO
## 8500 -0.004 3.268 1.041 MAY BE DIVERGING... INSPECT ELBO
## 8600 0.090 3.365 1.044 MAY BE DIVERGING... INSPECT ELBO
## 8700 -0.099 1.995 1.044 MAY BE DIVERGING... INSPECT ELBO
## 8800 0.176 2.047 1.564 MAY BE DIVERGING... INSPECT ELBO
## 8900 -0.040 2.513 1.909 MAY BE DIVERGING... INSPECT ELBO
## 9000 -0.038 2.464 1.909 MAY BE DIVERGING... INSPECT ELBO
## 9100 0.062 2.550 1.909 MAY BE DIVERGING... INSPECT ELBO
## 9200 0.012 2.646 1.909 MAY BE DIVERGING... INSPECT ELBO
## 9300 -0.157 2.748 1.909 MAY BE DIVERGING... INSPECT ELBO
## 9400 -0.108 2.492 1.607 MAY BE DIVERGING... INSPECT ELBO
## 9500 -0.001 10.351 1.607 MAY BE DIVERGING... INSPECT ELBO
## 9600 -0.084 10.345 1.607 MAY BE DIVERGING... INSPECT ELBO
## 9700 -0.233 10.218 1.564 MAY BE DIVERGING... INSPECT ELBO
## 9800 0.071 10.488 1.607 MAY BE DIVERGING... INSPECT ELBO
## 9900 -0.177 10.088 1.403 MAY BE DIVERGING... INSPECT ELBO
## 10000 -0.112 10.140 1.403 MAY BE DIVERGING... INSPECT ELBO
## Informational Message: The maximum number of iterations is reached! The algorithm may not have converged.
## This variational approximation is not guaranteed to be meaningful.
## Drawing a sample of size 1000 from the approximate posterior...
## COMPLETED.
## Finished in 0.2 seconds.
Fit linear model
model <- cmdstanr::cmdstan_model("golf_linear.stan")
fit_mcmc <- model$sample(data = stan_data)
## Running MCMC with 4 sequential chains...
##
## Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2 finished in 0.1 seconds.
## Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/var/folders/0h/kk_vktsd17q_p7867shj39zw0000gr/T/Rtmp8iFPcd/model-c5973dfa93a9.stan', line 16, column 2 to column 43)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.5 seconds.
fit_optim <- model$optimize(data = stan_data)
## Initial log joint probability = -6294.25
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 40 36.9545 7.19476e-05 0.0296863 0.6326 0.06326 67
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
fit_advi <- model$variational(data = stan_data)
## ------------------------------------------------------------
## EXPERIMENTAL ALGORITHM:
## This procedure has not been thoroughly tested and may be unstable
## or buggy. The interface is subject to change.
## ------------------------------------------------------------
## Gradient evaluation took 6e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Adjust your expectations accordingly!
## Begin eta adaptation.
## Iteration: 1 / 250 [ 0%] (Adaptation)
## Iteration: 50 / 250 [ 20%] (Adaptation)
## Iteration: 100 / 250 [ 40%] (Adaptation)
## Iteration: 150 / 250 [ 60%] (Adaptation)
## Iteration: 200 / 250 [ 80%] (Adaptation)
## Success! Found best value [eta = 1] earlier than expected.
## Begin stochastic gradient ascent.
## iter ELBO delta_ELBO_mean delta_ELBO_med notes
## 100 -68.846 1.000 1.000
## 200 -41.593 0.828 1.000
## 300 -15.776 1.097 1.000
## 400 -23.240 0.903 1.000
## 500 -18.341 0.776 0.655
## 600 -5.850 1.003 1.000
## 700 -12.690 0.936 0.655
## 800 -10.805 0.841 0.655
## 900 -22.455 0.805 0.539
## 1000 -2.340 1.584 0.655
## 1100 -17.364 1.571 0.655 MAY BE DIVERGING... INSPECT ELBO
## 1200 1.252 2.992 0.865 MAY BE DIVERGING... INSPECT ELBO
## 1300 0.884 2.870 0.539 MAY BE DIVERGING... INSPECT ELBO
## 1400 -1.031 3.024 0.865 MAY BE DIVERGING... INSPECT ELBO
## 1500 0.083 4.344 1.857 MAY BE DIVERGING... INSPECT ELBO
## 1600 -21.765 4.231 1.004 MAY BE DIVERGING... INSPECT ELBO
## 1700 0.930 6.617 1.857 MAY BE DIVERGING... INSPECT ELBO
## 1800 1.463 6.636 1.857 MAY BE DIVERGING... INSPECT ELBO
## 1900 -2.357 6.746 1.857 MAY BE DIVERGING... INSPECT ELBO
## 2000 1.297 6.168 1.857 MAY BE DIVERGING... INSPECT ELBO
## 2100 -2.254 6.239 1.857 MAY BE DIVERGING... INSPECT ELBO
## 2200 -7.059 4.821 1.621 MAY BE DIVERGING... INSPECT ELBO
## 2300 3.736 5.068 1.857 MAY BE DIVERGING... INSPECT ELBO
## 2400 3.881 4.886 1.621 MAY BE DIVERGING... INSPECT ELBO
## 2500 3.919 3.540 1.575 MAY BE DIVERGING... INSPECT ELBO
## 2600 2.084 3.527 1.575 MAY BE DIVERGING... INSPECT ELBO
## 2700 2.468 1.103 0.881 MAY BE DIVERGING... INSPECT ELBO
## 2800 2.339 1.072 0.881 MAY BE DIVERGING... INSPECT ELBO
## 2900 -4.360 1.064 0.881 MAY BE DIVERGING... INSPECT ELBO
## 3000 2.644 1.047 0.881 MAY BE DIVERGING... INSPECT ELBO
## 3100 4.381 0.929 0.681 MAY BE DIVERGING... INSPECT ELBO
## 3200 4.343 0.862 0.397 MAY BE DIVERGING... INSPECT ELBO
## 3300 -5.446 0.753 0.397 MAY BE DIVERGING... INSPECT ELBO
## 3400 -3.162 0.821 0.723 MAY BE DIVERGING... INSPECT ELBO
## 3500 3.032 1.025 0.881 MAY BE DIVERGING... INSPECT ELBO
## 3600 0.854 1.192 1.536 MAY BE DIVERGING... INSPECT ELBO
## 3700 -4.058 1.297 1.536 MAY BE DIVERGING... INSPECT ELBO
## 3800 4.128 1.490 1.797 MAY BE DIVERGING... INSPECT ELBO
## 3900 4.975 1.353 1.797 MAY BE DIVERGING... INSPECT ELBO
## 4000 -16.945 1.218 1.294 MAY BE DIVERGING... INSPECT ELBO
## 4100 4.484 1.656 1.797 MAY BE DIVERGING... INSPECT ELBO
## 4200 -10.772 1.797 1.797 MAY BE DIVERGING... INSPECT ELBO
## 4300 0.890 2.927 1.983 MAY BE DIVERGING... INSPECT ELBO
## 4400 4.767 2.936 1.983 MAY BE DIVERGING... INSPECT ELBO
## 4500 0.723 3.292 1.983 MAY BE DIVERGING... INSPECT ELBO
## 4600 3.653 3.117 1.416 MAY BE DIVERGING... INSPECT ELBO
## 4700 0.650 3.458 1.983 MAY BE DIVERGING... INSPECT ELBO
## 4800 -8.435 3.367 1.416 MAY BE DIVERGING... INSPECT ELBO
## 4900 4.225 3.650 2.996 MAY BE DIVERGING... INSPECT ELBO
## 5000 2.970 3.563 2.996 MAY BE DIVERGING... INSPECT ELBO
## 5100 3.648 3.103 1.416 MAY BE DIVERGING... INSPECT ELBO
## 5200 2.041 3.040 1.077 MAY BE DIVERGING... INSPECT ELBO
## 5300 -6.664 1.861 1.077 MAY BE DIVERGING... INSPECT ELBO
## 5400 3.897 2.050 1.306 MAY BE DIVERGING... INSPECT ELBO
## 5500 -17.473 1.613 1.223 MAY BE DIVERGING... INSPECT ELBO
## 5600 6.362 1.908 1.306 MAY BE DIVERGING... INSPECT ELBO
## 5700 5.035 1.472 1.223 MAY BE DIVERGING... INSPECT ELBO
## 5800 3.990 1.390 1.223 MAY BE DIVERGING... INSPECT ELBO
## 5900 5.616 1.120 0.788 MAY BE DIVERGING... INSPECT ELBO
## 6000 5.076 1.088 0.788 MAY BE DIVERGING... INSPECT ELBO
## 6100 5.530 1.078 0.788 MAY BE DIVERGING... INSPECT ELBO
## 6200 4.324 1.027 0.289 MAY BE DIVERGING... INSPECT ELBO
## 6300 4.603 0.902 0.279 MAY BE DIVERGING... INSPECT ELBO
## 6400 6.031 0.655 0.264 MAY BE DIVERGING... INSPECT ELBO
## 6500 5.712 0.538 0.262 MAY BE DIVERGING... INSPECT ELBO
## 6600 5.456 0.168 0.237
## 6700 5.500 0.143 0.107
## 6800 6.545 0.132 0.107
## 6900 1.213 0.543 0.107 MAY BE DIVERGING... INSPECT ELBO
## 7000 3.495 0.598 0.160 MAY BE DIVERGING... INSPECT ELBO
## 7100 2.679 0.620 0.237 MAY BE DIVERGING... INSPECT ELBO
## 7200 2.743 0.594 0.160 MAY BE DIVERGING... INSPECT ELBO
## 7300 4.992 0.633 0.237 MAY BE DIVERGING... INSPECT ELBO
## 7400 5.943 0.626 0.160 MAY BE DIVERGING... INSPECT ELBO
## 7500 6.226 0.625 0.160 MAY BE DIVERGING... INSPECT ELBO
## 7600 3.658 0.690 0.305 MAY BE DIVERGING... INSPECT ELBO
## 7700 6.410 0.732 0.429 MAY BE DIVERGING... INSPECT ELBO
## 7800 4.030 0.776 0.450 MAY BE DIVERGING... INSPECT ELBO
## 7900 5.679 0.365 0.429
## 8000 6.178 0.308 0.305
## 8100 -5.382 0.492 0.429
## 8200 2.027 0.855 0.450 MAY BE DIVERGING... INSPECT ELBO
## 8300 -17.463 0.922 0.591 MAY BE DIVERGING... INSPECT ELBO
## 8400 3.856 1.459 0.702 MAY BE DIVERGING... INSPECT ELBO
## 8500 -7.155 1.608 1.116 MAY BE DIVERGING... INSPECT ELBO
## 8600 6.542 1.747 1.539 MAY BE DIVERGING... INSPECT ELBO
## 8700 6.391 1.707 1.539 MAY BE DIVERGING... INSPECT ELBO
## 8800 5.342 1.667 1.539 MAY BE DIVERGING... INSPECT ELBO
## 8900 5.277 1.639 1.539 MAY BE DIVERGING... INSPECT ELBO
## 9000 6.287 1.647 1.539 MAY BE DIVERGING... INSPECT ELBO
## 9100 5.926 1.439 1.116 MAY BE DIVERGING... INSPECT ELBO
## 9200 -0.266 3.398 1.116 MAY BE DIVERGING... INSPECT ELBO
## 9300 2.899 3.395 1.092 MAY BE DIVERGING... INSPECT ELBO
## 9400 2.716 2.849 0.196 MAY BE DIVERGING... INSPECT ELBO
## 9500 5.058 2.742 0.196 MAY BE DIVERGING... INSPECT ELBO
## 9600 5.065 2.532 0.161 MAY BE DIVERGING... INSPECT ELBO
## 9700 6.640 2.554 0.196 MAY BE DIVERGING... INSPECT ELBO
## 9800 6.904 2.538 0.161 MAY BE DIVERGING... INSPECT ELBO
## 9900 -5.712 2.758 0.237 MAY BE DIVERGING... INSPECT ELBO
## 10000 5.302 2.949 0.463 MAY BE DIVERGING... INSPECT ELBO
## Informational Message: The maximum number of iterations is reached! The algorithm may not have converged.
## This variational approximation is not guaranteed to be meaningful.
## Drawing a sample of size 1000 from the approximate posterior...
## COMPLETED.
## Finished in 0.2 seconds.
(ggplot(data, aes(x = distance_feet, y = successes / tries))
+ geom_point(size = 3)
+ stat_function(fun = function(x) fit_optim$mle()['a'] + fit_optim$mle()['b'] * x, linewidth = 2)
+ xlim(0, 20) + ylim(0, 1)
+ theme_bw()
+ theme(plot.margin = margin(t = 0, b = 0),
text = element_text(size = 20))
#+ theme(axis.title.x = element_blank(),
#axis.title.y = element_blank())
)
bayesplot::mcmc_hist(fit_mcmc$draws(c('a', 'b', 'sigma')))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
fit_optim$mle(c('a', 'b', 'sigma'))
## a b sigma
## 0.8360840 -0.0408875 0.0867252
Fit logistic
model
![]()
model <- cmdstanr::cmdstan_model("golf_logistic.stan")
fit_mcmc <- model$sample(data = stan_data)
## Running MCMC with 4 sequential chains...
##
## Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2 finished in 0.0 seconds.
## Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3 finished in 0.0 seconds.
## Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.5 seconds.
fit_optim <- model$optimize(data = stan_data)
## Initial log joint probability = -55430.4
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 9 -3020.15 3.23616e-05 0.0772476 0.1508 1 17
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
fit_advi <- model$variational(data = stan_data)
## ------------------------------------------------------------
## EXPERIMENTAL ALGORITHM:
## This procedure has not been thoroughly tested and may be unstable
## or buggy. The interface is subject to change.
## ------------------------------------------------------------
## Gradient evaluation took 6e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Adjust your expectations accordingly!
## Begin eta adaptation.
## Iteration: 1 / 250 [ 0%] (Adaptation)
## Iteration: 50 / 250 [ 20%] (Adaptation)
## Iteration: 100 / 250 [ 40%] (Adaptation)
## Iteration: 150 / 250 [ 60%] (Adaptation)
## Iteration: 200 / 250 [ 80%] (Adaptation)
## Success! Found best value [eta = 1] earlier than expected.
## Begin stochastic gradient ascent.
## iter ELBO delta_ELBO_mean delta_ELBO_med notes
## 100 -295.589 1.000 1.000
## 200 -247.609 0.597 1.000
## 300 -210.296 0.457 0.194
## 400 -194.777 0.363 0.194
## 500 -209.129 0.304 0.177
## 600 -194.470 0.266 0.177
## 700 -206.721 0.236 0.080
## 800 -188.372 0.219 0.097
## 900 -191.304 0.196 0.080
## 1000 -198.204 0.180 0.080
## 1100 -187.972 0.086 0.075
## 1200 -189.699 0.067 0.069
## 1300 -191.455 0.050 0.059
## 1400 -203.690 0.048 0.059
## 1500 -188.661 0.049 0.059
## 1600 -222.633 0.057 0.059
## 1700 -189.531 0.069 0.060
## 1800 -188.416 0.060 0.054
## 1900 -188.408 0.058 0.054
## 2000 -194.233 0.058 0.054
## 2100 -190.617 0.054 0.030
## 2200 -188.691 0.054 0.030
## 2300 -188.201 0.053 0.030
## 2400 -199.053 0.053 0.030
## 2500 -209.863 0.050 0.030
## 2600 -200.509 0.040 0.030
## 2700 -196.438 0.024 0.021
## 2800 -188.093 0.028 0.030
## 2900 -208.323 0.038 0.044
## 3000 -201.550 0.038 0.044
## 3100 -188.338 0.043 0.047
## 3200 -202.754 0.049 0.052
## 3300 -191.540 0.055 0.055
## 3400 -188.049 0.051 0.052
## 3500 -205.658 0.055 0.059
## 3600 -190.818 0.058 0.070
## 3700 -189.054 0.057 0.070
## 3800 -193.479 0.054 0.070
## 3900 -188.588 0.047 0.059
## 4000 -195.256 0.047 0.059
## 4100 -190.392 0.043 0.034
## 4200 -189.660 0.036 0.026
## 4300 -188.380 0.031 0.026
## 4400 -195.633 0.033 0.026
## 4500 -192.527 0.026 0.026
## 4600 -190.556 0.019 0.023
## 4700 -203.891 0.025 0.026
## 4800 -195.729 0.027 0.026
## 4900 -192.966 0.026 0.026
## 5000 -188.345 0.025 0.025
## 5100 -190.367 0.023 0.016
## 5200 -192.293 0.024 0.016
## 5300 -204.072 0.029 0.025
## 5400 -201.161 0.027 0.016
## 5500 -188.340 0.032 0.025
## 5600 -189.432 0.031 0.025
## 5700 -193.915 0.027 0.023
## 5800 -188.760 0.026 0.023
## 5900 -188.964 0.024 0.023
## 6000 -189.290 0.022 0.014
## 6100 -188.275 0.021 0.014
## 6200 -188.788 0.021 0.014
## 6300 -188.317 0.015 0.006 MEDIAN ELBO CONVERGED
## Drawing a sample of size 1000 from the approximate posterior...
## COMPLETED.
## Finished in 0.1 seconds.
(ggplot(data, aes(x = distance_feet, y = successes / tries))
+ geom_point(size = 3)
+ stat_function(fun = function(x) boot::inv.logit(fit_optim$mle()['a'] + fit_optim$mle()['b'] * x), linewidth = 2)
+ xlim(0, 20) + ylim(0, 1)
+ theme_bw()
+ theme(plot.margin = margin(t = 0, b = 0),
text = element_text(size = 20))
#+ theme(axis.title.x = element_blank(),
# axis.title.y = element_blank())
)
Fit angle model
model <- cmdstanr::cmdstan_model("golf_angle.stan")
fit_mcmc <- model$sample(data = stan_data)
## Running MCMC with 4 sequential chains...
##
## Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2 finished in 0.0 seconds.
## Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3 finished in 0.0 seconds.
## Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.6 seconds.
fit_optim <- model$optimize(data = stan_data)
## Initial log joint probability = -7945.55
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 4 -2922.64 0.000121656 0.0275921 0.9523 0.9523 7
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
fit_advi <- model$variational(data = stan_data)
## ------------------------------------------------------------
## EXPERIMENTAL ALGORITHM:
## This procedure has not been thoroughly tested and may be unstable
## or buggy. The interface is subject to change.
## ------------------------------------------------------------
## Gradient evaluation took 6e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Adjust your expectations accordingly!
## Begin eta adaptation.
## Iteration: 1 / 250 [ 0%] (Adaptation)
## Iteration: 50 / 250 [ 20%] (Adaptation)
## Iteration: 100 / 250 [ 40%] (Adaptation)
## Iteration: 150 / 250 [ 60%] (Adaptation)
## Iteration: 200 / 250 [ 80%] (Adaptation)
## Success! Found best value [eta = 1] earlier than expected.
## Begin stochastic gradient ascent.
## iter ELBO delta_ELBO_mean delta_ELBO_med notes
## 100 -93.082 1.000 1.000
## 200 -91.471 0.509 1.000
## 300 -90.682 0.342 0.018
## 400 -94.401 0.266 0.039
## 500 -91.511 0.219 0.032
## 600 -91.193 0.183 0.032
## 700 -90.618 0.158 0.018
## 800 -93.417 0.142 0.030
## 900 -90.435 0.130 0.030
## 1000 -93.329 0.120 0.031
## 1100 -90.529 0.023 0.031
## 1200 -90.622 0.022 0.031
## 1300 -90.862 0.021 0.031
## 1400 -91.323 0.018 0.030
## 1500 -93.451 0.017 0.023
## 1600 -91.316 0.019 0.023
## 1700 -90.373 0.019 0.023
## 1800 -90.460 0.016 0.023
## 1900 -90.362 0.013 0.010
## 2000 -90.400 0.010 0.005 MEAN ELBO CONVERGED MEDIAN ELBO CONVERGED
## Drawing a sample of size 1000 from the approximate posterior...
## COMPLETED.
## Finished in 0.1 seconds.
r = 1.68 / 12 / 2
R = 4.25 / 12 / 2
(ggplot(data, aes(x = distance_feet, y = successes / tries))
+ geom_point(size = 3)
+ stat_function(fun = function(x)
2 * pnorm(asin((R - r) / x) / fit_optim$mle()['sigma']) - 1, linewidth = 2)
+ xlim(0, 20) + ylim(0, 1)
+ theme_bw()
+ theme(plot.margin = margin(t = 0, b = 0),
text = element_text(size = 20))
#+ theme(axis.title.x = element_blank(),
# axis.title.y = element_blank())
)
## Warning in asin((R - r)/x): NaNs produced
## Warning: Removed 1 row containing missing values (`geom_function()`).
bayesplot::mcmc_hist(fit_mcmc$draws('sigma_degrees'))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
