Pymc5 out of sample

Questions v5
5

Hi Jesse, sorry for the late reply. First of all thank you for your help
Unfortunately, I tried applying your code to my case but didn’t get satisfactory results. I attempted to simplify the number of features and controls, but the two models yield extremely different results. I’ll try to write both codes below and also attach the sample database.
PS I will definitely read the pymc marketing package :slight_smile:

My code:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import arviz as az
import pymc as pm
import pytensor.tensor as tt 



print(f"Running on PyMC v{pm.__version__}")


#decay
def geometric_adstock_tt(x_t, alpha=0, L=21, normalize=True):
    w = tt.as_tensor_variable(
        [tt.power(alpha, i) for i in range(L)]
    )
    
    xx = tt.stack(
        [tt.concatenate([
            tt.zeros(i),
            x_t[:x_t.shape[0]-i]
        ]) for i in range(L)]
    )
    
    if not normalize:
        y=tt.dot(w, xx)
    else:
        y=tt.dot(w/tt.sum(w),xx)
    
    return y


#logistic
def logistic_function(x_t, mu=0.1):
    return (1-np.exp(-mu*x_t))/(1+np.exp(-mu*x_t))




RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")



# import data
df = pd.read_csv ('scaled.csv', sep=";")
df['Period']=pd.to_datetime(df['Period'])
df = df.set_index('Period')

features = ['feature1', 'feature2']
control_vars = ['control1'] 
target = df['y'].to_numpy()

# bayes model
basic_model =pm.Model()

with basic_model:
    response_mean = []
    
    for feature in features:
        xx = df[feature].values
        print (f'Adding media: {feature}')
        beta = pm.HalfNormal(f'beta_{feature}', sigma = 2)
        decay = pm.Beta(f'decay_{feature}', alpha=3, beta=3)
        sat = pm.Gamma(f'sat_{feature}', alpha=3, beta=1)
        feature_contribution = pm.Deterministic(f'contribution_{feature}', (logistic_function(geometric_adstock_tt(xx, decay),sat)*beta))
        response_mean.append(feature_contribution)
        
        
    
    if control_vars:
        for control in control_vars:
            x = df[control].values
            
            print (f'Adding control: {control}')
            control_beta = pm.Normal(f'coeff_{control}', sigma = 2)
            
            control_contribution = pm.Deterministic(
            f'contribution_{control}', (control_beta * x))
            
            #channel_contrib = control_beta * x
            response_mean.append(control_contribution)
            
   
            
   
    sigma = pm.HalfNormal('sigma', sigma=1)
    likelihood =  pm.Normal("likelihood", mu=sum(response_mean), sigma=sigma, observed=target)
        

with basic_model:
    # draw 1000 posterior samples
    idata = pm.sample(return_inferencedata=True, tune= 1000)
    

with basic_model:
    post = pm.sample_posterior_predictive(idata)


arr = post.posterior_predictive['likelihood']
mean_y = np.mean(arr, axis=(0,1))
mean_y =mean_y.values

Here below your code adjusted by me:


import pandas as pd
import numpy as np
import arviz as az
import pymc as pm
import pytensor.tensor as tt 



print(f"Running on PyMC v{pm.__version__}")


#decay
def geometric_adstock_tt(x_t, alpha=0, L=21, normalize=True):
    w = tt.as_tensor_variable(
        [tt.power(alpha, i) for i in range(L)]
    )
    
    xx = tt.stack(
        [tt.concatenate([
            tt.zeros(i),
            x_t[:x_t.shape[0]-i]
        ]) for i in range(L)]
    )
    
    if not normalize:
        y=tt.dot(w, xx)
    else:
        y=tt.dot(w/tt.sum(w),xx)
    
    return y


#logistic
def logistic_function(x_t, mu=0.1):
    return (1-np.exp(-mu*x_t))/(1+np.exp(-mu*x_t))




RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")



# import data
df = pd.read_csv ('scaled.csv', sep=";")
[scaled.csv|attachment](upload://aVzn3kWmmO8eE3bhm2o6vEOynI7.csv) (10.1 KB)

df['Period']=pd.to_datetime(df['Period'])
df = df.set_index('Period')

features = ['feature1', 'feature2']
control_vars = ['control1'] 
target = df['y'].to_numpy()
df_features  = df[features]
n_obs, n_features = df_features.shape
coords = {'features':features, 
          'all_vars':features}


if control_vars:
    df_controls  = df[control_vars]
    _, n_controls = df_controls.shape
    coords.update({'controls':control_vars,
                   'all_vars':features + control_vars})

with pm.Model(coords=coords) as basic_model:
    X = pm.MutableData('feature_data', df_features.values)
    y = pm.MutableData('targets', df['y'].values)

    n_obs = X.shape[0]
    
    betas = pm.HalfNormal('beta', sigma = 2, dims=['features'])
    decays = pm.Beta('decay', alpha=3, beta=3, dims=['features'])
    sat = pm.Gamma('sat', alpha=3, beta=1, dims=['features'])
    contributions = tt.zeros((n_obs, n_features + n_controls))
    
    
    
    
    for i in range(n_features):
       x = logistic_function(geometric_adstock_tt(X[:, i], decays[i]),sat[i])*betas[i]
       contributions = tt.set_subtensor(contributions[:, i], x)
   
    if n_controls > 0:
        Z = pm.MutableData('control_data', df_controls.values)
        control_betas = pm.Normal('control_beta', sigma = 2, dims=['controls'])
        tt.set_subtensor(contributions[:, n_features:],  Z @ control_betas)
    
    pm.Deterministic("contributions", contributions)
     
    sigma = pm.HalfNormal('sigma', sigma=1)
    
    y_hat=  pm.Normal("y_hat", mu=contributions.sum(axis=-1), sigma=sigma, observed=y, shape=X.shape[0])
    

with basic_model:
    # draw 1000 posterior samples
    idata = pm.sample(idata_kwargs={'dims':{'contributions':[None, 'all_vars']}})
    

with basic_model:
    post = pm.sample_posterior_predictive(idata)
    
arr = post.posterior_predictive['y_hat']
mean_y = np.mean(arr, axis=(0,1))
mean_y =mean_y.values

Every suggestion is greatly appreciated.
All the best