Consider the following, simple model:
import pymc3 as pm import numpy as np with pm.Model() as model: pm.Binomial('scalar_binomial', n=1, p=0.25) pm.Binomial('arrayed_binomial', n=1, p=0.25, shape=10) pm.Binomial( 'arrayed_binomial_alt', n=np.ones(10), p=np.full(10, 0.25), shape=10) trace = pm.sample(draws=600, chains=2, tune=500) summary = pm.summary( trace, var_names=[ 'scalar_binomial', 'arrayed_binomial', 'arrayed_binomial_alt']) summary
The three RVs are identical binomials, except that the first scalar_binomial
is a scalar, the second arrayed_binomial has shape 10 with scalar arguments,
and the third arrayed_binomial_alt has shape 10 with (identical) arrayed
arguments. Note also that nothing in this model is constrained by observations.
The result of sampling is peculiar:
scalar_binomial is as expected. arrayed_binomial sees 7 of the 10 array
elements with no binomial samples, and means of 0.0. arrayed_binomial_alt
sees all 10 array elements with no binomial samples, and means of 0.0.
Is this a bug with pm.Binomial()? But how could a bug on something this basic escape notice? So perhaps I am misunderstanding.
The same model structure with normals instead of binomials produces expected
import pymc3 as pm import numpy as np with pm.Model() as model: pm.Normal('scalar_normal', mu=1, sigma=0.25) pm.Normal('arrayed_normal', mu=1, sigma=0.25, shape=10) pm.Normal( 'arrayed_normal_alt', mu=np.ones(10), sigma=np.full(10, 0.25), shape=10) trace = pm.sample( draws=600, chains=2, tune=500, step=pm.Metropolis()) # Metropolis so the sampling is identical to above summary = pm.summary( trace, var_names=[ 'scalar_normal', 'arrayed_normal', 'arrayed_normal_alt']) summary
Results, as expected: