PyStan3预测结果提取与可视化问题求助

在科研项目中，我已将R-Stan脚本转换为Python PyStan3代码，独立完成了Stan模型编写、PyStan调用及演示数据准备，模型可正常运行，但存在三个问题：

• 不知如何绘制观测变量与预测变量的XY图表
• 无法根据给定的x值提取对应的预测变量
• 将fit对象转为DataFrame后，对各列含义模糊，无法区分因变量、自变量/预测变量

import stan
import numpy as np

model_code_2s = '''
        data {
          int<lower=0> N; // number of measurements
        
          vector[N] h;
          vector[N] Q_meas;
          vector<lower=0>[N] Q_error;
          
          real mu_a1;
          real sigma_a1;
          real mu_b1;
          real sigma_b1;
          real mu_c1;
          real sigma_c1;
          real mu_k1;
          real sigma_k1;
          real mu_a2;
          real sigma_a2;
          real mu_c2;
          real sigma_c2;
          
          int<lower=0> N_ext;
          vector[N_ext] h_ext;
        }
        
        parameters {
          real<lower=0> a1;
          real<upper=min(h)> b1;
          real<lower=0> c1;
          
          real<lower=0> a2;
          real<lower=0> c2;
          
          real<lower=b1> k1;
          
          real<lower=0> sigma;
        }
        
        transformed parameters {
        
          real b2 = k1 - pow(a1/a2,1/c2)*pow(k1-b1,c1/c2);
          
          vector[N] log_mu;
          
          for (i in 1:N) {
            if (h[i] <= k1) {
              log_mu[i] = log(a1) + c1*log(h[i]-b1);
            } else {
              log_mu[i] = log(a2) + c2*log(h[i]-b2);
            }
          }
          
          vector[N] mu = exp(log_mu);
          vector[N] sigma_mu = mu*sigma;
          
          vector[N_ext] mu_ext;
          
          for (i in 1:N_ext) {
            if (h_ext[i] <= k1) {
              mu_ext[i] = a1*(h_ext[i] - b1)^c1;
            } else {
              mu_ext[i] = a2*(h_ext[i] - b2)^c2;
            }
          }
        
        }
        
        model {
        
          // priors
          a1 ~ normal(mu_a1, sigma_a1);
          b1 ~ normal(mu_b1, sigma_b1);
          c1 ~ normal(mu_c1, sigma_c1);
          
          k1 ~ normal(mu_k1, sigma_k1);
          
          a2 ~ normal(mu_a2, sigma_a2);
          c2 ~ normal(mu_c2, sigma_c2);
          
          sigma ~ uniform(0,10000);
          
          // likelihood
          Q_meas ~ normal(mu, ((sigma_mu)^2 + Q_error^2)^0.5);
        
        }
        
        generated quantities {
          
          vector[N] Q_res; // residuals Q_res = Q_meas - Q_mu (accounting for uncertainty in Q_meas)
        
          for (i in 1:N) {
            Q_res[i] = normal_rng(Q_meas[i],Q_error[i]) - normal_rng(mu[i],sigma_mu[i]);
          }
          
          vector[N_ext] Q_ext;
          
          for (i in 1:N_ext) {
            if (h_ext[i] <= k1) {
              Q_ext[i] = normal_rng(a1*(h_ext[i]-b1)^c1,mu_ext[i]*sigma);
            } else {
              Q_ext[i] = normal_rng(a2*(h_ext[i]-b2)^c2,mu_ext[i]*sigma);
            }
          }
          
          
        }
'''
h = [3600.0, 2000.0, 1400.0, 1300.0]
Q_meas = [1110.2, 194.4, 90.6, 69.1]
Q_error = [61.3, 8.7, 4.8, 3.7]

N = len(h)
h_ext = np.arange(1300.0, 3600.0, 10) # min(h) to max(h) with 0.1 step
N_ext = len(h_ext)

