1

我现在正在学习 Stan，想实现一个简单的混合模型。

在参考手册（stan-reference-2.14.0）中已经有一个解决方案：

data {
  int<lower=1> K; // number of mixture components
  int<lower=1> N; // number of data points
  real y[N]; // observations
}
parameters {
  simplex[K] theta; // mixing proportions
  real mu[K]; // locations of mixture components
  real<lower=0> sigma[K]; // scales of mixture components
}
model {
  real ps[K]; // temp for log component densities
  sigma ~ cauchy(0, 2.5);
  mu ~ normal(0, 10);
  for (n in 1:N) {
    for (k in 1:K) {
      ps[k] = log(theta[k])
      + normal_lpdf(y[n] | mu[k], sigma[k]);
    }
    target += log_sum_exp(ps);
  }
}

下一页描述了外部循环的矢量化是不可能的。但是，我想知道内部循环的并行化是否仍然存在。

所以我尝试了以下模型：

data {
  int<lower=1> K; // number of mixture components
  int<lower=1> N; // number of data points
  real y[N]; // observations
}

parameters {
  simplex[K] theta; // mixing proportions
  vector[K] mu; // locations of mixture components
  vector<lower=0>[K] sigma; // scales of mixture components
}

model {
  vector[K] ps;//[K]; // temp for log component densities
  vector[K] ppt;
  sigma ~ cauchy(0, 2.5);
  mu ~ normal(0, 10);
  for (n in 1:N) {
    ppt = log(theta);
    /*
    for (k in 1:K) {
      ps[k] = ppt[k] + //log(theta[k])
        normal_lpdf(y[n] | mu[k], sigma[k]);
    }
    */
    ps = ppt + normal_lpdf(y[n] | mu, sigma);
    target += log_sum_exp(ps);
  }
}

...并且该模型做出错误的估计（与原始模型相反）。

data("faithful")
erupdata <- list(
  K = 2,
  N = length(faithful$eruptions),
  y = faithful$eruptions
)

fiteruptions <- stan(file = 'mixturemodel.stan', data = erupdata, iter = 1000, chains = 1)

我想知道，我对模型规范的理解是错误的。我想了解语法提供的区别（以及vector[K]和之间的区别real[K]），也许对 Stan 有更深入的了解。

4

2 回答 2

2

第二个程序定义了不同的密度。 normal_lpdf返回单个标量值，它是数据/参数容器上的日志 pdf 的总和。

手册中有关于矩阵/向量与数组的一章。

您想将ppt更高的定义拉高以提高效率。

于 2017-02-04T18:27:05.783 回答
0

我不是 Stan 用户，我不解释 和 之间的vector[K]区别real[K]。但是根据一些资源（这里和这里），混合位置的先验分布应该被定义为，ordered[K]而不是vector[K]避免不可识别/可交换先验的问题。

为了更深入地理解上述问题，强烈推荐识别贝叶斯混合模型博客。从该博客复制的工作模型：

data {
 int<lower = 0> N;
 vector[N] y;
}

parameters {
  ordered[2] mu;
  real<lower=0> sigma[2];
  real<lower=0, upper=1> theta;
}

model {
 sigma ~ normal(0, 2);
 mu ~ normal(0, 2);
 theta ~ beta(5, 5);
 for (n in 1:N)
   target += log_mix(theta,
                     normal_lpdf(y[n] | mu[1], sigma[1]),
                     normal_lpdf(y[n] | mu[2], sigma[2]));
}
于 2018-12-13T08:35:39.723 回答