PyMC 5.10.0中MMM自变量滞后参数化遇错求助

错误信息 #

Cell In[15], line 75
     72 # Creating a tensor to store transformed values
     73 adstocked = tt.zeros_like(x)
---> 75 for i in range(max_lag, x.shape[0]):
     76     weights = tt.power(rate, tt.arange(max_lag + 1))
     77     adstocked = tt.set_subtensor(adstocked[i], tt.dot(x[i-max_lag:i+1][::-1], weights))

TypeError: 'TensorVariable' object cannot be interpreted as an integer

可复现代码 #

## Create a simple MMM data 
import pandas as pd
from random import randint
import numpy as np
import pytensor.tensor as tt
import pytensor as pt
import pymc as pm
import pymc.sampling.jax as pmjax
import arviz as az


# Generate date range
dates = pd.date_range(start="2021-01-01", end="2022-01-01")

data = {
    "date": dates,
    "gcm_direct_Impressions": [randint(10000, 20000) for _ in dates],
    "display_direct_Impressions" :[randint(100000,150000) for _ in dates],
    "tv_grps": [randint(30, 50) for _ in dates],
    "tiktok_direct_Impressions": [randint(10000, 15000) for _ in dates],
    "sell_out_quantity": [randint(150, 250) for _ in dates]
}
df = pd.DataFrame(data)
m = max(df['sell_out_quantity'].values)

print(f"Max sales Volume {m}")

channel_columns = [col for col in df.columns if 'Impressions' in col or 'grps' in col]

transform_variables = channel_columns


delay_channels = channel_columns

media_channels = channel_columns

target = 'sell_out_quantity'

### Transform each channel variable

data_transformed = df.copy()

numerical_encoder_dict = {}


for feature in transform_variables:
    # Extracting the original values of the feature.
    original = df[feature].values

    # Calculating the maximum value of the feature.
    max_value = original.max()

    # Dividing each value in the feature by the maximum value.
    transformed = original / max_value

    # Storing the transformed data back into the 'data_transformed' DataFrame.
    data_transformed[feature] = transformed

    # Storing the maximum value used for scaling in the dictionary.
    # This will be used for reversing the transformation if needed.
    numerical_encoder_dict[feature] = max_value



def adstock_transform(x, rate,max_lag):
    """ Apply adstock transformation with PyTensor.
    :param x: PyTensor tensor, original data for the channel
    :param rate: PyTensor tensor, decay rate of the adstock transformation
    :param max_lag: int, maximum lag to consider for the adstock effect
    :return: PyTensor tensor, transformed data
    """
    # Creating a tensor to store transformed values
    adstocked = tt.zeros_like(x)
    
    for i in range(max_lag, x.shape[0]):
        weights = tt.power(rate, tt.arange(max_lag + 1))
        adstocked = tt.set_subtensor(adstocked[i], tt.dot(x[i-max_lag:i+1][::-1], weights))
    
    return adstocked

### Create a model
response_mean = []

with pm.Model() as model_2:
    # Looping through each channel in the list of delay channels.
    for channel_name in delay_channels:
        print(f"Delay Channels: Adding {channel_name}")

        # Extracting the transformed data for the current channel.
        x = data_transformed[channel_name].values

        # Defining Bayesian priors for the adstock, gamma, and alpha parameters for the current channel.
        adstock_param = pm.Beta(f"{channel_name}_adstock", 2, 2)
        saturation_gamma = pm.Beta(f"{channel_name}_gamma", 2, 2)
        saturation_alpha = pm.Gamma(f"{channel_name}_alpha", 3, 1)
        rate = pm.Beta(f'{channel_name}_rate', alpha=1, beta=1)
        lag = pm.DiscreteUniform(f"{channel_name}_lag",lower=0,upper=17)
        
        transformed_X1 = adstock_transform(x,rate,max_lag=lag)
        transformed_X2 = tt.zeros_like(x)
        for i in range(1,len(x)):
            transformed_X2 = tt.set_subtensor(transformed_X2[i],(transformed_X1[i]**saturation_alpha)/(transformed_X1[i]**saturation_alpha+saturation_gamma**saturation_alpha))
        channel_b = pm.HalfNormal(f"{channel_name}_media_coef", sigma = m)
        response_mean.append(transformed_X2 * channel_b)

    intercept = pm.Normal("intercept",mu = np.mean(data_transformed[target].values), sigma = 3)
    sigma = pm.HalfNormal("sigma", 4)
    likelihood = pm.Normal("outcome", mu = intercept + sum(response_mean), sigma = sigma,
                           observed = data_transformed[target].values)

with model_2:
    trace = pmjax.sample_numpyro_nuts(1000, tune=1000, target_accept=0.95)
    
    trace_summary = az.summary(trace)

已尝试的无效方法 #

