代码之家  ›  专栏  ›  技术社区  ›  M. Regan

为什么这个分层泊松模型与生成数据中的真实参数不匹配?

  •  3
  • M. Regan  · 技术社区  · 10 年前

    我正在尝试拟合分层泊松回归,以估计每个组和全局的时间延迟。我对pymc是否自动将日志链接函数应用于 mu 或者我必须明确这样做:

    with pm.Model() as model:
        alpha = pm.Gamma('alpha', alpha=1, beta=1)
        beta = pm.Gamma('beta', alpha=1, beta=1)
    
        a = pm.Gamma('a', alpha=alpha, beta=beta, shape=n_participants)
    
        mu = a[participants_idx]
        y_est = pm.Poisson('y_est', mu=mu, observed=messages['time_delay'].values)
    
        start = pm.find_MAP(fmin=scipy.optimize.fmin_powell)
        step = pm.Metropolis(start=start)
        trace = pm.sample(20000, step, start=start, progressbar=True)
    

    以下跟踪图显示了 a 。您可以看到组估计值介于0和750之间。

    traceplot

    当我使用平均值绘制超参数伽马分布时,我开始感到困惑 alpha beta 作为参数。下面的分布显示支持率大约在0到5之间,这不符合我的预期,同时查看组估计 在上面什么是 代表它是 log(a) 还是其他什么?

    enter image description here

    谢谢你的指点。


    按照注释中的要求添加使用假数据的示例:这个示例只有一个组,因此应该更容易看到超参数是否能够合理地生成组的泊松分布。

    test_data = []
    model = []
    
    for i in np.arange(1):
        # between 1 and 100 messages per conversation
        num_messages = np.random.uniform(1, 100)
        avg_delay = np.random.gamma(15, 1)
        for j in np.arange(num_messages):
            delay = np.random.poisson(avg_delay)
    
            test_data.append([i, j, delay, i])
    
        model.append([i, avg_delay])
    
    model_df = pd.DataFrame(model, columns=['conversation_id', 'synthetic_mean_delay'])
    test_df = pd.DataFrame(test_data, columns=['conversation_id', 'message_id', 'time_delay', 'participants_str'])
    test_df.head()
    
    # Estimate parameters of model using test data
    # convert categorical variables to integer
    le = preprocessing.LabelEncoder()
    test_participants_map = le.fit(test_df['participants_str'])
    test_participants_idx = le.fit_transform(test_df['participants_str'])
    n_test_participants = len(test_df['participants_str'].unique())
    
    with pm.Model() as model:
        alpha = pm.Gamma('alpha', alpha=1, beta=1)    
        beta = pm.Gamma('beta', alpha=1, beta=1)
    
        a = pm.Gamma('a', alpha=alpha, beta=beta, shape=n_test_participants)
    
        mu = a[test_participants_idx]
    
        y = test_df['time_delay'].values
        y_est = pm.Poisson('y_est', mu=mu, observed=y)
    
        start = pm.find_MAP(fmin=scipy.optimize.fmin_powell)
        step = pm.Metropolis(start=start)
        trace = pm.sample(20000, step, start=start, progressbar=True)
    

    enter image description here

    我不知道下面的超参数是如何产生参数介于13和17之间的泊松分布的。

    enter image description here

    1 回复  |  直到 10 年前
        1
  •  3
  •   M. Regan    10 年前

    回答:pymc使用不同于scipy的参数来表示Gamma分布。scipy使用alpha&scale,而pymc使用alpha和beta。以下模型按预期工作:

    with pm.Model() as model:
        alpha = pm.Gamma('alpha', alpha=1, beta=1)    
        scale = pm.Gamma('scale', alpha=1, beta=1)
    
        a = pm.Gamma('a', alpha=alpha, beta=1.0/scale, shape=n_test_participants)
    
        #mu = T.exp(a[test_participants_idx])
        mu = a[test_participants_idx]
    
        y = test_df['time_delay'].values
        y_est = pm.Poisson('y_est', mu=mu, observed=y)
    
        start = pm.find_MAP(fmin=scipy.optimize.fmin_powell)
        step = pm.Metropolis(start=start)
        trace = pm.sample(20000, step, start=start, progressbar=True)
    

    enter image description here

    enter image description here

    推荐文章