data_to_fit = {
        'N': N,
        'h': h,
        'Q_meas': Q_meas,
        'Q_error': Q_error,
        'mu_a1': 20, 
        'sigma_a1': 10,
        'mu_b1': 0,
        'sigma_b1': 5,
        'mu_c1': 2,
        'sigma_c1': 0.5,
        'mu_k1': 4,
        'sigma_k1': 2,
        'mu_a2': 20,
        'sigma_a2': 10,
        'mu_c2': 2,
        'sigma_c2': 0.5,      
        'N_ext': N_ext,
        'h_ext': h_ext
        }

posterior = stan.build(model_code_2s, data = data_to_fit, random_seed = 42)
fit = posterior.sample(num_chains = 4, num_samples = 5000, num_warmup = 1000)

1. 绘制观测变量与预测变量的XY图表 #

使用matplotlib绘制观测值、拟合值及扩展预测区间，代码如下：

import matplotlib.pyplot as plt

# 提取观测点拟合均值的统计量
mu_mean = fit['mu'].mean(axis=0)
mu_low = np.percentile(fit['mu'], 2.5, axis=0)
mu_high = np.percentile(fit['mu'], 97.5, axis=0)

# 提取扩展预测的统计量
q_ext_mean = fit['Q_ext'].mean(axis=0)
q_ext_low = np.percentile(fit['Q_ext'], 2.5, axis=0)
q_ext_high = np.percentile(fit['Q_ext'], 97.5, axis=0)

plt.figure(figsize=(10,6))
# 绘制原始观测点
plt.scatter(h, Q_meas, color='red', label='观测值')
# 绘制观测点的拟合置信区间与均值线
plt.fill_between(h, mu_low, mu_high, alpha=0.2, color='blue')
plt.plot(h, mu_mean, color='blue', label='观测点拟合均值')
# 绘制扩展预测的置信区间与均值线
plt.fill_between(h_ext, q_ext_low, q_ext_high, alpha=0.2, color='green')
plt.plot(h_ext, q_ext_mean, color='green', label='扩展预测均值')

plt.xlabel('h（自变量）')
plt.ylabel('Q（因变量）')
plt.legend()
plt.show()

2. 根据给定x值提取对应预测变量 #

通过匹配目标x值在h_ext中的位置，提取对应预测值的采样数据：

# 设定目标h值
h_target = 2500.0

# 找到h_ext中最接近目标值的索引
target_idx = np.argmin(np.abs(h_ext - h_target))

# 提取该位置的所有预测采样值
q_target_samples = fit['Q_ext'][:, target_idx]
# 计算均值与95%置信区间
q_target_mean = q_target_samples.mean()
q_target_ci = np.percentile(q_target_samples, [2.5, 97.5])

print(f"h={h_target}时，Q的预测均值：{q_target_mean:.2f}，95%置信区间：{q_target_ci}")

3. 理解fit转DataFrame后的列含义 #

先将fit对象转为DataFrame，再按Stan代码块区分列的含义：

import pandas as pd

df = fit.to_frame()
# 查看所有列名
print(df.columns.tolist())

列名对应规则：

• 模型参数列：对应Stan代码parameters块中的变量，如a1、b1、k1、sigma等
• 转换参数列：对应transformed parameters块中的变量，如b2、mu[i]（第i个观测点的拟合均值）、mu_ext[j]（第j个扩展h值的预测均值）
• 生成量列：对应generated quantities块中的变量，如Q_res[i]（第i个观测点的残差）、Q_ext[j]（第j个扩展h值的预测采样值）

变量类型区分：

• 自变量：你输入的h、h_ext（不属于fit的列，是外部数据）
• 因变量观测值：你输入的Q_meas（外部数据）
• 因变量预测值：mu（观测点拟合均值）、Q_ext（扩展点预测采样）
• 模型参数：a1、b1、c1等用于构建模型的参数

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