Thanks @fehiepsi

I want to define this hierarchical dependency

b[i,j] \sim N(b_\text{mean}[j], b_\text{scale}[j]) \;\; \forall j \in \{0, 1\} \;\; \forall i \in \{0, 1, 2\}

b_\text{mean}[j] \sim N(0, 1) \;\; \forall j \in \{0, 1\}

b_\text{scale}[j] \sim HN(2) \;\; \forall j \in \{0, 1\}

Currently, I use the following code in my numpyro model:

```
b_mean = numpyro.sample('b_mean', dist.Normal(0, 1), sample_shape=(2,))
b_scale = numpyro.sample('b_scale', dist.HalfNormal(2), sample_shape=(2,))
b = numpyro.sample('b', dist.Normal(b_mean, b_scale), sample_shape=(3,))
```

Here, `b_mean `

, `b_scale`

, and `b`

get the shapes `(3,)`

, `(3,)`

, and `(3, 2)`

which satisfies my requirements.

But it seems from the replies that this is wrong or does not work as intended?

And instead, I should do it like this?

```
b_mean = numpyro.sample('b_mean', dist.Normal(0, 1).expand((2,)))
b_scale = numpyro.sample('b_scale', dist.HalfNormal(2).expand((2,)))
b = numpyro.sample('b', dist.Normal(b_mean, b_scale).expand((3,2)))
```

Could you please confirm if the second code chunk is correct? Also, is there a way to see or test the difference between the two chunk of codes? In my code, I’m getting (exactly) the same estimates for these parameters, though.

@martinjankowiak I will also appreciate your input on this. Thank you both!

Also, I’m not sure what the equivalent code would be in plate notation.