PyMC构建分层贝叶斯时间序列模型时遇形状不兼容ValueError

我正在用Python的PyMC 5.4构建分层贝叶斯时间序列模型，分析四家门店的销售数据。所有门店存在季节性成分，后续计划用傅里叶级数建模，当前先从线性模型入手（注释部分为后续傅里叶级数代码）。

模型变量说明 #

• stores：array([2, 5, 6, 7])，门店编号
• coords：字典，{'fourier_features': array([0,1,...,19]), 'stores': array([2,5,6,7])}
• t：形状为(4,59)的数组，代表4家门店各59个月度时间观测值
• y_sales：形状为(4,59)的数组，代表4家门店各59个月度销售观测值

原模型代码 #

with pm.Model(check_bounds=False, coords=coords) as partial_pooled_linear:
    # Hyperpriors
    mu_alpha = pm.Normal("mu_alpha", mu=0, sigma=0.5)
    mu_beta = pm.Normal("mu_beta", mu=0, sigma=0.5)
    mu_sigma = pm.HalfNormal("mu_sigma", sigma=0.1)
    #β_fourier = pm.Normal("β_fourier", mu=0, sigma=0.01, dims="fourier_features")
    #seasonality = pm.Deterministic("seasonality", pm.math.dot(β_fourier,
    #              fourier_features.to_numpy().T))

    # Store-level priors
    alpha = pm.Normal("alpha", mu=mu_alpha, sigma=1, dims=("stores",))
    beta = pm.Normal("beta", mu=mu_beta, sigma=1, dims=("stores",))
    sigma = pm.HalfNormal("sigma", sigma=mu_sigma, dims=("stores",))
    
    # Define trend for each store
    trend = pm.Deterministic("trend", alpha + (beta * t).T)
    #mu = pm.Deterministic("mu", alpha[stores] + trend, dims=("stores"))

    # Likelihood for each store
    y = pm.Normal("likelihood", mu=trend, observed=y_sales)

运行错误信息 #

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:441, in Elemwise.get_output_info(self, dim_shuffle, *inputs)
    439 try:
    440     out_shapes = [
--> 441         [
    442             get_most_specialized_shape(shape)
    443             for shape in zip(*[inp.type.shape for inp in inputs])
    444         ]
    445     ] * shadow.nout
    446 except ValueError:

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:442, in (.0)
    439 try:
    440     out_shapes = [
    441         [
--> 442             get_most_specialized_shape(shape)
    443             for shape in zip(*[inp.type.shape for inp in inputs])
    444         ]
    445     ] * shadow.nout
    446 except ValueError:

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:434, in Elemwise.get_output_info..get_most_specialized_shape(shapes)
    433 if len(shapes) > 1:
--> 434     raise ValueError
    435 return tuple(shapes)[0]

ValueError: 

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
Cell In[66], line 17
     14 sigma = pm.HalfNormal("sigma", sigma=mu_sigma, dims=("stores",))
     16 # Define trend for each store
--> 17 trend = pm.Deterministic("trend", alpha + (beta * t).T)
     18 #mu = pm.Deterministic("mu", alpha[stores] + trend, dims=("stores"))
     19 
     20 # Likelihood for each store
     21 y = pm.Normal("likelihood", mu=trend, observed=y_sales)

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/var.py:133, in _tensor_py_operators.__mul__(self, other)
    129 def __mul__(self, other):
    130     # See explanation in __add__ for the error caught
    131     # and the return value in that case
    132     try:
--> 133         return at.math.mul(self, other)
    134     except (NotImplementedError, TypeError):
    135         return NotImplemented

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/pytensor/graph/op.py:295, in Op.__call__(self, *inputs, **kwargs)
    253 r"""Construct an `Apply` node using :meth:`Op.make_node` and return its outputs.
    254 
    255 This method is just a wrapper around :meth:`Op.make_node`.
   (...)
    292 
    293 """
    294 return_list = kwargs.pop("return_list", False)
--> 295 node = self.make_node(*inputs, **kwargs)
    297 if config.compute_test_value != "off":
    298     compute_test_value(node)

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:485, in Elemwise.make_node(self, *inputs)
    479 """
    480 If the inputs have different number of dimensions, their shape
    481 is left-completed to the greatest number of dimensions with 1s
    482 using DimShuffle.
    483 """
    484 inputs = [as_tensor_variable(i) for i in inputs]
--> 485 out_dtypes, out_shapes, inputs = self.get_output_info(DimShuffle, *inputs)
    486 outputs = [
    487     TensorType(dtype=dtype, shape=shape)()
    488     for dtype, shape in zip(out_dtypes, out_shapes)
    489 ]
    490 return Apply(self, inputs, outputs)

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:447, in Elemwise.get_output_info(self, dim_shuffle, *inputs)
    440     out_shapes = [
    441         [
    442             get_most_specialized_shape(shape)
    443             for shape in zip(*[inp.type.shape for inp in inputs])
    444         ]
    445     ] * shadow.nout
    446 except ValueError:
--> 447     raise ValueError(
    448         f"Incompatible Elemwise input shapes {[inp.type.shape for inp in inputs]}"
    449     )
    451 # inplace_pattern maps output idx -> input idx
    452 inplace_pattern = self.inplace_pattern

ValueError: Incompatible Elemwise input shapes [(1, 4), (4, 59)]

错误原因分析 #

错误出现在beta * t的运算中：

• beta是维度为(4,)的张量（对应4家门店的斜率参数）
• t是形状为(4,59)的时间序列数组
  两者直接相乘时，PyTensor无法自动广播匹配维度，导致形状不兼容。

修正后的代码 #

核心是将alpha和beta扩展维度，使其能与t进行广播运算：

with pm.Model(check_bounds=False, coords=coords) as partial_pooled_linear:
    # Hyperpriors
    mu_alpha = pm.Normal("mu_alpha", mu=0, sigma=0.5)
    mu_beta = pm.Normal("mu_beta", mu=0, sigma=0.5)
    mu_sigma = pm.HalfNormal("mu_sigma", sigma=0.1)
    #β_fourier = pm.Normal("β_fourier", mu=0, sigma=0.01, dims="fourier_features")
    #seasonality = pm.Deterministic("seasonality", pm.math.dot(β_fourier,
    #              fourier_features.to_numpy().T))

    # Store-level priors
    alpha = pm.Normal("alpha", mu=mu_alpha, sigma=1, dims=("stores",))
    beta = pm.Normal("beta", mu=mu_beta, sigma=1, dims=("stores",))
    sigma = pm.HalfNormal("sigma", sigma=mu_sigma, dims=("stores",))
    
    # Define trend for each store：扩展维度实现广播
    trend = pm.Deterministic("trend", alpha[:, None] + beta[:, None] * t)
    #mu = pm.Deterministic("mu", alpha[stores] + trend, dims=("stores"))

    # Likelihood for each store
    y = pm.Normal("likelihood", mu=trend, observed=y_sales)

修正说明 #

• alpha[:, None]将(4,)的截距数组转换为(4,1)的形状
• beta[:, None]将(4,)的斜率数组转换为(4,1)的形状
• 这样两者都能和(4,59)的t进行广播运算，最终得到(4,59)的trend数组，与y_sales的形状完全匹配，可直接作为似然函数的均值输入。

内容的提问来源于stack exchange，提问作者vbrei