<jupyter_start><jupyter_text>WBIC This notebook gives a tutorial on how to use Watanabe-Bayesian information criterion (WBIC) for feature selection. The WBIC is an information criterion. Similarly to other criteria (AIC, BIC, DIC) the WBIC endeavors to find the most parsimonious model, i.e., the model that balances fit with complexity. In other words a model (or set of features) that optimizes WBIC should neither over nor under fit the available data. In this tutorial a data set is simulated using the damped linear trend (DLT) model. This data set is then used to fit DLT models with varying number of features as well as a global local trent model (GLT), and a Error-Trend-Seasonal (ETS) model. The WBIC criteria is then show to find the true model.<jupyter_code>from datetime import timedelta

import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
%matplotlib inline

import orbit
from orbit.models import DLT,ETS, KTRLite, LGT
from orbit.utils.simulation import make_trend, make_regression
print(orbit.__version__)<jupyter_output>1.1.1dev<jupyter_text>Data Simulation This block of code creats random data set (365 observations with 10 features) assuming a DLT model. Of the 10 features 5 are effective regressors; i.e., they are used in the true model to create the data.As an exerise left to the user once you have run the code once try changing the NUM_OF_EFFECTIVE_REGRESSORS (line 2), the SERIES_LEN (line 3), and the SEED (line 4) to see how it effects the results.<jupyter_code>NUM_OF_REGRESSORS = 10
NUM_OF_EFFECTIVE_REGRESSORS = 5
SERIES_LEN = 365
SEED = 1
# sample some coefficients
COEFS = np.random.default_rng(SEED).uniform(-1, 1, NUM_OF_EFFECTIVE_REGRESSORS)
trend = make_trend(SERIES_LEN, rw_loc=0.01, rw_scale=0.1)
x, regression, coefs = make_regression(series_len=SERIES_LEN, coefs=COEFS)

# combine trend and the regression
y = trend + regression
y = y - y.min()


x_extra = np.random.normal(0, 1, (SERIES_LEN, NUM_OF_REGRESSORS - NUM_OF_EFFECTIVE_REGRESSORS))
x = np.concatenate([x, x_extra], axis=-1)

x_cols = [f"x{x}" for x in range(1, NUM_OF_REGRESSORS + 1)]
response_col = "y"
dt_col = "date"
obs_matrix = np.concatenate([y.reshape(-1, 1), x], axis=1)
# make a data frame for orbit inputs
df = pd.DataFrame(obs_matrix, columns=[response_col] + x_cols)
# make some dummy date stamp
dt = pd.date_range(start='2016-01-04', periods=SERIES_LEN, freq="1W")
df['date'] = dt
print(df.shape)
print(df.head())<jupyter_output>(365, 12)
          y        x1        x2        x3        x4        x5        x6  \
0  3.192673  0.172792  0.000000  0.165219 -0.000000  0.452678 -0.719291   
1  3.951436  0.223187 -0.000000  0.290559  0.182286  0.147066  1.223912   
2  4.122889  0.014211  0.273356 -0.368227 -0.081455 -0.241060  0.385733   
3  1.669311  0.299423  0.019861 -0.146228 -0.390954 -0.128596 -0.074326   
4  3.567194  0.004071 -0.137801  0.647032  0.000000 -0.000000  0.268331   

         x7        x8        x9       x10       date  
0 -0.167722 -1.066560  2.008808  2.479992 2016-01-10  
1  0.415947  0.731192 -1.029236  0.334582 2016-01-17  
2  0.517487  0.018241 -1.095658 -1.548805 2016-01-24  
3  0.008068  0.253795 -0.710884 -1.144773 2016-01-31  
4  0.469875 -0.571233 -0.410984  1.714833 2016-02-07<jupyter_text>WBIC In this section, we use DLT model as an example. Different DLT models (the number of features used changes) are fitted and their WBIC values are calculated respectively.<jupyter_code>WBIC_ls = []
for k in range(0, NUM_OF_REGRESSORS):
    regressor_col = x_cols[:k + 1]
    dlt_mod = DLT(
        response_col=response_col,
        date_col=dt_col,
        regressor_col=regressor_col,
        seed=2022,
        # fixing the smoothing parameters to learn regression coefficients more effectively
        level_sm_input=0.01,
        slope_sm_input=0.01,
        num_warmup=4000,
        num_sample=4000,
        )
    WBIC_temp = dlt_mod.fit_wbic(df=df) 
    print(f"WBIC value with {k} regressors: {WBIC_temp}")
    print('------------------------------------------------------------------')
    WBIC_ls.append(WBIC_temp)<jupyter_output>WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.
To run all diagnostics call pystan.check_hmc_diagnostics(fit)<jupyter_text>It is also interesting to see if WBIC can distinguish between model types; not just do feature selection for a given type of model. To that end the next block fits an LGT and ETS model to the data; the WBIC values for both models are then calculated.<jupyter_code>lgt = LGT(response_col=response_col,
          date_col=dt_col,
          regressor_col=regressor_col,
          seasonality=52,
          estimator='stan-mcmc',
          seed=8888)
WBIC_lgt = lgt.fit_wbic(df=df) 
print(f"WBIC value for LGT model: {WBIC_lgt}")

ets = ETS(
        response_col=response_col,
        date_col=dt_col,
        seed=2020,
        # fixing the smoothing parameters to learn regression coefficients more effectively
        level_sm_input=0.01,
    )

WBIC_ets = ets.fit_wbic(df=df) 
print(f"WBIC value for ETS model: {WBIC_ets}")

WBIC_ls.append(WBIC_lgt)
WBIC_ls.append(WBIC_ets)<jupyter_output>WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.
To run all diagnostics call pystan.check_hmc_diagnostics(fit)<jupyter_text>The plot below shows the WBIC vs the number of features / model type (blue line). The true model is indicated by the vertical red line. The horizontal gray line shows the minimum (optimal) value. The minimum is at the true value.<jupyter_code>num_reg_ls = list(range(1, NUM_OF_REGRESSORS+1)) + ['LGT','ETS']
plt.plot(num_reg_ls, WBIC_ls, )
plt.axhline(y=min(WBIC_ls), color='grey', linewidth = 1)
plt.axvline(x=NUM_OF_EFFECTIVE_REGRESSORS-1, color='red', linewidth = 3) # because list type is mixed index from 1
plt.ylabel("WBIC")
plt.xlabel("# number of Features / model type")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0992136523339468, 'fraction_numerical': 0.09084825163125314, 'mean_word_length': 3.634108527131783, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '592477', 'license': 'Apache-2.0', 'n_tokens_mistral': 2459, 'n_tokens_neox': 2106, 'n_words': 652, 'path': 'docs/tutorials/wbic.ipynb', 'repo_name': 'juanitorduz/orbit', 'size': '5977'}