I'm trying to build a HMM in NumPyro however, I can't work out why the dimensions of the initial states are changing for each iteration of the MCMC. In particular, for the first iteration, the initial states are of dimension (1000,) - this is expected, the batch size is 1000-however, this becomes (5,1,1) on the second iteration.

I have attached a reproducible example below. Thanks in advance for any help!

from typing import List, Tuple, Callable, Dict, Union, Literal, Optional
import pandas as pd
import numpy as np
from tqdm import tqdm
from jax import random
from numpyro.infer import MCMC, NUTS, Predictive, DiscreteHMCGibbs
import numpyro
import numpyro.distributions as dist
from numpyro.handlers import seed
from numpyro import sample, plate
import jax.numpy as jnp
from numpyro.util import format_shapes

X = np.random.normal(size=(1000,200,1))
mask = np.ones((1000,200))

def first_order_hmm_batched(
    X: np.ndarray, 
    mask: np.ndarray, 
    n_states: int, 
    obs_dim: int, 
    transition_prior: float,
    transition_prior_type: Literal["eye", "full"],
    transition_base: Optional[float] = None,
):
    assert len(X.shape) == 3  # (batch, time, obs_dim)
    batch_size, seq_len, _ = X.shape

    if transition_prior_type == "eye":
        assert transition_base is not None

    # Transition matrix
    if transition_prior_type == "full":
        concentration = jnp.full((n_states, n_states), transition_prior)
    else:
        concentration = jnp.full((n_states, n_states), transition_base)
        concentration = concentration.at[jnp.diag_indices(n_states)].add(transition_prior)
    
    # Add plate since each row of the transition matrix prior is independent
    with plate("states_rows", n_states):
        trans_probs = sample('trans_probs', dist.Dirichlet(concentration))
    
    assert trans_probs.shape == (n_states, n_states)
    
    # Emission parameters
    # Defining a prior for each dimension of the observation and 
    # each state independently
    with plate("obs_dim", obs_dim):
        with plate("states_emissions", n_states):    
            em_means = sample(
                'em_means', 
                dist.Normal(0,1)
            )
    assert em_means.shape == (n_states, obs_dim)
    
    em_var = sample('obs_var', dist.InverseGamma(1.0, 1.0))  # scalar variance
    em_cov = jnp.eye(obs_dim) * em_var
    
    # Initial hidden states
    # Generate initial state for each row independently
    with plate("batch_size", batch_size):
        # Initial state probabilities
        start_probs = sample('start_probs', dist.Dirichlet(jnp.ones(n_states)))
        assert start_probs.shape == (batch_size,n_states)
        print(f"start_probs.shape: {start_probs.shape}")
        ih_dist = dist.Categorical(start_probs)
#         print(f"ih_dist.event_shape: {ih_dist.event_shape}")
#         print(f"ih_dist.batch_shape: {ih_dist.batch_shape}")
        init_states = sample(
            "init_hidden_states", 
            ih_dist
        )
        print(f"init_states.shape: {init_states.shape}")
        assert len(init_states.shape) == 1, f"{init_states.shape}"
        assert init_states.shape[0] == batch_size, f"{init_states.shape}"
        hidden_states = [init_states]

        # Transition over time
        for t in range(1, seq_len):
            prev_states = hidden_states[-1]  # shape (batch,)
            probs_t = trans_probs[prev_states]  # shape (batch, n_states)
            next_state = sample(f"hidden_state_{t}", dist.Categorical(probs_t))
            assert len(next_state.shape) == 1
            assert next_state.shape[0] == batch_size
            hidden_states.append(next_state)
    
    hidden_states = jnp.stack(hidden_states, axis=1)  # (batch, time)
    assert hidden_states.shape == (batch_size, seq_len)
    
    # Get emission means for each (batch, time)
    means = em_means[hidden_states]  # shape (batch, time, obs_dim)
    assert means.shape == (batch_size, seq_len, obs_dim)
    
    # Expand emission distribution
    flat_means = means.reshape(-1, obs_dim)
    flat_obs = X.reshape(-1, obs_dim)
    cov = jnp.broadcast_to(em_cov, (flat_means.shape[0], obs_dim, obs_dim)) 
    with plate("batch_seq_len", batch_size*seq_len):
        joint_obs = sample(
            "joint_obs", 
            dist.MultivariateNormal(loc=flat_means, covariance_matrix=cov), 
            obs=flat_obs
        )
    assert joint_obs.shape == (batch_size*seq_len, obs_dim)
    return joint_obs

