Заказанная оценка Пробита в Стане - PullRequest
Новая
       112

Заказанная оценка Пробита в Стане

0 голосов
/ 27 ноября 2018

Я пытаюсь воспроизвести упорядоченную модель JAGS пробита в книге Джона Крушке «Выполнение байесовского анализа» (стр. 676) в Stan:

Модель JAGS: 

model {
    for ( i in 1:Ntotal ) {
      y[i] ~ dcat( pr[i,1:nYlevels] )
      pr[i,1] <- pnorm( thresh[1] , mu , 1/sigma^2 )
      for ( k in 2:(nYlevels-1) ) {
        pr[i,k] <- max( 0 ,  pnorm( thresh[ k ] , mu , 1/sigma^2 )
                           - pnorm( thresh[k-1] , mu , 1/sigma^2 ) )
      }
      pr[i,nYlevels] <- 1 - pnorm( thresh[nYlevels-1] , mu , 1/sigma^2 )
    }
    mu ~ dnorm( (1+nYlevels)/2 , 1/(nYlevels)^2 )
    sigma ~ dunif( nYlevels/1000 , nYlevels*10 )
    for ( k in 2:(nYlevels-2) ) {  # 1 and nYlevels-1 are fixed, not stochastic
      thresh[k] ~ dnorm( k+0.5 , 1/2^2 )
    }
  }

Пока чтоУ меня есть следующее, которое работает, но не дает те же результаты, что и в книге.Модель Стэна: 

data{
  int<lower=1> n; // number of obs
  int<lower=3> n_levels; // number of categories

  int y[n]; // outcome var 
}

parameters{
  real mu; // latent mean
  real<lower=0> sigma; // latent sd
  ordered[n_levels] thresh; // thresholds

}

model{
  vector[n_levels] pr[n];

  mu ~ normal( (1+n_levels)/2 , 1/(n_levels)^2 );
  sigma ~ uniform( n_levels/1000 , n_levels*10 );


  for ( k in 2:(n_levels-2) ) // 1 and nYlevels-1 are fixed, not stochastic
    thresh[k] ~ normal( k+0.5 , 1/2^2 );

  for(i in 1:n) {

    pr[i, 1] = normal_cdf(thresh[1], mu, 1/sigma^2);

    for (k in 2:(n_levels-1)) {
      pr[i, k] = max([0.0, normal_cdf(thresh[k], mu, 1/sigma^2) - normal_cdf(thresh[k-1], mu, 1/sigma^2)]);
    }

    pr[i, n_levels] = 1 - normal_cdf(thresh[n_levels - 1], mu, 1/sigma^2);

    y[i] ~ categorical(pr[i, 1:n_levels]);
  }

}

Данные здесь: 

list(n = 100L, n_levels = 7, y = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 6L, 6L, 7L))

Должны восстановить mu 1,0 и sigma 2,5.Вместо этого я получаю mu 3,98 и sigma 1,25.

Я уверен, что делаю что-то не так в модели Стэна, но я очень новичок и не уверенчто делать дальше.Спасибо!

1 Ответ

0 голосов
/ 01 декабря 2018

Обновление: после поиска в Интернете (особая благодарность Конор Гольд ), я пришел к этому решению, которое довольно близко копирует результаты книги.Конечно, любая обратная связь / лучшая факторизация модели принесут приемлемый ответ! 

data {
  real L;                     // Lower fixed thresholds
  real<lower=L> U;            // Upper fixed threshold

  int<lower=2> J;             // Number of outcome levels

  int<lower=0> N;             // Data length

  int<lower=1,upper=J> y[N];  // Ordinal responses
}

transformed data {
  real<lower=0> diff;         // difference between upper and lower fixed thresholds
  int<lower=1> K;             // Number of thresholds

  K = J - 1;
  diff = U - L;
}

parameters {
  simplex[K - 1] thresh_raw;      // raw thresholds
  real mu; // latent mean
  real<lower=0> sigma; // latent sd
}

transformed parameters {
  ordered[K] thresh;     // new thresholds with fixed first and last

  thresh[1] = L;
  thresh[2:K] = L + diff * cumulative_sum(thresh_raw);
  thresh[K] = U; // Overwrite last value to fix it
}

model{
  vector[J] theta;                  // local parameter for ordinal categories

  //priors
  mu ~ normal( (1+J)/2.0 , J );
  sigma ~ uniform( J/1000.0 , J*10 );

  for (i in 2:K-2)
    thresh[i] ~ normal(i + .05, 2);

  // likelihood statement
  for(n in 1:N) {

    // probability of ordinal category definitions
    theta[1] = normal_cdf( thresh[1] , mu, sigma );

    for (l in 2:K)
      theta[l] = fmax(0, normal_cdf(thresh[l], mu, sigma ) - normal_cdf(thresh[l-1], mu, sigma));

    theta[J] = 1 - normal_cdf(thresh[K] , mu, sigma);

    y[n] ~ categorical(theta);
  }
}
...