使用PyMC进行时间序列分层建模

 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)


 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()


 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)


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


 # 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()


|
2月前
|

89 3
|
2月前
|

26 1
|
2月前

59 4
|
9月前
|

70 0
|

1153 0
|
9月前
|

245 2
|
10月前
|

388 1
|
10月前
|

R-建模 randomForest

88 0
|
10月前
|

R语言-建模(广义)线性(加性、混合)模型

408 0
|

⭐本专栏针对FPGA进行入门学习，从数电中常见的逻辑代数讲起，结合Verilog HDL语言学习与仿真，主要对组合逻辑电路与时序逻辑电路进行分析与设计，对状态机FSM进行剖析与建模。
88 0