Stan-缺失数据的处理

缺失数据的处理

1 missing

在统计研究中，经常遇到缺失数据。stan能够较好地对缺失数据进行处理，在data和transformed data块中声明的变量均视为已知的，而将parameters块中的视为未知。
则涉及缺失的示例编码如下：

data {
  int<lower=0> N_obs;  //可观测数量
  int<lower=0> N_mis;  //缺失数量
  real y_obs[N_obs];   //已观察的数据y_obs
}
parameters {
  real mu;            //待估计的参数
  real<lower=0> sigma;
  real y_mis[N_mis];  //未观测的（未知）
}
model {
  y_obs ~ normal(mu, sigma); //观测数据模型
  y_mis ~ normal(mu, sigma); //未观测数据模型
  //假设数据均服从正态分布
}

上述将观察到和缺失的数据点的数量分别记为N_obs和N_mis，观察到的数据为y_obs，缺失数据记为y_mis，被估计的参数为$\mu$和$\sigma$以及y_mis（缺失未知）。

2 Partially Known Parameters

在一些情况下，例如多元概率密度函数具有部分可观测的数据和参数时，创建一个关于已知（data)和未知的（parameter）的混合向量是非常有必要的。这将在stan中的transformed parameters块中表述和声明。
以下例子涉及一个二元协方差矩阵，其中方差是已知的，但是协方差是未知的：

data {
  int<lower=0> N; //样本量
  vector[2] y[N];  // ？？？
  real<lower=0> var1;  //方差1   
  real<lower=0> var2;  //方差2
}
//对数据进行转换——transformed data块中声明
transformed data {
  real<lower=0> max_cov = sqrt(var1 * var2); //最大协方差，确定协方差上界
  real<upper=0> min_cov = -max_cov;          //最小协方差，确定协方差下界
}
parameters {
  vector[2] mu;     
  real<lower=min_cov, upper=max_cov> cov; //协方差，确定了上界和下界
}
//对参数进行转换——transformed parameters块中声明
transformed parameters {
  matrix[2, 2] Sigma;  //协方差矩阵，2*2维的
  Sigma[1, 1] = var1;     Sigma[1, 2] = cov;
  Sigma[2, 1] = cov;      Sigma[2, 2] = var2;
}
model {
  y ~ multi_normal(mu, Sigma);//多元正态分布模型
}

其中关于协方差矩阵的声明约束条件，确保协方差矩阵为正的。向量化的形式使之处理较为效率。转换后的参数sigma可以包含在model块中的局部变量，如果其不需要在采样器的输出结果中显示的话。

sliced missing data

如果丢失的数据是某个较大数据结构的一部分，那么通常可以使用索引数组和切片对其进行有效的重新组装。下面是时间序列数据的一个示例，其中只观察到序列中的一些项。stan编译代码如下所示：

data {
  int<lower = 0> N_obs;     //可观测样本量
  int<lower = 0> N_mis;     //缺失样本量
  int<lower = 1, upper = N_obs + N_mis> ii_obs[N_obs]; //确定了上界和下界，ii_obs
  int<lower = 1, upper = N_obs + N_mis> ii_mis[N_mis];
  real y_obs[N_obs]; 
}
transformed data {
  int<lower = 0> N = N_obs + N_mis;
}
parameters {
  real y_mis[N_mis];
  real<lower=0> sigma;
}
transformed parameters {
  real y[N];
  y[ii_obs] = y_obs;
  y[ii_mis] = y_mis;
}
model {
  sigma ~ gamma(1, 1);
  y[1] ~ normal(0, 100);
  y[2:N] ~ normal(y[1:(N - 1)], sigma);
}