Introduction to Bayesian Data Analysis - Applied Version¶

GOOD_REVIEWS = np.random.binomial(666, 0.52, 100)

def build_model(observed):
    """
    Function for building a PyMC3 model
    
    Arguments:
        observed: value or array of values, data to fit the model on
    
    Returns:
        estimation_model: A PyMC3 model
    """
    with pm.Model() as esitmation_model:

        # prior distribution of the paramter to estimate
        p = pm.Beta(
            alpha=3,
            beta=3,
            name='proportion'
        )

        # prior distribution of the paramter to estimate
        n = pm.Poisson(
            mu=666,
            name='total reviews'
        )

        # the outocme distribution
        outcome = pm.Binomial(
            p=p,
            observed=GOOD_REVIEWS,
            n=n,
            name='good reviews'
        )
    
    return esitmation_model

with build_model(GOOD_REVIEWS):
    
    # we sample from the prior distribution of 
    # the paramters
    prior_checks = pm.sample_prior_predictive(
        samples=1000
    )

# sklearn will do "bayesian" fitting 
# but will return point estimate for the parameters
from sklearn.linear_model import BayesianRidge

ridge_model = BayesianRidge(
    alpha_1=0.01, # prior Gamma
    alpha_2=0.01, # prior Gamma
    lambda_1=0.01, # prior Gamma
    lambda_2=0.01 # prior Gamma
)

X = np.random.random(100)
y = np.random.random(100)
# PyMC3 will return distribution for the parameters
with pm.Model() as ridge_model:
    # I decide that my intercpet comes from a normal distribution
    intercept = pm.Normal(
        'Intercept',
        mu=0.,
        sigma=1.
    )
    # I decide which type of prior I want to put on the slope
    slope_sigma = pm.Gamma(
        'Slope Sigma',
        alpha=0.01,
        beta=0.01
    )
    slope = pm.Normal(
        'Slope',
        mu=0,
        sd=slope_sigma ,
    )
    mu = intercept + slope*X
    sigma = pm.HalfCauchy(
        'sigma',
        beta=10.
    )
    # I decide that my outcome is normally distributed
    outcome = pm.Normal(
        'y',
        mu=mu,
        sd=sigma,
        observed=y
    )

from modules.visuals import visualize_bivariate_relationship

X, y, coef = make_regression(
    n_samples=100,
    n_features=1,
    bias=200,
    noise=50,
    coef=True,
    random_state=56
)
X = (X - X.mean()) / X.std()
X = X.flatten()

prediction = ar_1_model.predict(X_ts, y_ts, verbose=False)

visualize_time_series(
    X=[i for i in range(len(process))], 
    y=process, 
    prediction=prediction, 
    prediction_point=36, 
    X_label='Months', 
    y_label='Number of Visits',
    figsize=(8, 8)
)

from modules.utils.data_utils import generate_game_difficulty_data

df, true_parameters = generate_game_difficulty_data(
    players=10,
    levels=5,
    n_sessions=100
)
df.head(10)

visualize_bivariate_parameter_grid(
    parameter_1=true_parameters['level_difficulty'], 
    parameter_2=true_parameters['player_ability'],
    parameter_1_name='Level Difficulty', 
    parameter_2_name='Player Ability',
    height=7
)

from modules.visuals import visualize_heatmap

visualize_heatmap(
    df=true_parameters, 
    pivot_varaibles=('player_ability', 'level_difficulty', 'probability_success'), 
    rounding=4, 
    figsize=(8, 8),
)

visualize_bivariate_relationship(
    X=true_parameters['delta'],
    y=true_parameters['probability_success'],
    X_label='$\Delta$(Difficulty, Ability)',
    y_label='Probability Success',
    title=f'Probability Succeding Level as Function Difficulty - Ability',
    c=true_parameters['probability_success'],
    cmap='rocket'
)

from modules.stats.models.graphical_models import GameDifficultyModel

PLAYERS_ID = df['player_id'].values
LEVELS_ID = df['level'].values
ATTEMPTS = df['num_attempts'].values
SUCCESSES = df['num_success'].values

game_difficulty_model = GameDifficultyModel(
    players_id=PLAYERS_ID,
    levels_id=LEVELS_ID,
    num_attempts=ATTEMPTS,
    num_success=SUCCESSES
)

game_difficulty_model.show_plate()

game_difficulty_model.fit(
    MAP=False,
    draws=3000, 
    tune=2000,
    chains=2,
    cores=1,
    target_accept=.90
)

game_difficulty_model.show_posterior_summary(
    parameters_name=[
        'player_ability_mu',
        'hyper_sigma',
        'level_difficulty_mu',
        'player_ability',
        'level_difficulty'
    ],
    compact=True
)

game_difficulty_model.show_forest_plot(
    parameters_name=['player_ability']
)

game_difficulty_model.show_forest_plot(
    parameters_name=['level_difficulty']
)

from modules.utils.data_utils import generate_polynomial_data
from modules.stats.models.linear_models import PolynomialRegression


df, true_parameters = generate_polynomial_data(
    X=[i for i in range(100)],
    degree=2,
    noise_ratio=0.1
)

X = df['X'].values
y = df['y'].values

visualize_bivariate_regression(
    X=X,
    y=y,
    X_label='Predictor',
    y_label='Outcome',
    title=f'Quadratic Relationship',
)

models = {
    'Linear': BivariateNormalRegression(
        intercept_prior=(0, 1),
        slope_prior=(0, 1),
        likelihood_sigma_prior=1,
        X=X,
        y=y
    ),
    'Quadratic': PolynomialRegression(
        X=X,
        y=y,
        cubic=False
    ),
    'Cubic': PolynomialRegression(
        X=X,
        y=y,
        cubic=True
    )
} 
for name, model in models.items():
    
    print(f'Sampling from {name}')
    model.fit(
        MAP=False, 
        draws=1000, 
        tune=1000,
        chains=2,
        cores=1,
        target_accept=.90,
        init='adapt_diag'
    )

df_results = pm.compare(
    {name: model.get_traces() for name, model in models.items()},
    ic='loo',
    scale='log'
)
df_results

pm.plot_compare(
    df_results,
    figsize=(10, 5)
)