Andrew Gelman, Jennifer Hill, Aki Vehtari 2021-04-20
Tidyverse version by Bill Behrman.
Demonstration of using Stan for optimization. See Appendix A in Regression and Other Stories.
# Packages
library(tidyverse)
library(rstan)
# Parameters
# Common code
file_common <- here::here("_common.R")
#===============================================================================
# Run common code
source(file_common)
Plot of y = 15 + 10x − 2x2.
v <-
tibble(
x = seq_range(c(-2, 5)),
y = 15 + 10 * x - 2 * x^2
)
v %>%
ggplot(aes(x, y)) +
geom_line() +
geom_point(data = tibble(x = 2.5, y = 27.5)) +
labs(
title = expression(paste("Plot of ", y == 15 + 10 * x - 2 * x^2)),
subtitle = "Maximum is at x = 2.5, y = 27.5"
)
Stan model for function 15 + 10x − 2x2.
model_code =
"
parameters {
real x;
}
model {
target += 15 + 10 * x - 2 * x^2;
}
"
Compile the Stan function and optimize it.
set.seed(245)
fit <-
stan_model(model_code = model_code) %>%
optimizing()
fit
#> $par
#> x
#> 2.5
#>
#> $value
#> [1] 27.5
#>
#> $return_code
#> [1] 0
#>
#> $theta_tilde
#> x
#> [1,] 2.5
The output labeled par
is the value of x
at which the function is
optimized, and the output labeled value
is the target function at the
optimum. A return_code
of 0 corresponds to the optimizer running with
no problem, and you can ignore theta_tilde
.