# Combining two components in a single model, with hierarchy

**URL:** <https://discourse.pymc.io/t/combining-two-components-in-a-single-model-with-hierarchy/11154>\
**Category:** v5\
**Tags:** modeling, hierarchical\
**Created:** [January 10, 2023, 7:17pm UTC](https://discourse.pymc.io/t/combining-two-components-in-a-single-model-with-hierarchy/11154 "2023-01-10T19:17:13Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![davipatti](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/davipatti/32/5146_2.png) [@davipatti](https://discourse.pymc.io/u/davipatti)\
**Post date:** [January 10, 2023, 7:17pm UTC](https://discourse.pymc.io/t/combining-two-components-in-a-single-model-with-hierarchy/11154/1 "2023-01-10T19:17:13Z")

</div>

Hello, I have a model with two distinct components that I’m having trouble knitting together.

* * *

**Component 1** is a measurement model that infers a value of interest from raw data:

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/c/c221459e9feb842d161afbd61e670a8d7fb4ce03.png)

y\_i - data is normally distributed with mean coming from a function of a\_i, predictors x\_i, and other (scalar) parameters (...).  
a\_i - value of interest; inferred as part of a basic hierarchical model, with hyperpriors \mu\_a and \sigma\_a.

* * *

**Component 2** is a process model which tries to explain how values of a are generated:

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/8/844f118b48d9175fd574a6654cb018c7b9a31df0.png)

- a\_t is related to its previous value in time (a\_{t-1}) and a function of e\_t as well as other (scalar) parameters (...).

* * *

Initially I had set this up as two separate models:

```python
# component 1 (simplified code)
with pm.Model() as model1:

    # hyperpriors
    mu_a = pm.Normal("mu_a", 0, 1)
    sigma_a = pm.Exponential("sigma_a", 1)

    # priors
    a = pm.Normal("a", mu_a, sigma_a, shape=shape)
    sigma = pm.Exponential("sigma", 1)

    # other parameters
    ...

    # likelihood
    pm.Normal("lik1", f(x, a, ...), sigma=sigma, observed=observed)

    trace1 = pm.sample()

```

Then, I would compute the posterior mean of a:

```python
a_observed = trace1.posterior["a"].mean(["chain", "draw"])

```

and use that as the observed in the second component:

```python
# component 2 - simplified code
with pm.Model() as model2:

    # priors
    p = pm.Beta("p", 1, 1)
    e = pm.Bernoulli("e", p, shape=shape_e)

    # other parameters
    ....

    a = aesara.scan(
        fun,
        sequences=e,
        outputs_info=at.zeros(n),
        non_sequences=[...]
    )
    sigma = pm.Exponential("sigma", 1)

    # likelihood
    pm.Normal("lik2", a, sigma=sigma, observed=a_observed)
    trace2 = pm.sample()

```

This works fine, but ideally uncertainty in a should flow from component 1 to component 2, rather than simply taking the mean (distributions of a are variable, sometimes they are tight and symmetrical, sometimes they are wide and skewed).

However, I don’t see a way of combining these two components without losing the hierarchical inference of a.

E.g. it could be something like:

```python
# combined
with pm.Model() as model3:

   # (from component 2)
    p = pm.Beta("p", 1, 1)
    e = pm.Bernoulli("e", p, shape=shape_e)
    a = aesara.scan(
        fun,
        sequences=e,
        outputs_info=at.zeros(n),
        non_sequences=[...]
    )

    # (from component 1)
    mu_a = pm.Normal("mu_a", 0, 1) # mu_a is now superfluous!
    sigma_a = pm.Exponential("sigma_a", 1)

    # priors
    a = pm.Normal("a", a, sigma_a)
    sigma = pm.Exponential("sigma", 1)
    pm.Normal("lik1", f(x, a, ...), sigma=sigma, observed=observed)

```

But, `mu_a` has become superfluous - we’ve lost the hierarchical sampling of a.

Even if this is the ‘correct’ way to do this, pragmatically the chains don’t mix at all - traces contain the same value repeated.

* * *

I’ve tried various different things (e.g. different formulations, compound steps in sampling, using `pm.Potential`, avoiding using `aesara.scan`) but none that sample.

Does anyone have any suggestions? Perhaps I’m thinking about this problem the wrong way?

Thanks in advance  
David