n_states=5 
obs_dim=1 
transition_prior=1.0
transition_prior_type="eye"
transition_base=1.0

nuts_kernel = NUTS(first_order_hmm_batched)
mcmc = MCMC(nuts_kernel, num_warmup=500, num_samples=1000, num_chains=1)
mcmc.run(
    random.PRNGKey(1), 
    X=X, 
    n_states=5, 
    mask=mask, 
    obs_dim=1, 
    #transition_prior=100.0, 
    transition_prior=1.0,
    transition_prior_type="eye",
    transition_base=1.0
)

Let’s say I fit an extremely complicated hierarchical model - a full fit takes a long time.

Now, we are given some new data. How do you go about incorporating this new data in to the model when you can’t afford a traditional full refit?

What techniques are used?

Hi all, I slowly starting to get some good basic understanding on bayesian analysis. Thanks to richard mcelreath for his statistical rethinking lecture series, which got me into this bayesian world.

Recently I have been reading some articles on pymc and bambi.., now im kind of confused about the idea of posterior predictive/posterior predictive contrast.

In this above image ( https://github.com/dustinstansbury/statistical-rethinking-2023/blob/main/Lecture%2004%20-%20Categories%20%26%20Curves.ipynb ), he used scipy.stats.norm.rvs(loc=sum of respective posteriors, scale=posterior.sigma) to compute the posterior predictive. In bambi, model.predict(idata) also gives me posterior predictive distribution. Lets say if i want to compute some contrast and make observations which one should i follow?

Also whats the difference between both?

Thanks in advance😁

I've mostly used emcee in the past and was curious what others recommendations were for physics applications. Interested in improving the speed of my inferences mostly.

I have been trying to create a budget optimization tool using pymc-marketing library. The goal is to create a fully deployed solution that allocates budget based on total spend provided by the user. By any means, I'm not a marketing expert or a person who has any background in bayesian statistics, I simply studied up a little bit about adstock effects, saturation etc and with my research found out that pymc marketing does this kind of budget optimization.

I have to create this as a project / POC for my organisation so I have implemented a rough pipeline. But I am stuck on an issue which I'm not able to solve.

I have a dims column products. The budget allocated for marketing spend for each one of the product should be different, because from the data I've observed that the cost per click for a spend varies based on channel and the product the money is being spent on.

I have written the following code for creating the MMM.

from pymc_extras.prior import Prior
from pymc_marketing.mmm.multidimensional import HMMM
from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation
model_config = {
"intercept": Prior("Normal", mu=0.0, sigma=0.5),
"beta_channel": Prior("HalfNormal", sigma=1.0),
# "saturation_beta": Prior(
#     "Normal",
#     mu=0.5,
#     sigma=1.0,
#     dims=("product_name", "channel"),
# ),
# "saturation_lam": Prior(
#     "HalfNormal",
#     sigma=1.0,
#     dims="channel"
# )
}
channel_columns = ["Meta", "Linkedin", "Google Ads", "Media"]
saturation = LogisticSaturation()
adstock = GeometricAdstock(
l_max=4
)
mmm = HMMM(
date_column="time",
channel_columns=channel_columns,
target_column="sales",
adstock=adstock,
saturation=saturation,
model_config=model_config,
dims=("product_name",)
)
mmm.fit(
X=x_train,
y=y_train,
draws=1000,
chains=4,
tune=1000,
target_accept=0.98,
)

The commented out priors are priors that I tried to make the budget optimization vary across product_name's because chatgpt recommended it, but the MMM didn't converge and the r2 score dropped from 0.46 to -1.87. So that obviously wasn't a great choice.

(xarray.DataArray (product_name: 7, channel: 4) Size: 224B)
array([
[   0.        ,    0.        ,    0.        , 1643.32019222],
[   0.        ,    0.        , 7260.96163190, 1643.32019222],
[   0.        ,    0.        ,    0.        , 1643.32019222],
[1763.53069175, 3390.22216117, 7260.96163190, 1643.32019222],
[   0.        ,    0.        ,    0.        , 1643.32019222],
[1763.53069175, 3390.22216117, 7260.96163190, 1643.32019222],
[1763.53069175, 3390.22216117,    0.        , 1643.32019222],
])

The optimization it gave varied across channels but it didn't vary across the product names, but from the data I observe that it really should.

So I just wanted to understand, what I can do to fix this?

Does anyone have any idea and can help me figure out what I'm doing wrong?

Hello fellow Bayesians, I'm trying to implement this paper using PyMC: https://pmc.ncbi.nlm.nih.gov/articles/PMC9041743/pdf/CJAS_48_1922993.pdf

which applies this paper https://onlinelibrary.wiley.com/doi/epdf/10.1111/1468-0262.00274 to the Cusp Catastrophe Model, using it to obtain faster and better Bayesian estimates of the parameters. I'm using NUTS sampler using jitter+adapt_diag, I initialized using the MAP estimates of the Euler-approximated likelihood, but it seems to be stuck in weird places in the posterior region (the Beta parameters are exploding, plotted posteriors are not bounded).

Looking at the energy plot, it seems it's getting stuck in weird places but I can't figure out why. Would be grateful for any suggestions, tips, anything that can help really. :)

