# How to estimate state specific parameters in PyMC3 model
**URL:** https://discourse.pymc.io/t/how-to-estimate-state-specific-parameters-in-pymc3-model/8447
**Category:** Questions
**Created:** [December 13, 2021, 8:36pm UTC](https://discourse.pymc.io/t/how-to-estimate-state-specific-parameters-in-pymc3-model/8447 "2021-12-13T20:36:19Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@michnag](https://discourse.pymc.io/u/michnag)
#### Post date: [December 13, 2021, 8:36pm UTC](https://discourse.pymc.io/t/how-to-estimate-state-specific-parameters-in-pymc3-model/8447/1 "2021-12-13T20:36:19Z")
I have the following bayesian regression
where I want to estimate the parameters beta which depend on the state `s=[1, 2, 3]` which are given by an array `S`. `y_it` is of shape `100x90` and `X_it` is of shape `100x90x20` and `beta` is of shape `20x1` while `N=90` refers to the number of individuals, `T=90` to the time increments and `K=20` to the number of predictors.
For the priors I stick to weakly informative priors such a standard normal distributions.
In absence of the state variable, my model would like as follows:
```auto
with pm.Model() as model:
alpha = pm.Normal(name="alpha", mu=0.0, sigma=1.0, shape=())
beta = pm.Normal(name="beta", mu=0.0, sigma=1.0, shape=K)
eps = pm.Gamma(name="eps", alpha=9.0, beta=4.0, shape=())
y_hat = pm.Deterministic(name="y_hat", var=pm.math.dot(X, beta) + alpha)
y_like = pm.Normal(name="y_like", mu=y_hat, sigma=eps, observed=y)
```
How can I accommodate the state variable in the model at hand?
Thank you very much!
---
### Author:  [@Martin\_Ingram](https://discourse.pymc.io/u/Martin_Ingram)
#### Post date: [December 13, 2021, 9:16pm UTC](https://discourse.pymc.io/t/how-to-estimate-state-specific-parameters-in-pymc3-model/8447/2 "2021-12-13T21:16:10Z")
Hi there,
Just to clarify: do I understand you correctly that the state s\_{it} is known for each observation? Or do you need to estimate which state each observation belongs to, too?
---
### Author:  [@michnag](https://discourse.pymc.io/u/michnag)
#### Post date: [December 14, 2021, 7:34am UTC](https://discourse.pymc.io/t/how-to-estimate-state-specific-parameters-in-pymc3-model/8447/3 "2021-12-14T07:34:20Z")
Hey Martin,
thanks for your reply.
In first place the state is observed for every observation.
However, I would also be keen how the model looks like if one has to estimate the states, too.
---
### Author:  [@Martin\_Ingram](https://discourse.pymc.io/u/Martin_Ingram)
#### Post date: [December 14, 2021, 10:31am UTC](https://discourse.pymc.io/t/how-to-estimate-state-specific-parameters-in-pymc3-model/8447/4 "2021-12-14T10:31:29Z")
Hi again,
I had a go at this. Here’s the code:
```python
# Make up some random data
K = 20
N = 100
T = 90
S = 3
X = np.random.randn(N, T, K)
y = np.random.randn(N, T)
# Number of states
# States, zero-indexed
s = np.random.randint(0, S, size=(N, T))
with pm.Model() as model:
# Note the added dimension for alpha and beta
alpha = pm.Normal(name="alpha", mu=0.0, sigma=1.0, shape=(S))
beta = pm.Normal(name="beta", mu=0.0, sigma=1.0, shape=(S, K))
eps = pm.Gamma(name="eps", alpha=9.0, beta=4.0, shape=())
# This picks out the relevant alphas for each observation
# It should have shape (N, T)
rel_alpha = alpha[s]
# This picks out the relevant betas for each observation
# It should have shape (N, T, K)
rel_beta = beta[s]
# The sum here computes the dot product over the last axis; maybe there is a nicer/faster way.
y_hat = pm.Deterministic(name="y_hat", var=pm.math.sum(X * rel_beta, axis=2) + rel_alpha)
# This is the same as before
y_like = pm.Normal(name="y_like", mu=y_hat, sigma=eps, observed=y)
```
The trick is basically to use the matrix of states `s` to index `alpha` and `beta` and pull out the relevant coefficients. A slight complication is that we can’t seem to use `pm.math.dot` with this because there’s an extra axis now. I got around this using the `pm.math.sum` expression, but there’s probably a nicer way, maybe using `tensordot`. This should work though.
I’ve tried to comment each step; let me know if it makes sense! To work out the shapes, it helped me to play around with things in `numpy` first:
```python
sample_beta = np.random.randn(S, K)
sample_alpha = np.random.randn(S)
sample_beta[s].shape, X.shape, sample_alpha[s].shape
```
This returns:
```auto
((100, 90, 20), (100, 90, 20), (100, 90))
```
If you want to estimate the states, too, things get a bit more complicated. Basically, you’ll want to fit something quite similar to a mixture model, so maybe looking at some tutorials for that (maybe e.g. [here](https://docs.pymc.io/en/stable/pymc-examples/examples/variational_inference/gaussian-mixture-model-advi.html)) could give you some initial pointers.
---
### Author:  [@michnag](https://discourse.pymc.io/u/michnag)
#### Post date: [December 15, 2021, 4:58pm UTC](https://discourse.pymc.io/t/how-to-estimate-state-specific-parameters-in-pymc3-model/8447/5 "2021-12-15T16:58:51Z")
Works fine.
Thanks a lot!