Hi all,

I’ve been trying to create a hierarchical model using 3 dimensional data from an xarray.

I’m having trouble connected the final level to the likelihood.

In normal cases I would just have done

`mu = exercise_intercept[exercise_idx] + exercise_slope[exercise_idx]*velocity`

.

This is however giving me an error message and I’m lost on how to solve it:

`Input dimension mis-match. (input[0].shape[2] = 5, input[1].shape[2] = 15)`

.

The code for generating the test data and model is as follows:

```
import pymc3 as pm
import numpy as np
import xarray as xr
import pandas as pd
# Generate test data
n_exercises = 5
n_dates = 10
n_observations = 15
exercises = [f'Exercise {i + 1}' for i in np.arange(5)]
dates = pd.date_range(start = '2020-01-01', periods = n_dates)
observations = np.arange(n_observations)
velocity_test_data = np.arange(0.15 + 0.05*n_observations, 0.15, -0.05)
velocities = np.tile(velocity_test_data, n_exercises*n_dates)
loads = list(map(lambda x: 100 + -50*x + np.random.rand()*10, velocities))
matrix_shape = (n_exercises, n_dates, n_observations)
velocity_matrix = np.reshape(velocities, matrix_shape)
load_matrix = np.reshape(loads, matrix_shape)
# Create xarray
dims = ['exercise', 'date', 'observation']
coords = {'exercise': exercises,
'date': dates,
'observation': observations}
dataset = xr.Dataset({'velocity': (dims, velocity_matrix),
'load': (dims, load_matrix)},
coords = coords)
# Create PyMC model
with pm.Model(coords = coords) as model:
# Inputs
exercise_idx = pm.Data('exercise_idx', np.arange(n_exercises), dims = 'exercise')
date_idx = pm.Data('date_idx', np.arange(n_dates), dims = 'date')
obs_idx = pm.Data('obs_idx', np.arange(n_observations), dims = 'observation')
velocity = pm.Data('velocity', dataset['velocity'], dims = dims)
load = pm.Data('load', dataset['load'], dims = dims)
# Global parameters
global_intercept = pm.Normal('global_intercept', mu = 100, sigma = 10)
global_intercept_sd = pm.Exponential('global_intercept_sd', lam = 1)
global_slope = pm.Normal('global_slope', mu = -50, sigma = 10)
global_slope_sd = pm.Exponential('global_slope_sd', lam = 1)
# Exercise parameters
exercise_intercept = pm.Normal('exercise_intercept', mu = global_intercept, sigma = global_intercept_sd, dims = 'exercise')
exercise_slope = pm.Normal('exercise_slope', mu = global_slope, sigma = global_slope_sd, dims = 'exercise')
# Model error
sigma = pm.Exponential('sigma', lam = 1)
# Expected value
mu = exercise_intercept[exercise_idx] + exercise_slope[exercise_idx]*velocity
# Likelihood of observed values
likelihood = pm.Normal('likelihood',
mu = mu,
sigma = sigma,
observed = load,
dims = dims)
```

@OriolAbril perhaps you know how to do it? :