Problem with broadcasting variable in new dimensions

Questions v5
1

I have two variables v1 and v2, one with dimension (586, 2, 2, 2), and the other with dimension (586,). I want to feed them into a new variable, with dimension (586, 2, 2, 2), as follows:

v3 = pymc.Normal("v3", mu=v1 + v2, dims=(586, 2, 2, 2),)

Obviously, this has the wrong dimensionality, but I cannot for the life of me figure out how to broadcast the values of v2 in the 3 missing dimensions, so that v3[0, 0, 0, 0] is informed by v1[0, 0, 0, 0] and v2[0,], and so forth. Can someone give a pointer? I haven’t been able to find a good example in the documentation. Thanks!

2

This might be a little hacky, but have you tried doing the following?

v3 = pymc.Normal("v3", mu=v1 + v2[:, np.newaxis, np.newaxis, np.newaxis], dims=(586, 2, 2, 2),)
1 Like
3

Not hacky, this is how to do it. In numpy broadcasting semantics (which pytensor follows), you need to insert size-1 dimensions to signal that it is safe to broadcast. See here for a full rundown.

2 Likes
4

Thanks for the quick responses. At the same time, I’m kind of having trouble understanding which functions to use to perform basic mathematical operations on variables, especially choosing where I can use simply numpy functions versus having to resort to pytensor. Here is an MWE something that has been puzzling me.

import pymc
import numpy as np
m = pymc.Model()
arrival_rate_idx = np.array([1,4,2,2,0,3,4])
departure_rate_idx = np.array([1,3,3,1,0,5,2])
with m:
    arrival_rate_effect = pymc.Dirichlet("arrival_rate_effect", a=np.ones(4))
    departure_rate_effect = pymc.Dirichlet("departure_rate_effect", a=np.ones(4))
    arrival_rate_impact = pymc.Normal("arrival_rate_impact", 0, 1)
    departure_rate_impact = pymc.Normal("departure_rate_impact", 0, 1)
    const = pymc.Normal("const", -5, 1)
    departure_rate_effect_vector = np.array([np.sum(departure_rate_effect[0:i-1]) for i in range(5)])
    arrival_rate_effect_vector = np.array([np.sum(arrival_rate_effect[0:i-1]) for i in range(5)])
    logits = pymc.Deterministic("logits",
                  const[np.newaxis] + 
                  arrival_rate_effect_vector[arrival_rate_idx] * arrival_rate_impact + 
                  departure_rate_effect_vector[departure_rate_idx] * departure_rate_impact                               
                 )

This is a component of a binomial GLM where there are ordered effects based on arrival and departure rates, using a product of a Dirichlet and Normal function. Sometimes, numpy operations work, but in this case evaluating departure_rate_effect_vector returns a TypeError: _tensor_py_operators.sum() got an unexpected keyword argument 'out' . How does one fix this, and is there a rule of thumb for when to use numpy vs. pytensor operations?

5

Never use numpy