• 在adstock_transform函数中将max_lag替换为int(max_lag.eval())
• 将x = data_transformed[channel_name].values替换为x = pm.ConstantData("data_{channel_name}",data_transformed[channel_name].values)

错误原因 #

报错核心问题：

1. lag是PyMC定义的离散随机变量（pm.DiscreteUniform），属于TensorVariable类型，Python原生range()无法识别该类型，必须传入整数。
2. 用Python循环实现adstock变换属于静态代码逻辑，无法融入PyTensor动态计算图——模型定义阶段lag无具体数值，循环次数无法动态调整。

已尝试方法无效原因：

• max_lag.eval()仅能在计算图构建完成后获取值，模型定义阶段lag只是变量，无具体数值，会直接报错。
• pm.ConstantData包装数据是正确方向，但未解决循环依赖离散变量的核心问题。

解决方案 #

1. 重写adstock变换为向量化实现 #

利用PyTensor张量操作替代Python循环，由于lag取值范围是0-17（有限值），预先计算所有可能滞后的adstock结果，再根据采样的lag选择对应值：

def adstock_transform(x, rate, max_lag_possible):
    """向量化实现adstock变换，支持动态选择滞后值"""
    x_tensor = tt.as_tensor_variable(x)
    n_obs = x_tensor.shape[0]
    adstock_results = []
    
    # 预先计算所有可能滞后（0到max_lag_possible）的结果
    for lag in range(max_lag_possible + 1):
        # 补零处理开头的滞后窗口
        padded_x = tt.concatenate([tt.zeros(lag), x_tensor])
        # 生成所有时间步的滞后窗口
        windows = tt.stack([padded_x[i:i+n_obs] for i in range(lag+1)], axis=1)[:, ::-1]
        # 计算滞后权重并归一化
        weights = tt.power(rate, tt.arange(lag+1))
        weights = weights / tt.sum(weights)
        # 计算当前滞后的adstock值
        adstocked = tt.dot(windows, weights)
        adstock_results.append(adstocked)
    
    # 将所有结果堆叠为张量，方便后续通过索引选择
    return tt.stack(adstock_results, axis=0)

2. 修改模型定义 #

使用向量化adstock变换，同时将饱和度变换改为向量化操作：

with pm.Model() as model_2:
    response_mean = []
    max_lag_upper = 17  # 和DiscreteUniform的upper一致
    
    for channel_name in delay_channels:
        print(f"Delay Channels: Adding {channel_name}")
        
        # 用pm.ConstantData包装数据，确保PyTensor能追踪
        x = pm.ConstantData(f"data_{channel_name}", data_transformed[channel_name].values)
        
        # 定义先验
        saturation_gamma = pm.Beta(f"{channel_name}_gamma", 2, 2)
        saturation_alpha = pm.Gamma(f"{channel_name}_alpha", 3, 1)
        rate = pm.Beta(f'{channel_name}_rate', alpha=1, beta=1)
        lag = pm.DiscreteUniform(f"{channel_name}_lag", lower=0, upper=max_lag_upper)
        
        # 获取所有可能滞后的adstock结果，根据lag选择对应值
        adstock_stack = adstock_transform(x, rate, max_lag_upper)
        transformed_X1 = adstock_stack[lag]
        
        # 向量化实现饱和度变换，替代Python循环
        transformed_X2 = (transformed_X1 ** saturation_alpha) / (transformed_X1 ** saturation_alpha + saturation_gamma ** saturation_alpha)
        # 保持原逻辑：第一个元素设为0
        transformed_X2 = tt.set_subtensor(transformed_X2[0], 0.0)
        
        channel_b = pm.HalfNormal(f"{channel_name}_media_coef", sigma=m)
        response_mean.append(transformed_X2 * channel_b)
    
    intercept = pm.Normal("intercept", mu=np.mean(data_transformed[target].values), sigma=3)
    sigma = pm.HalfNormal("sigma", 4)
    likelihood = pm.Normal("outcome", mu=intercept + sum(response_mean), sigma=sigma,
                           observed=data_transformed[target].values)

3. 更换采样器 #

numpyro的pmjax.sample_numpyro_nuts仅支持连续变量的NUTS采样，无法处理离散变量。改用PyMC默认采样器，它会自动对离散变量使用Metropolis采样，连续变量使用NUTS：

with model_2:
    trace = pm.sample(1000, tune=1000, target_accept=0.95, cores=4)
    trace_summary = az.summary(trace)

问题回答 #

1. 需要更换采样器：numpyro的NUTS不支持离散变量，PyMC默认的pm.sample()是处理离散-连续混合模型最简便的选择。
2. 离散均匀分布的实现：你当前使用的pm.DiscreteUniform是正确的实现方式，问题不在于分布本身，而是要通过预先计算所有可能结果再索引选择的方式，适配PyTensor动态计算图的逻辑。

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