2017-02-03 17 views
1

Я изучаю Stan сейчас и хочу реализовать простую модель смеси.Mixture Models in Stan - векторный дизайн

В справочном руководстве (Стан-справочно-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]) и, возможно, получить более глубокое понимание Стэна.

ответ

1

Вторая программа определяет разную плотность. normal_lpdf возвращает одно скалярное значение, которое представляет собой сумму лог-файлов pdf над контейнерами данных/параметров.

В руководстве содержится глава о матрицах/векторах против массивов.

Вы хотите получить определение ppt выше по эффективности.