Modelling sales for two cities

Questions v5
,
8

From what you are saying and the data you are showing, I think that you mean that you want the random fluctuations of both cities to be correlated with each other. To do this, you need to change the last three lines of your model.

Your current model says:

    sigma = pm.HalfNormal("sigma", sigma=200)

    # fit two cities
    pm.Normal(name="l_city_a", mu=mu_city_a, sigma=sigma, observed=sales_city_a)
    pm.Normal(name="l_city_b", mu=mu_city_b, sigma=sigma, observed=ftr_ios)

This basically means that you assume that the observations for city a and city b come from independent normal distributions. This means that, without knowing, you baked into your model the assumption that the fluctuations across cities should be uncorrelated!

To fix this, you should change your observed distribution to something that can handle correlations across observations. I suggest that you use an MvNormal and also take the opportunity to actually learn how correlated each city is.

    # Your old sigma now is replaced by an LKJ prior for the correlation
    # plus independent sigmas for each city, both of which are assumed to come
    # from a HalfNormal with 200 standard deviation
    chol, corr, sigma = pm.LKJCholeskyCov(
        "packed_L", n=2, eta=1.0, sd_dist=pm.HalfNormal.dist(200), compute_corr=True
    )
    
    # fit two cities
    pm.MvNormal(
        name="l_city_a",
        mu=[mu_city_a, mu_city_b],
        chol=chol,
        observed=np.stack([sales_city_a, ftr_ios], axis=-1)
    )

Let us know how this works out. I can’t test it without the actual data.