stan_user_guide

最新推荐文章于 2024-04-01 12:28:47 发布

柚儿娣

最新推荐文章于 2024-04-01 12:28:47 发布

阅读量351

点赞数

原文链接：https://mc-stan.org/

版权

线性回归模型分析

1.1简单的线性回归模型

1.一元线性回归模型，其中y和x均为n维向量，待估计的参数为： $\alpha$ , $\beta$ , $\sigma$

data {
  int<lower=0> N;
  vector[N] x;
  vector[N] y;
}
parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;
}
model {
  y ~ normal(alpha + beta * x, sigma);
}

2.当x为n*k维矩阵时，待估计的参数 $\beta$ 为k维向量：

data {
  int<lower=0> N;   // number of data items
  int<lower=0> K;   // number of predictors
  matrix[N, K] x;   // predictor matrix
  vector[N] y;      // outcome vector
}
parameters {
  real alpha;           // intercept
  vector[K] beta;       // coefficients for predictors
  real<lower=0> sigma;  // error scale
}
model {
  y ~ normal(x * beta + alpha, sigma);  // likelihood
}

用于计算矩估计，可应用
其中表示矩阵每一行的元素为x[n],

for (n in 1:N)
    y[n] ~ normal(x[n] * beta, sigma); //由于beta为一个列向量，则x[n]*beta为尺度类型为real。

3.设计矩阵X中包含截距项， $\beta$ 包含截距项的值

1.2QR的再参数化

前述为：the linear predictor can be written as η=xβη=xβ, where η is a N-vector of predictions, x is a N×Katrix, and β is a K-vector of coefficients.
设计矩阵的分解：分解为一个正交阵Q和一个上三角矩阵R; $X = Q R$

data{
  int <lower=0> N;   // number of data items
  int<lower=0> K;   // number of predictors
  matrix[N, K] x;   // predictor matrix
  vector[N] y;      // outcome vector
}
transformed data {
  matrix[N, K] Q_ast;
  matrix[K, K] R_ast;
  matrix[K, K] R_ast_inverse;
  // thin and scale the QR decomposition
  Q_ast = qr_thin_Q(x) * sqrt(N - 1);
  R_ast = qr_thin_R(x) / sqrt(N - 1);
  R_ast_inverse = inverse(R_ast);        **//求逆inverse**
}
parameters {
  real alpha;           // intercept
  vector[K] theta;      // coefficients on Q_ast
  real<lower=0> sigma;  // error scale
}
model {
  y ~ normal(Q_ast * theta + alpha, sigma);  // likelihood
}
generated quantities {
  vector[K] beta;
  beta = R_ast_inverse * theta; // coefficients on x
}

下面的图片是三者之间的关系，通过引入中间变量 $\theta$ 来估计 $\beta$ 在这里插入图片描述
其中 $\beta$ 是在genreated quantities中产生的。最终得到的 $\beta$ 才是最终估计的参数。
值得注意的是在线性和广义线性模型中当设计矩阵x的维数大于1时，一般使用此QR重新参数化，其应用效果更好。

1.3系数和标度的先验

回归模型中参数先验的设置

1.4 稳健的噪声模型

一般在线性回归的实际应用中，均假设残差项是服从正态分布的，而在实际研究中并不总是如此，可以采用T分布来修改其以上形式：

data {
  ...
  real<lower=0> nu;
}
...
model {
  y ~ student_t(nu, alpha + beta * x, sigma);
}

自由度常数nu被指定为数据，可在数据中声明。

1.5 logistic和Probit回归模型

对于logistic回归模型:

data {
  int<lower=0> N;  //样本量
  vector[N] x;
  int<lower=0,upper=1> y[N];
}
parameters {
  real alpha;
  real beta;
}
model {
  y ~ bernoulli_logit(alpha + beta * x);
}

其它连接函数也可以用其它形式给定，对于采用正态分布的累积概率函数。

y[n] ~ bernoulli(Phi(alpha + beta * x[n]));    //标准的累积概率分布

y[n] ~ bernoulli(Phi_approx(alpha + beta * x[n]));  //近似的累积概率分布

1.6多元logistic回归模型

对于每一个输出变量 $y_n$ 有K个可能的取值，并给定 $\beta$ 的先验为N(0,0.5)

