Issue with observed data in joint distribution

Doing p_dow = pm.Multinomial('P_dow',num_dow, prior_dow, shape=num_dow) should suppress the error, however this is a bit odd as I remember you dont need to explicitly write down the shape.

Also, there are two problem of

pm.Potential('mylike', p_sd * p_ed * p_length * p_dow * p_hr, observed=data)
  1. You cannot assign observed to Potential. If your potential function evaluate on some input data it is better to write it as a DensityDist
  2. p_sd, p_ed, p_length, p_dow, p_hr all has different shape - it will get broadcasted in a weird way. are you sure you want to do this?

Also, you should do

trace = pm.sample(1000, tune=1000)

instead of

    step = pm.NUTS()
    trace = pm.sample(20000, step=step)

as passing step directly to sample is not recommanded, and 20000 is likely too many samples.