使用PyMC進(jìn)行時(shí)間序列分層建模

在統(tǒng)計(jì)建模領(lǐng)域，理解總體趨勢(shì)的同時(shí)解釋群體差異的一個(gè)強(qiáng)大方法是分層(或多層)建模。這種方法允許參數(shù)隨組而變化，并捕獲組內(nèi)和組間的變化。在時(shí)間序列數(shù)據(jù)中，這些特定于組的參數(shù)可以表示不同組隨時(shí)間的不同模式。

今天，我們將深入探討如何使用PyMC(用于概率編程的Python庫(kù))構(gòu)建分層時(shí)間序列模型。

讓我們從為多個(gè)組生成一些人工時(shí)間序列數(shù)據(jù)開始，每個(gè)組都有自己的截距和斜率。

import numpy as np
 import matplotlib.pyplot as plt
 import pymc as pm
 
 # Simulating some data
 np.random.seed(0)
 n_groups = 3  # number of groups
 n_data_points = 100  # number of data points per group
 x = np.tile(np.linspace(0, 10, n_data_points), n_groups)
 group_indicator = np.repeat(np.arange(n_groups), n_data_points)
 slope_true = np.random.normal(0, 1, size=n_groups)
 intercept_true = np.random.normal(2, 1, size=n_groups)
 y = slope_true[group_indicator]*x + intercept_true[group_indicator] + np.random.normal(0, 1, size=n_groups*n_data_points)

我們生成了三個(gè)不同組的時(shí)間序列數(shù)據(jù)。每組都有自己的時(shí)間趨勢(shì)，由唯一的截距和斜率定義。

colors = ['b', 'g', 'r']  # Define different colors for each group
 
 plt.figure(figsize=(10, 5))
 
 # Plot raw data for each group
 for i in range(n_groups):
     plt.plot(x[group_indicator == i], y[group_indicator == i], 'o', color=colors[i], label=f'Group {i+1}')
 
 plt.title('Raw Data with Groups')
 plt.xlabel('Time')
 plt.ylabel('Value')
 plt.legend()
 plt.show()

下一步是構(gòu)建層次模型。我們的模型將具有組特定的截距(alpha)和斜率(beta)。截距和斜率是從具有超參數(shù)mu_alpha、sigma_alpha、mu_beta和sigma_beta的正態(tài)分布中繪制的。這些超參數(shù)分別表示截距和斜率的組水平均值和標(biāo)準(zhǔn)差。

with pm.Model() as hierarchical_model:
     # Hyperpriors
     mu_alpha = pm.Normal('mu_alpha', mu=0, sigma=10)
     sigma_alpha = pm.HalfNormal('sigma_alpha', sigma=10)
     mu_beta = pm.Normal('mu_beta', mu=0, sigma=10)
     sigma_beta = pm.HalfNormal('sigma_beta', sigma=10)
   
     # Priors
     alpha = pm.Normal('alpha', mu=mu_alpha, sigma=sigma_alpha, shape=n_groups)  # group-specific intercepts
     beta = pm.Normal('beta', mu=mu_beta, sigma=sigma_beta, shape=n_groups)  # group-specific slopes
     sigma = pm.HalfNormal('sigma', sigma=1)
 
     # Expected value
     mu = alpha[group_indicator] + beta[group_indicator] * x
 
     # Likelihood
     y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)
 
     # Sampling
     trace = pm.sample(2000, tune=1000)

現(xiàn)在我們已經(jīng)定義了模型并對(duì)其進(jìn)行了采樣。讓我們檢查不同參數(shù)的模型估計(jì):

# Checking the trace
 pm.plot_trace(trace,var_names=['alpha','beta'])
 plt.show()

最后一步是將原始數(shù)據(jù)和模型預(yù)測(cè)可視化:

# Posterior samples
 alpha_samples = trace.posterior['alpha'].values
 beta_samples = trace.posterior['beta'].values
 
 # New x values for predictions
 x_new = np.linspace(0, 10, 200)
 
 plt.figure(figsize=(10, 5))
 
 # Plot raw data and predictions for each group
 for i in range(n_groups):
     # Plot raw data
     
     plt.plot(x[group_indicator == i], y[group_indicator == i], 'o', color=colors[i], label=f'Group {i+1} observed')
     x_new = x[group_indicator == i]
     # Generate and plot predictions
     alpha = trace.posterior.sel(alpha_dim_0=i,beta_dim_0=i)['alpha'].values
     beta = trace.posterior.sel(alpha_dim_0=i,beta_dim_0=i)['beta'].values
     y_hat = alpha[..., None] + beta[..., None] * x_new[None,:]
     y_hat_mean = y_hat.mean(axis=(0, 1))
     y_hat_std = y_hat.std(axis=(0, 1))
     plt.plot(x_new, y_hat_mean, color=colors[i], label=f'Group {i+1} predicted')
     plt.fill_between(x_new, y_hat_mean - 2*y_hat_std, y_hat_mean + 2*y_hat_std, color=colors[i], alpha=0.3)
 
 plt.title('Raw Data with Posterior Predictions by Group')
 plt.xlabel('Time')
 plt.ylabel('Value')
 plt.legend()
 plt.show()

從圖中可以看出，分層時(shí)間序列模型很好地捕獲了每組中的單個(gè)趨勢(shì),而陰影區(qū)域給出了預(yù)測(cè)的不確定性。

層次模型為捕獲時(shí)間序列數(shù)據(jù)中的組級(jí)變化提供了一個(gè)強(qiáng)大的框架。它們?cè)试S我們?cè)诮M之間共享統(tǒng)計(jì)數(shù)據(jù)，提供部分信息池和對(duì)數(shù)據(jù)結(jié)構(gòu)的細(xì)微理解。使用像PyMC這樣的庫(kù)，實(shí)現(xiàn)這些模型變得相當(dāng)簡(jiǎn)單，為健壯且可解釋的時(shí)間序列分析鋪平了道路。