Broadcasting mathematical operations with PYMC3/Theano
<p>I think this issue boils down to my lack of understanding with <code>Theano</code> works. I'm in a situation where I want to create a variable that is the result of a subtraction between a distribution and a numpy array. This works fine when I specify the shape parameter as <code>1</code></p>

<pre><code>import pymc3 as pm
import numpy as np
import theano.tensor as T

X = np.random.randint(low = -10, high = 10, size = 100)

with pm.Model() as model:
    nl = pm.Normal('nl', shape = 1)
    det = pm.Deterministic('det', nl - x)

nl.dshape
(1,)
</code></pre>

<p>However, this breaks when I specify shape > 1</p>

<pre><code>with pm.Model() as model:
    nl = pm.Normal('nl', shape = 2)
    det = pm.Deterministic('det', nl - X)

ValueError: Input dimension mis-match. (input[0].shape[0] = 2, input[1].shape[0] = 100)

nl.dshape
(2,)

X.shape
(100,)
</code></pre>

<p>I tried transposing X to make it broadcastable</p>

<pre><code>X2 = X.reshape(-1, 1).transpose()

X2.shape
(1, 100)
</code></pre>

<p>But now it declares a mismatch at <code>.shape[1]</code> instead of <code>.shape[0]</code></p>

<pre><code>with pm.Model() as model:
    nl = pm.Normal('nl', shape = 2)
    det = pm.Deterministic('det', nl - X2)

ValueError: Input dimension mis-match. (input[0].shape[1] = 2, input[1].shape[1] = 100)
</code></pre>

<p>I can make this work if I loop over the elements of the distribution</p>

<pre><code>distShape = 2
with pm.Model() as model:
    nl = pm.Normal('nl', shape = distShape)

    det = {}
    for i in range(distShape):
        det[i] = pm.Deterministic('det' + str(i), nl[i] - X)

det
{0: det0, 1: det1}
</code></pre>

<p>However this feels inelegant and constrains me to use loops for the rest of the model. I was wondering if there was a way to specify this operation so that it could work the same as with distributions.</p>

<pre><code>distShape = 2
with pm.Model() as model:
    nl0 = pm.Normal('nl1', shape = distShape)
    nl1 = pm.Normal('nl2', shape = 1)

    det = pm.Deterministic('det', nl0 - nl1)
</code></pre>

<p>You can do</p>

<pre><code>X = np.random.randint(low = -10, high = 10, size = 100)
X = x[:,None] # or x.reshape(-1, 1)
</code></pre>

<p>and then</p>

<pre><code>with pm.Model() as model:
    nl = pm.Normal('nl', shape = 2)
    det = pm.Deterministic('det', nl - X)
</code></pre>

<p>In this case the shapes of nl and X will be ((2, 1), (100,)), respectively and then broadcastable.</p>

<p>Notice we get the same behavior with two NumPy arrays (not only one Theano tensor and one NumPy array)</p>

<pre><code>a0 = np.array([1,2])
b0 = np.array([1,2,3,5])
a0 = a0[:,None]  # comment/uncomment this line
print(a0.shape, b0.shape)
b0-a0
</code></pre>

python theano bayesian pymc3