devindg / pybuc

Bayesian Structural Time Series / Unobserved Components
BSD 3-Clause "New" or "Revised" License
29 stars 6 forks source link

Adding Priors Into Model #28

Closed tim-mcwilliams closed 1 year ago

tim-mcwilliams commented 1 year ago

Hi @devindg !!!

First off, such an awesome package that is easy to use and understand.

I was wondering how would one go about adding priors into a multivariate bsts model. Would you use the state_var_scale_prior variable in the modeling set up? Could you provide an example, please?!

devindg commented 1 year ago

Hi @tim-mcwilliams . Happy that you like the package! I'll actually (hopefully) be releasing a new version soon that will allow you to see the posterior distributions and Monte Carlo sampling traces for each of the parameters. But to answer your question: what do you mean by "multivariate"? What parameters do you want to put priors on?

In general, priors will be passed to the sample() method. Here is an example:

from pybuc.buc import BayesianUnobservedComponents as BUC

mod = BUC(response=y_train,
                        level=True,
                        trend=True, damped_trend=True,
                        lag_seasonal=(4, 12),
                        seed=123)
mod.sample(num_samp=1000, 
                        level_var_shape_prior=1, 
                        level_var_scale_prior=0.01, 
                        trend_var_shape_prior=1, 
                        trend_var_scale_prior=0.01, 
                        lag_season_var_shape_prior=(1, 2), 
                        lag_season_var_scale_prior=(0.01, 0.1),
                        damped_trend_coeff_mean_prior=[0.95], 
                        damped_trend_coeff_prec_prior=[4])
tim-mcwilliams commented 1 year ago

Hi @devindg!

Thanks for the quick reply here! Looking forward to those updates!

When I say "multivariate" I mean that the model has predictors in it. For example,

bayes_uc = buc.BayesianUnobservedComponents(
    response=y_train,
    predictors=x_train,
    level=True,
    stochastic_level=True,
    trend=True,
    stochastic_trend=True,
    trig_seasonal=((12, 0),),
    stochastic_trig_seasonal=(True,),
    seed=123
)

It may be my lack of fully understanding the Bayesian way (I'm relatively new to it), but is there way to add priors to the predictors or are those priors only for the response variable?

devindg commented 1 year ago

Good morning, @tim-mcwilliams. Yes, there's a way to add a prior for the predictor coefficients. For example, let's suppose you have two predictors, $x_1$ and $x_2$, and your mean and covariance priors are, respectively, [2, 3] and [[0.5, 0], [0, 1]]. That is, your prior for the $x_1$ variable is $N(2, 0.5)$, and your prior for the $x_2$ variable is $N(3, 1)$. Since your covariance matrix is a diagonal matrix (zeros everywhere on the off-diagonal), the corresponding precision matrix, which is just the inverse of your covariance matrix, is [[1 / 0.5, 0], [0, 1]]. The precision matrix, as opposed to the covariance matrix, is what will be passed to sample(). The code is as follows:

bayes_uc = buc.BayesianUnobservedComponents(
    response=y_train,
    predictors=x_train,
    level=True,
    stochastic_level=True,
    trend=True,
    stochastic_trend=True,
    trig_seasonal=((12, 0),),
    stochastic_trig_seasonal=(True,),
    seed=123
)
post = bayes_uc.sample(num_samp=1000, reg_coeff_mean_prior=[2, 3], reg_coeff_prec_prior=[[2, 0], [0, 1]])

Note that you aren't constrained to lists for the mean and precision priors; you can also use tuples or Numpy arrays. Also, in general, if you prefer to think in terms of covariance instead of precision, you'll have to invert your covariance matrix first to get the corresponding precision matrix. Alternatively, if you don't have or want to provide a specific precision matrix, you can access the zellner_prior_obs argument, which will dictate how much weight you give your mean prior. While pybuc does not exactly follow the implementation of Zellner's g-prior (because it's not strictly regression), the principle is the same and can be read about here. pybuc follows the same approach as bsts for the implementation of Zellner's g-prior. You can read about it more here on pg. 8.

Finally, just as an FYI, you don't need to explicitly pass True for the stochastic_<component> arguments. These are True by default. Hope this helps.

tim-mcwilliams commented 1 year ago

Good morning, @devindg,

Really appreciate the detailed explanation here as well as providing some extra documents to digest. Going to work through my example and will let you know if I have any further questions.

Roger that on the stochastic_<component> tip. Just had those there for my learning purpose.

tim-mcwilliams commented 1 year ago

Closing!