Simple way to making a prediction using PyMC model

Questions
1

Hi,

I am new to PyMC but am trying to define my first model. The problem I would like to solve is as follows:

Buyers are viewing the products, but the potential purchase of the viewed products may be postponed for a few days (maximum 20 days). Moreover, the amount they pay when buying may be different than the price of the product at the time of viewing (because they later decide to buy more products at once or buy the same product on sale, etc.). I’d like to model the seller’s total revenue based on today’s views.

My training data includes past transaction and revenues data, but I would like to predict the total revenue for today.

I defined the simple model below:

# Vector of 0s and 1s - 1s when view ends with purchase (transaction), 0 otherwise
t = train['was_transaction'] 

# Vector of differences between revenue and price (defined only for views with purchase)
residues = (train['Revenue'] - train['Price'])[train['was_transaction']]

coords = {"view_index": np.arange(len(t)), "transaction_index": np.arange(len(residues))}

with pm.Model(coords=coords) as pooled_model:
    observed_transactions = pm.Data("observed_transactions", t, dims="view_index")
    observed_revenues = pm.Data("observed_revenue", residues, dims="transaction_index")

    transaction_probability = pm.Normal("transaction_probability", mu=0.01, sd=0.1)
    was_transaction = pm.Bernoulli(
        "was_transaction",
        p=transaction_probability,
        dims="view_index",
        observed=observed_transactions
    )

    mu = pm.Normal("mu", mu=70, sd=30)
    revenue_residue = pm.Normal(
        "revenue_residue",
        mu=mu,
        sigma=100,
        dims="transaction_index",
        observed=observed_revenues
    )

After sampling, I have estimates for transaction_probability and mu. But now I would like to predict the overall revenue for today. I did this by defining a few additional variables:

with pooled_model:
    viewed_product_prices = pm.Data("prices", test['Price'])

    was_transaction_predict = pm.Bernoulli(
        "was_transaction_predict",
        p=transaction_probability,
        shape=len(test)
    )

    revenue_residue_predict = pm.Normal(
        "revenue_residue_predict",
        mu=mu,
        sigma=100,
        shape=len(test)
    )

    overall_revenue = pm.Deterministic("overall_revenue", tt.sum(was_transaction_predict * (revenue_residue_predict + viewed_product_prices)))

    posterior_predictive = pm.sample_posterior_predictive(
        traces, var_names=["overall_revenue"], random_seed=999
    )

This way I have posterior for overall_revenue.

I don’t think this is the best way to make prediction. Is there any simpler way (without the need to define additional variables “was_transaction_predict” and “revenue_residue_predict”, which are duplicates of “was_transaction” and “revenue_residue”)? Maybe I can formulate the model differently?

Every remark is valuable to me.

Thanks in advance!

2 Likes
2

Welcome!

A couple of preliminary points. First, transaction_probability is normally distributed (with support from negative infinity to positive infinity), but this can cause problems when used in Bernoulli. Conventionally, one would use a beta distribution on p so that support would be on p \in [0,1].

Second, your prediction code seems to add things to your model that may not be entirely needed (e.g., you add a bunch of variability to revenue_residue_predict for reasons that I can’t quite grasp). It seems like it would be more straightforward to just use the samples from your posterior directly. Something along these lines:

with pm.Model(coords=coords) as pooled_model:
    #model definition goes here

    # sample
    idata = pm.sample(return_inferencedata=True)

# unpack sampled parameter values
p = idata.posterior['transaction_probability'].values.ravel()
mu = idata.posterior['mu'].values.ravel()

# combine sampled values with prices
expected_rev = np.sum([p * (mu * price) for price in test['Price']])

There are definitely other ways of doing it, but maybe that will help?

3

Thank you for yours response!

Thank you also for your remark about Beta distribution, it definitely improves the model.

I would like to make sure - if I need the distibution of revenue, not just expected revenue, I can do the following:

# p.shape = (2000, 1)
p = idata.posterior['transaction_probability'].values.ravel()[:, np.newaxis]
# mu.shape = (2000, 1)
mu = idata.posterior['mu'].values.ravel()[:, np.newaxis]

# overall_revenue_samples.shape = (2000,)
overall_revenue_samples = (p * (mu + test['Price'].values[np.newaxis, :])).sum(axis=1)

Then, I can visualize overall_revenue_samples with, for example, arviz plot_posterior. Am I right?

1 Like
4

Sure. Or you could just sum across prices, but keep the sum separate for each sample in your posterior:

posterior_rev = np.sum([p * (mu * price) for price in test['Price']], axis=0)

Then I would probably feed this to arviz’s plot_dist().