Last active
April 23, 2018 00:42
-
-
Save alexhallam/5fa6a3475b5347e8f5d56305c929b48d to your computer and use it in GitHub Desktop.
Debug some stan code
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# load libraries | |
library(purrr) | |
library(tidyverse) | |
# positive_round function | |
positive_round <- function(...) round(pmax(..., 0), 0) | |
# adstock function | |
adstock<-function(x,rate=0){ | |
return(as.numeric(stats::filter(x=x,filter=rate,method="recursive"))) | |
} | |
my_adstock <- function(x, rate){ | |
y <- NULL | |
my_vector <- vector("numeric") | |
for(i in seq(x)) { | |
# Adstock: y[i] = x[i] + f[1]*y[i-1] + … + f[p]*y[i-p] | |
y[1] <- x[i] | |
y[i + 1] <- x[i + 1] + rate * y[i + 1 - 1] | |
# store as vector | |
my_vector[i] <- y[i] | |
} | |
# return results of vector | |
my_vector | |
} | |
# base and ad constants | |
base <- 100 | |
n_weeks <- 104 | |
avg <- 20 | |
sdev <- 10 | |
# adstock rates | |
ad1_rate <- .7 | |
ad2_rate <- .4 | |
ad3_rate <- .5 | |
set.seed(42) | |
# 3 ads with mean | |
fake_df <- rep(avg, 3) %>% | |
# normal distribution with 2 years of observations | |
map(~ rnorm(n = n_weeks, mean = .x, sd = sdev)) %>% | |
# make sure there are no 0s and make values integers | |
map(~ positive_round(0,.x)) %>% | |
# convert lists to data frame | |
map_dfc(~as_tibble(.x)) %>% | |
# sensible names | |
setNames(c("ad1", "ad2", "ad3")) %>% | |
# sales = sum of adstock and some noise | |
mutate(sales = base + | |
adstock(ad1, ad1_rate) + | |
adstock(ad2, ad2_rate) + | |
adstock(ad3, ad3_rate) + | |
rnorm(n(), sd = 5)) %>% | |
# round sales to whole number | |
mutate(sales = round(sales, 0)) | |
library(brms) | |
prior1 <- prior(normal(.5, .5), nlpar = "ad1rate") + | |
prior(normal(.5, .5), nlpar = "ad2rate") + | |
prior(normal(.5, .5), nlpar = "ad3rate") + | |
prior(normal(1, .5), nlpar = "b1") + | |
prior(normal(1, .5), nlpar = "b2") + | |
prior(normal(1, .5), nlpar = "b3") + | |
prior(normal(100, 50), nlpar = "basesales") | |
my_stan_adstock <- " | |
vector my_stan_adstock(vector x, real rate, int N){ | |
vector[N] y; | |
vector[N] my_vector; | |
for(i in 1:N){ | |
y[0] = x[i]; | |
y[i] = x[i] + rate * y[i - 1]; | |
my_vector[i] = y[i]; | |
} | |
return my_vector; | |
} | |
" | |
fit1 <- brm(bf(sales ~ basesales + | |
b1 * my_stan_adstock(ad1, ad1rate, 104) + | |
b2 * my_stan_adstock(ad2, ad2rate, 104) + | |
b3 * my_stan_adstock(ad3, ad3rate, 104), | |
basesales + b1 + b2 + b3 + ad1rate + ad2rate + ad3rate ~ 1, | |
nl = TRUE), | |
stan_funs = my_stan_adstock, | |
data = fake_df, | |
prior = prior1) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment