Skip to content

Instantly share code, notes, and snippets.

@chrishanretty
Created November 16, 2021 19:45
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save chrishanretty/c8705539d312d5f22cbe8eaaad4c8ff7 to your computer and use it in GitHub Desktop.
Save chrishanretty/c8705539d312d5f22cbe8eaaad4c8ff7 to your computer and use it in GitHub Desktop.
Likely faulty Ordered Dirichlet R/Stan code
library(tidyverse)
library(cmdstanr)
library(MCMCpack)
### log PDF is based on eqn. 3 in https://www2.seas.gwu.edu/~dorpjr/Publications/Bookchapter/MarcelDekker2003.pdf
stan_code <-
"
functions {
real orddir_lpdf(vector y, real[] alpha) {
int K = size(alpha);
real alpham1[K];
for (i in 1:K) {
alpham1[i] = alpha[i] - 1.0;
}
real deltas[K];
deltas[1] = y[1];
for (i in 2:K) {
deltas[i] = y[i] - y[i-1];
}
real retval = log(tgamma(sum(alpha))) -
sum(log(tgamma(alpha))) +
sum(log(pow(to_array_1d(deltas), alpham1)));
return(retval);
}
}
data {
int K;
int<lower=1> nobs;
simplex[K] y[nobs];
}
parameters {
real<lower=0> alpha[K];
}
model {
for (i in 1:nobs) {
y[i] ~ orddir(alpha);
}
}
"
writeLines(stan_code, con = "custom.stan")
mod <- cmdstan_model("custom.stan")
nobs <- 10000
truth <- c(2, 5, 10, 100)
y <- rdirichlet(n = nobs, alpha = truth)
y <- t(apply(y, 1, sort, decreasing = FALSE))
stan_data <- list(K = 4,
y = y,
nobs = nobs)
fit <- mod$sample(
data = stan_data,
seed = 123,
chains = 4,
parallel_chains = 4,
refresh = 500
)
s <- fit$summary()
###
alpha <- s %>%
filter(grepl("alpha", variable)) %>%
pull(mean)
alpha <- alpha / sum(alpha)
### Can't recover proportions that match the data generating process
cbind(alpha, truth / sum(truth))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment