2017-02-03 65 views
1

我现在正在学习斯坦,想要实现一个简单的混合模型。斯坦的混合模型 - 矢量化

在参考手册(斯坦参考-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); 
    } 
} 

下一个页描述了外环的向量化是不可能的。但是,我想知道内部循环的parallization仍然是。

所以我尝试了以下型号:

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]之间的区别),也许可以深入了解斯坦。

回答

1

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

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

你想拉高效率的ppt的定义。