data {
  int K;  //k个取值
  int N;
  int D;
  int y[N];  //？？？
  matrix[N, D] x;   //设计矩阵
}
parameters {
  matrix[D, K] beta;
}
model {
  matrix[N, K] x_beta = x * beta;

  to_vector(beta) ~ normal(0, 5);

  for (n in 1:N)
    y[n] ~ categorical_logit(x_beta[n]');
}

1.7向量的中心化参数

定义一个参数向量 $\beta$ 其累加和为0,这样的参数向量可以用来识别一个多logit回归参数向量(详细信息请参阅多logit部分)，或者可以用来识别IRT模型中的能力或难度参数(但不能两者都使用)(详细信息请参阅项目响应模型部分)。
自由度为K-1;

1.11 Item -response theory models

缺失的数据声明

data {
  int<lower=1> J;              // number of students
  int<lower=1> K;              // number of questions
  int<lower=1> N;              // number of observations
  int<lower=1,upper=J> jj[N];  // student for observation n
  int<lower=1,upper=K> kk[N];  // question for observation n
  int<lower=0,upper=1> y[N];   // correctness for observation n
}

4Truncated or Censored Data

可以按照各自的概率模型在Stan中对已被截断或检查的测量数据进行编码

4.1截断的分布

Stan中的截断仅限于单变量分布，对于该变量，相应的对数累积分布函数（CDF）和对数互补累积分布（CCDF）函数可用。有关截断的分布，CDF和CCDF的更多信息，请参见有关截断的分布的参考手册部分。
截断数据：积分上限为确定，下限无穷；积分下限确定，积分上限为无穷；积分区间。
截断的数据可以在Stan中建模，若假设截断的点为U=300,截断数据 $y_n$ <300.在Stan中对于可观察到的数据建立截断分布。

data {
  int<lower=0> N;
  real U;   //截断的范围
  real<upper=U> y[N];  //y_n的上限是U
}
parameters {
  real mu;   
  real<lower=0> sigma;
}
model {
  for (n in 1:N)
    y[n] ~ normal(mu, sigma) T[,U]; //积分上限确定的形式；
 for (n in 1:N)
  y[n] ~ normal(mu, sigma) T[L,U];   //积分上下限分别确定的形式，需修改对应的声明进行修改。
       **parameters {
            real<lower=L,upper=U> y[N];**
   for (n in 1:N)
  y[n] ~ normal(mu, sigma) T[L,];   //积分下限确定的形式；需要对以下的声明进行修改
   **parameters {
            real<upper=min(y)> L; // L < y[n]
             real<lower=fmax(L, max(y))> U; // L < U; y[n] < U
}**

该模型将一个上界U声明为数据，并约束y的数据以满足约束条件;这将在采样开始前将数据加载到模型中时进行检查。
该模型在尺度和位置参数上隐式地使用了一个不合适的平面先验;这些可以在模型中使用抽样语句给出先验。

未知的截断点

data {
  int<lower=1> N;
  real y[N];
}
parameters {
  real<upper = min(y)> L;    //在参数部分声明下限
  real<lower = max(y)> U;  //在参数部分声明上限
  real mu;
  real<lower=0> sigma;
}
model {
  L ~ ...;  //给定L的先验
  U ~ ...;  //给定U的先验
  for (n in 1:N)
    y[n] ~ normal(mu, sigma) T[L,U];  //模型，截断上下限
}

4.3 Censored Data(截尾数据）

截尾数据与截断数据有所差异，censored和trunched的区别

censored data 与truncated data的区别在于y在范围以外的个体的x看不看的到，看得到是censor，看不到是truncate
censored data与sample selection的区别在于，censored data虽然看不到某些个体的具体的y，但是其y的范围是知道的（left or right），而sample selection中，y看不到，也不知道其范围。
truncated data 与sample selection都是样本选择问题，区别在于truncated data的样本选择是完全依赖于y的，只要y在某个范围以内就可以看到，否则完全（包括x）看不到；而sample selection是能看到所有的x，只不过有一部分人的y看不到。
估计截尾数据：

data {
  int<lower=0> N_obs;    //已观察到的数据点个数
  int<lower=0> N_cens;  //截尾的数据点个数
  real y_obs[N_obs];       //
  real<lower=max(y_obs)> U;
}
parameters {
  real<lower=U> y_cens[N_cens];
  real mu;
  real<lower=0> sigma;
}
model {
  y_obs ~ normal(mu, sigma);
  y_cens ~ normal(mu, sigma);
}