Has anyone developed a custom tweedie distribution for sampling? How to implement once developed?

Questions
6

The problem is that you’re trying to pass a dictionary as observed data. I could suggest two things. Firstly, not using numpy on your likelihood function, but using pytensor instead. Secondly, use a pm.Potential instead of pm.DensityDist. Here’s the idea:

import pymc as pm
import pytensor.tensor as pt

# Define the log probability function for the Tweedie distribution
def tweedie_logp(name, value, mu, p, phi):
    # Ensure positive values for mu, p, and phi
    mu = pt.maximum(mu, 1e-10)
    p = pt.maximum(p, 1e-10)
    phi = pt.maximum(phi, 1e-10)

    # Calculate the log probability using the provided parameters and observed value
    logp_value = (
        (1 - p) * pt.log(value / mu) - (value / (mu ** p * phi)) - pt.log(phi)
        + pt.log((1 - p) * mu ** (1 - p))
    )

    return logp_value.sum()


with pm.Model() as model:
    # Prior distributions for the unknown parameters

    # Prior distributions for the unknown parameters
    mu = pm.Normal('mu', mu=0, sigma=10)
    p = pm.Normal('p', mu=1.5, sigma=1)
    phi = pm.Gamma('phi', alpha=1, beta=1)

    # Priors for the binary variables
    offline_campaign_prior = pm.Bernoulli('MISX_campaign_prior', p=0.5)
    discounts_prior = pm.Bernoulli('PA_prior', p=0.5)
    promotion_events_prior = pm.Bernoulli('promotion_events_prior', p=0.5)

    # Compute the expected sales based on the variables
    expected_sales = (storedata['day_of_week'].values + storedata['month'].values + storedata['year'].values + storedata['MISX'].values * offline_campaign_prior +
                      storedata['Status'].values * discounts_prior +
                      storedata['event'].values * promotion_events_prior)
    

    # Define the observed sales as a data variable
    sales_observed = storedata['Rslr_Sales_Amount'].values


    # Define the likelihood using the custom Tweedie distribution
    sales = pm.Potential("sales", lp("sales_lp", sales_observed, expected_sales, p, phi)) 
    

    # Define the model
    trace = pm.sample(1000, tune=1000)

I haven’t tested it, as I don’t have the data or appropriate simulations for it. But hopefully helps.

1 Like