# What is the proper way to do einsum in Pymc3?

**URL:** <https://discourse.pymc.io/t/what-is-the-proper-way-to-do-einsum-in-pymc3/9294>\
**Category:** version agnostic\
**Created:** [April 27, 2022, 6:36pm UTC](https://discourse.pymc.io/t/what-is-the-proper-way-to-do-einsum-in-pymc3/9294 "2022-04-27T18:36:17Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![roesta07](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/roesta07/32/6060_2.png) [@roesta07](https://discourse.pymc.io/u/roesta07)\
**Post date:** [April 27, 2022, 6:36pm UTC](https://discourse.pymc.io/t/what-is-the-proper-way-to-do-einsum-in-pymc3/9294/1 "2022-04-27T18:36:17Z")

</div>

I am trying to convert the following code (written in tensorflow) to Pymc3 which is from the [Book](https://bayesiancomputationbook.com/markdown/chp_06.html); _Bayesian Modeling and Computation in Python, Martin, Kumar & Lao, chapter 6, listing 6.3._

```auto
    y_hat = (intercept[..., None] +
             tf.einsum("ij,...->...i", trend, trend_coeff) +
             tf.einsum("ij,...j->...i", seasonality, seasonality_coeff))

```

I am especially having trouble writing tf.einsum into Pymc3 inside the pymc model …  
However, the full code block is as follows

```auto
tfd = tfp.distributions
root = tfd.JointDistributionCoroutine.Root

@tfd.JointDistributionCoroutine
def ts_regression_model():
    intercept = yield root(tfd.Normal(0., 100., name="intercept"))
    trend_coeff = yield root(tfd.Normal(0., 10., name="trend_coeff"))

    seasonality_coeff = yield root(
        tfd.Sample(tfd.Normal(0., 1.),
                   sample_shape=seasonality.shape[-1],
                   name="seasonality_coeff"))
    noise = yield root(tfd.HalfCauchy(loc=0., scale=5., name="noise_sigma"))
    y_hat = (intercept[..., None] +
             tf.einsum("ij,...->...i", trend, trend_coeff) +
             tf.einsum("ij,...j->...i", seasonality, seasonality_coeff))
    observed = yield tfd.Independent(
        tfd.Normal(y_hat, noise[..., None]),
        reinterpreted_batch_ndims=1,
        name="observed")

```

Also, I have heard Theano’s tensor\_dot can help; but is there a another way around?

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [April 28, 2022, 12:22am UTC](https://discourse.pymc.io/t/what-is-the-proper-way-to-do-einsum-in-pymc3/9294/2 "2022-04-28T00:22:28Z")

</div>

I will let @aloctavodia @RavinKumar or @junpenglao chime in, but in the meantime, [this](https://ajcr.net/Basic-guide-to-einsum/) is a pretty thorough discussion of einsum.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [April 28, 2022, 1:03pm UTC](https://discourse.pymc.io/t/what-is-the-proper-way-to-do-einsum-in-pymc3/9294/3 "2022-04-28T13:03:46Z")

</div>

aesara/pymc does not have `einsum` implemented yet, you can use [tensordot](https://aesara.readthedocs.io/en/latest/library/tensor/basic.html#aesara.tensor.tensordot) instead (might need to transpose the input)

---

<div class="post-metadata">

**Author:** ![roesta07](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/roesta07/32/6060_2.png) [@roesta07](https://discourse.pymc.io/u/roesta07)\
**Post date:** [May 21, 2022, 3:13pm UTC](https://discourse.pymc.io/t/what-is-the-proper-way-to-do-einsum-in-pymc3/9294/4 "2022-05-21T15:13:37Z")

</div>

Ok after standardizing the outputs and features I did something like this; but the posteriors were wrong!!

```auto
with pm.Model() as model_1:
  intercept = pm.Normal("intercept",0,0.2)
  trend_coeff=pm.Normal('trend_coeff',0.0,0.5)
  seasonality_coeff=pm.Normal('seasonality_coeff',0.0,0.5,shape=seasonality.shape[-1])
  noise=pm.Exponential('noise',1)
  y_hat=intercept[...,None] +trend@trend_coeff +tensor.tensordot(seasonality_coeff,seasonality,(0,1))
  y= pm.Normal('y',mu=y_hat,sigma=noise,observed=co2_by_month_training_data.CO2.values)
  trace_1=pm.sample(cores=2)

```

Is there a problem with how I specified the model?