Replies: 4 comments 7 replies
-
|
Hi @JustStas Yes, we are still working on a forecast functionality that will need some validation as you mention, but we still need some more time. We still don't have it all figured out so any ideas, feedback, or learnings you've already acquired and want to share, are more than welcome. |
Beta Was this translation helpful? Give feedback.
-
|
@laresbernardo for now I am just trying to put together a script that will be matching the predictions from the model. I am writing it in Python, so not sure it will be helpful, but I will still put the code here :) def adstock(x_ser, theta):
if len(x_ser) == 1:
return 0
else:
return x_ser.iloc[0] + theta * adstock(x_ser.iloc[1:], theta)
class Media:
def __init__(self, media_name, alpha, gamma, theta, beta):
self.name = media_name
self.alpha = alpha
self.gamma = gamma
self.theta = theta
self.beta = beta
def s_curve(self, x_ser):
# TO BE OPTIMIZED (CACHE PRIORS)
decay = adstock(x_ser, self.theta)
return self.beta * decay**self.alpha / (decay**self.alpha + self.gamma**self.alpha)
class Context:
def __init__(self, context_name, beta):
self.name = context_name
self.beta = beta
def lin(self, x):
return self.beta * x
class ProphetPredictor:
def __init__(self, wd_betas, trend_beta, trend_base_date, trend_base_value):
self.wd_betas = wd_betas
self.trend_beta = trend_beta
self.trend_base_date = pd.to_datetime(trend_base_date)
self.trend_base_value =trend_base_value
def date_diff(self, ds):
return (pd.to_datetime(ds) - self.trend_base_date).days
def wd_forecast(self, ds):
wd = pd.to_datetime(ds).day_of_week
return self.wd_betas[wd]
def trend_forecast(self, ds):
days_from_base = self.date_diff(ds)
return self.trend_base_value + days_from_base*self.trend_beta
def full_forecast(self, ds):
return self.wd_forecast(ds) + self.trend_forecast(ds)
class RobynPredictor:
def __init__(self, media_params_df, aggregated_report_df, alldecomp_df, model_v, max_date=None):
self.medias = []
self.contexts = []
if max_date:
alldecomp_df = alldecomp_df.loc[pd.to_datetime(alldecomp_df['ds'])<=max_date]
self.model_v = model_v
self.media_names = [c[:-7] for c in media_params_df.columns if c[-7:]=='_alphas']
self.media_params = media_params_df[media_params_df['solID']==model_v].iloc[0]
self.aggregated_report_df = aggregated_report_df[aggregated_report_df['solID']==model_v]
self.intercept = self.aggregated_report_df.loc[self.aggregated_report_df['rn']=='(Intercept)', 'coef'].iloc[0]
self.context_names = [i for i in self.aggregated_report_df['rn'] if i[-5:]=='_comp' or i[-8:]=='_context']
self.alldecomp_df = alldecomp_df[alldecomp_df['solID']==model_v]
self.x_prior = self.alldecomp_df.loc[:, self.media_names+self.context_names+['ds']].sort_values(by='ds', ascending=False)
self.ingest_all_medias()
self.ingest_all_contexts()
self.ingest_prophet()
def ingest_one_media(self, media_name):
alpha = self.media_params[media_name+'_alphas']
gamma = self.media_params[media_name+'_gammas']
theta = self.media_params[media_name+'_thetas']
beta = self.aggregated_report_df.loc[self.aggregated_report_df['rn']==media_name, 'coef']
media = Media(media_name, alpha, gamma, theta, beta)
self.medias += [media]
def ingest_all_medias(self):
for media_name in self.media_names:
self.ingest_one_media(media_name)
def ingest_one_context(self, context_name):
beta = self.aggregated_report_df.loc[self.aggregated_report_df['rn']==context_name, 'coef']
context = Context(context_name, beta)
self.contexts += [context]
def ingest_all_contexts(self):
for context_name in self.context_names:
self.ingest_one_context(context_name)
def ingest_prophet(self):
last_week = self.alldecomp_df.sort_values(by='ds')[:-7]
last_week['wd'] = last_week['ds'].dt.weekday
trend_base_value = last_week['trend'].iloc[-1]
trend_beta = last_week['trend'].iloc[-1] - last_week['trend'].iloc[-2]
trend_base_date = last_week['ds'].iloc[-1]
wd_betas = dict(zip(last_week['wd'], last_week['weekday']))
self.prophet = ProphetPredictor(wd_betas, trend_beta, trend_base_date, trend_base_value)
def predict_one(self, x):
x = x.sort_values(by='ds', ascending=False) # Row 0 will be the value to predict
s_curves_component = np.sum([m.s_curve(x[m.name]) for m in self.medias])
prophet_component = self.prophet.full_forecast(x['ds'].iloc[0])
context_component = np.sum([c.lin(x[c.name].iloc[0]) for c in self.contexts])
for i, c in enumerate([self.intercept, s_curves_component, prophet_component, context_component]):
print(i, '>', c.round(2))
return self.intercept + s_curves_component + prophet_component + context_component
MODEL_NAME = '1_16_5'
media_params_df = pd.read_csv(os.path.join(MODEL_DIR, 'report_hyperparameters.csv'))
aggregated_report_df = pd.read_csv(os.path.join(MODEL_DIR, 'report_aggregated.csv'))
alldecomp_df = pd.read_csv(os.path.join(MODEL_DIR, 'report_alldecomp_matrix.csv'), parse_dates=['ds'])
p = RobynPredictor(media_params_df, aggregated_report_df, alldecomp_df, MODEL_NAME, max_date='2022-12-31')
full_df = alldecomp_df.loc[alldecomp_df['solID']==MODEL_NAME].sort_values(by='ds', ascending=False)
y_hat = p.predict_one(full_df)For now I have a bug somewhere - my predictions with this method are 7.6 times larger than the model's :) for now I only managed to confirm that the Prophet component matches the model's one, but I have no idea how to compare the s-curves and context ones, since their outputs are not in the CSVs generated by Robyn... |
Beta Was this translation helpful? Give feedback.
-
|
@laresbernardo is there a way to output the adstock-adjusted values from Robyn so that I can compare them with the ones I generate? |
Beta Was this translation helpful? Give feedback.
-
Update on the code:
def adstock(x_ser, theta, d=0):
try:
if len(x_ser) == 1:
return 0
else:
return x_ser.iloc[0] + theta * adstock(x_ser.iloc[1:], theta, d+1)
except RecursionError:
print(d)
class Media:
def __init__(self, media_name, alpha, gamma, theta, beta):
self.name = media_name
self.alpha = alpha
self.gamma = gamma
self.theta = theta
self.beta = beta
def s_curve(self, x_ser):
# TO BE OPTIMIZED (CACHE PRIORS)
# Includes Hill saturation
decay = adstock(x_ser, self.theta)
return self.beta * decay**self.alpha / (decay**self.alpha + self.gamma**self.alpha)
class Context:
def __init__(self, context_name, beta):
self.name = context_name
self.beta = beta
def lin(self, x):
return self.beta * x
class ProphetPredictor:
def __init__(self, wd_betas, trend_beta, trend_base_date, trend_base_value):
self.wd_betas = wd_betas
self.trend_beta = trend_beta
self.trend_base_date = pd.to_datetime(trend_base_date)
self.trend_base_value = trend_base_value
def date_diff(self, ds):
return (pd.to_datetime(ds) - self.trend_base_date).days
def wd_forecast(self, ds):
wd = pd.to_datetime(ds).day_of_week
return self.wd_betas[wd]
def trend_forecast(self, ds):
days_from_base = self.date_diff(ds)
return self.trend_base_value + days_from_base*self.trend_beta
def full_forecast(self, ds):
return self.wd_forecast(ds) + self.trend_forecast(ds)
class RobynPredictor:
def __init__(
self,
media_params_df,
aggregated_report_df,
# alldecomp_df,
pareto_decomp_df,
pareto_clusters,
model_v,
raw_data,
max_date=None
):
self.medias = []
self.contexts = []
# display(alldecomp_df.head())
if max_date:
# alldecomp_df = alldecomp_df.loc[pd.to_datetime(alldecomp_df['ds'])<max_date]
pareto_decomp_df = pareto_decomp_df.loc[pd.to_datetime(
pareto_decomp_df['ds']) < max_date]
raw_data = raw_data.loc[pd.to_datetime(raw_data['ds']) < max_date]
# display(alldecomp_df.head())
self.model_v = model_v
# self.raw_data = raw_data.loc[:, [self.media_names+self.context_names]]
self.pareto_clusters = pareto_clusters[pareto_clusters['solID'] == model_v].iloc[0]
self.media_names = [c[:-7]
for c in media_params_df.columns if c[-7:] == '_alphas']
self.media_params = media_params_df[media_params_df['solID']
== model_v].iloc[0]
self.aggregated_report_df = aggregated_report_df[aggregated_report_df['solID'] == model_v]
self.intercept = self.aggregated_report_df.loc[self.aggregated_report_df['rn']
== '(Intercept)', 'coef'].iloc[0]
self.context_names = [i for i in self.aggregated_report_df['rn']
if i[-5:] == '_comp' or i[-8:] == '_context']
self.pareto_decomp_df = pareto_decomp_df[pareto_decomp_df['solID'] == model_v]
self.raw_data = raw_data.loc[:, self.media_names +
self.context_names+['ds']].sort_values(by='ds', ascending=False)
# self.alldecomp_df = alldecomp_df[alldecomp_df['solID']==model_v]
self.ingest_all_medias()
self.ingest_all_contexts()
self.ingest_prophet()
def ingest_one_media(self, media_name):
alpha = self.media_params[media_name+'_alphas']
gamma = self.media_params[media_name+'_gammas']
theta = self.media_params[media_name+'_thetas']
beta = self.aggregated_report_df.loc[self.aggregated_report_df['rn']
== media_name, 'coef'].iloc[0]
# beta = self.pareto_clusters[media_name]
media = Media(media_name, alpha, gamma, theta, beta)
self.medias += [media]
def ingest_all_medias(self):
for media_name in self.media_names:
self.ingest_one_media(media_name)
def ingest_one_context(self, context_name):
beta = self.aggregated_report_df.loc[self.aggregated_report_df['rn']
== context_name, 'coef']
context = Context(context_name, beta)
self.contexts += [context]
def ingest_all_contexts(self):
for context_name in self.context_names:
self.ingest_one_context(context_name)
def ingest_prophet(self):
last_week = self.pareto_decomp_df.sort_values(by='ds')[:-7]
last_week['wd'] = last_week['ds'].dt.weekday
trend_base_value = last_week['trend'].iloc[-1]
trend_beta = last_week['trend'].iloc[-1] - last_week['trend'].iloc[-2]
trend_base_date = last_week['ds'].iloc[-1]
wd_betas = dict(zip(last_week['wd'], last_week['weekday']))
self.prophet = ProphetPredictor(
wd_betas, trend_beta, trend_base_date, trend_base_value)
def predict_one(self, x):
# Row 0 will be the value to predict
x = x.sort_values(by='ds', ascending=False)
for m in self.medias:
print(m.name, m.beta, '>', m.s_curve(x[m.name]).round(2))
s_curves_component = np.sum([m.s_curve(x[m.name])
for m in self.medias])
prophet_component = self.prophet.full_forecast(x['ds'].iloc[0])
context_component = np.sum(
[c.lin(x[c.name].iloc[0]) for c in self.contexts])
for i, c in enumerate([self.intercept, s_curves_component, prophet_component, context_component]):
print(i, '>', c.round(2))
return self.intercept + s_curves_component + prophet_component + context_component
MODEL_DIR = '/Users/stanislav.nosulenko/Documents/MMM/WIP/robyn_modelling/20220421_v2/2022-04-21 11.52 init'
MODEL_NAME = '1_24_6'
media_params_df = pd.read_csv(os.path.join(
MODEL_DIR, 'pareto_hyperparameters.csv'))
aggregated_report_df = pd.read_csv(
os.path.join(MODEL_DIR, 'pareto_aggregated.csv'))
pareto_decomp_df = pd.read_csv(os.path.join(
MODEL_DIR, 'pareto_alldecomp_matrix.csv'), parse_dates=['ds'])
raw_data = pd.read_csv(
'/Users/stanislav.nosulenko/Documents/MMM/WIP/robyn_modelling/20220421_robyn_df.csv', parse_dates=['ds'])
pareto_clusters = pd.read_csv(os.path.join(
MODEL_DIR, 'pareto_clusters.csv'))I match all components of yhat in pareto_all_decomp_matrix.csv, except the s-curves. If anyone has ideas on what is going wrong in my code - I will be grateful for the advice. |
Beta Was this translation helpful? Give feedback.
-
|
OK, new update: |
Beta Was this translation helpful? Give feedback.
-
|
OK, new update (I think I am getting closer :)). The issue is probably with the gamma - I am using the one from pareto_hyperparameters.csv, while the formula from the link above refers to a "transformed gamma". |
Beta Was this translation helpful? Give feedback.
-
|
I updated the s-curve computation: class Adstock:
def __init__(self, theta):
self.theta = theta
self.decays = []
def compute(self, x_ser, d=0):
try:
if len(x_ser) == 1:
return 0
else:
decay = self.compute(x_ser.iloc[1:], d+1)
self.decays += [decay]
return x_ser.iloc[0] + self.theta * decay
except RecursionError:
print(d)
def transform_gamma(self, gamma):
return np.quantile(self.decays, gamma)
def s_curve(adstocked, alpha, gamma):
return adstocked**alpha / (adstocked**alpha + gamma**alpha)
class Media:
def __init__(self, media_name, alpha, gamma, theta, beta):
self.name = media_name
self.alpha = alpha
self.gamma = gamma
self.theta = theta
self.beta = beta
def transform(self, x_ser):
# TO BE OPTIMIZED (CACHE PRIORS)
# Includes Hill saturation
self.adstock = Adstock(self.theta)
ads_value = self.adstock.compute(x_ser)
self.gamma_transformed = self.adstock.transform_gamma(self.gamma)
adstocked = adstock(x_ser, self.theta)
s_curved = s_curve(adstocked, self.alpha, self.gamma_transformed)
# print('\nbeta', self.beta)
# print('alpha', self.alpha)
# print('gamma', self.gamma)
# print('gamma_transformed', self.gamma_transformed)
# print('decay', adstocked)
# print('decay alpha', adstocked**self.alpha)
# print('s curve', s_curved)
# print('return', self.beta * s_curved)
return self.beta * s_curved
class Context:
def __init__(self, context_name, beta):
self.name = context_name
self.beta = beta
def transform(self, x):
return self.beta * x
class ProphetPredictor:
def __init__(self, wd_betas, trend_beta, trend_base_date, trend_base_value):
self.wd_betas = wd_betas
self.trend_beta = trend_beta
self.trend_base_date = pd.to_datetime(trend_base_date)
self.trend_base_value = trend_base_value
def date_diff(self, ds):
return (pd.to_datetime(ds) - self.trend_base_date).days
def wd_forecast(self, ds):
wd = pd.to_datetime(ds).day_of_week
return self.wd_betas[wd]
def trend_forecast(self, ds):
days_from_base = self.date_diff(ds)
return self.trend_base_value + days_from_base*self.trend_beta
def full_forecast(self, ds):
return self.wd_forecast(ds) + self.trend_forecast(ds)with a transformation of Gamma (as I understand it). Unfortunately, the s-curve results still don't match the corresponding components of the model... |
Beta Was this translation helpful? Give feedback.
-
|
All right, I think I made it work: class Adstock:
def __init__(self, theta):
self.theta = theta
self.decays = []
def compute(self, x_ser, d=0):
try:
if len(x_ser) == 1:
return 0
else:
decay = self.compute(x_ser.iloc[1:], d+1)
self.decays += [decay]
return x_ser.iloc[0] + self.theta * decay
except RecursionError:
print(d)
def transform_gamma(self, gamma):
min_d = np.min(self.decays)
max_d = np.max(self.decays)
lin_sp = np.linspace(min_d, max_d, 100)
# return np.quantile(self.decays[:100], gamma)
return round(np.quantile(lin_sp, gamma), 4)
def s_curve(adstocked, alpha, gamma):
return adstocked**alpha / (adstocked**alpha + gamma**alpha)
class Media:
def __init__(self, media_name, alpha, gamma, theta, beta):
self.name = media_name
self.alpha = alpha
self.gamma = gamma
self.theta = theta
self.beta = beta
def transform(self, x_ser):
# TO BE OPTIMIZED (CACHE PRIORS)
# Includes Hill saturation
self.adstock = Adstock(self.theta)
ads_value = self.adstock.compute(x_ser)
self.gamma_transformed = self.adstock.transform_gamma(self.gamma)
# adstocked = self.adstock(x_ser, self.theta)
s_curved = s_curve(ads_value, self.alpha, self.gamma_transformed)
# print('\nbeta', self.beta)
# print('alpha', self.alpha)
# print('gamma', self.gamma)
# print('gamma_transformed', self.gamma_transformed)
# print('decay', adstocked)
# print('decay alpha', adstocked**self.alpha)
# print('s curve', s_curved)
# print('return', self.beta * s_curved)
return self.beta * s_curved
class Context:
def __init__(self, context_name, beta):
self.name = context_name
self.beta = beta
def transform(self, x):
return self.beta * x
class ProphetPredictor:
def __init__(self, wd_betas, trend_beta, trend_base_date, trend_base_value):
self.wd_betas = wd_betas
self.trend_beta = trend_beta
self.trend_base_date = pd.to_datetime(trend_base_date)
self.trend_base_value = trend_base_value
def date_diff(self, ds):
return (pd.to_datetime(ds) - self.trend_base_date).days
def wd_forecast(self, ds):
wd = pd.to_datetime(ds).day_of_week
return self.wd_betas[wd]
def trend_forecast(self, ds):
days_from_base = self.date_diff(ds)
return self.trend_base_value + days_from_base*self.trend_beta
def full_forecast(self, ds):
return self.wd_forecast(ds) + self.trend_forecast(ds)for some reason, the gamma transformation in the initial code looks like this: |
Beta Was this translation helpful? Give feedback.
-
|
Hey @JustStas I'm Cristian Nozzi from facebook. I haven't had the chance to test this out yet but In the meantime I wanted to share a thought with you. The idea is super cool and I see how this can be really helpful but in your specific case, the model you build with 24 months of data and validate with the 25th month of data will only be good for retro-active analysis. If you provide the 25th month of data using the Refresh function you'll end up with a new set of data. Hope I was clear enough, what do you think about it? |
Beta Was this translation helpful? Give feedback.
-
|
Hi @Tsurty !
The reason why I think it would work is because on each iteration we are doing incremental updates to a model that proves itself in the "real" world (aka validation data). Thus, hopefully, each consecutive month will have a strong base to build and improve upon (Side note: I would keep the refresh iterations # to a relatively low level to avoid overtraining for each month). The inspiration for such an approach comes from evolutionary neural networks: https://www.youtube.com/watch?v=aeWmdojEJf0 Hopefully, such an approach would be less prone to overfitting - otherwise, the choice of a model using solely our understanding of the business seems like a significant reduction of the value addition of the model... What do you think? |
Beta Was this translation helpful? Give feedback.
-
I mean the concept of validating over an holdout set is something I 100% agree with as we do that everytime with other kinds of models. Anyway it's something to simply test out and validate but the idea is really cool :) |
Beta Was this translation helpful? Give feedback.
-
Hi everyone!
Could you please tell me what is the best way to do out-of-sample validation/testing of Robyn models? I would like to see how a trained model performs on a holdout set.
I saw the refresh functionality, but if I understand correctly, it does not show how the non-updated model performs on the new data...
Is there a way to make predictions based on the CSVs generated by the model?
Beta Was this translation helpful? Give feedback.
All reactions