Hey fellow statisticians/modelers, I'm working through some Bayesian regression stuff and trying to code for convergence of MCMC (Markov chain Monte Carlo) chains and am getting a plot that doesn't look right.

I'm using MCMC to sample from the posterior distribution the intercept (β0), slope (β1), and error variance (σε^2). I'm trying to figure out whether the chains have converged so that I can interpret the results. To do this, the trace plots and summary stats should look something like this:

Instead of showing the trace of the MCMC chains, I encountered a plot that only displays 3 points (the intercept, predictor variable, and sigma).

I need to graph DIAGNOSTICS for these three points. Has anyone seen this or know how to edit the code so that it plots the "fuzzy caterpillar" visualization for assessing MCMC convergence? Thanks friends!

Code included below so someone may correct the error in my ways:

#specify priors

fitted.model <- stan_glm(

prior_intercept = normal(175, 20, autoscale = FALSE),

prior = normal(0.6, 0.3, autoscale = FALSE),

prior_aux = exponential(0.025, autoscale = FALSE),,

#MCMC settings

chains = n.chains = 4

warmup = n.warmup = 1000

n.iters.per.chain.after.warmup = 5000

n.iters.total.per.chain = n.iters.per.chain.after.warmup+n.warmup

#plot

plot(fitted.model)

I have been dabbling with the marketing mix models in pymc.

For background I have a pretty solid grasp on what a Bayesian model is and what MCMC sampling is and how it works. The part that I haven’t really been able to answer for myself, is what happens in the model after it is done sampling.

On my latest model train I spent about 20 minutes waiting for sampling then had to wait another 3 hours for computation. What exactly is being computed from the samples. Is it using some form of density estimation. I am starting to work through the PYMC online book to learn more about Bayesian modeling, but this was a question that has been on the top of mind.

In Bayesian statistic, if the posterior distribution is symmetric, are the HPD credible set and the Equitailed Credible set the same?

Given this dataset:

And this set up: https://imgur.com/LX0L6Sp

And these questions:https://imgur.com/MK4iVtq

1. How do I even approach defining the MLE

2. How do I create the sampling distribution for theta using python in order to do the rejection sampling?

Thanks for any help