BIDA: Density Estimation (Finite Mixture Model)
Kazuki Yoshida
2018-11-30
References
- Books
- Web sites
- jSoftware
Background
Here we have 82 data points in the galaxy from a 1-dimensional unknown distribution. The aim is to fit a normal finite mixture model with a pre-specified number of latent clusters.
Assessment thus-far
Each chain seems to converge to something, but the mixing of several chains are very poor. Is there still a label switching issue although the means are ordered in the prior?
Load packages
library(tidyverse)
library(rstan)
## devtools::install_github('jburos/biostan', build_vignettes = TRUE, dependencies = TRUE)
## library(biostan)
library(DPpackage)
set.seed(732565397)Prepare data
data(galaxy, package = "DPpackage")
galaxy <- galaxy %>%
as_data_frame() %>%
mutate(log_speed = log(speed),
k_speed = speed / 1000)
galaxy## # A tibble: 82 x 3
## speed log_speed k_speed
## <dbl> <dbl> <dbl>
## 1 9172 9.12 9.17
## 2 9350 9.14 9.35
## 3 9483 9.16 9.48
## 4 9558 9.17 9.56
## 5 9775 9.19 9.78
## 6 10227 9.23 10.2
## 7 10406 9.25 10.4
## 8 16084 9.69 16.1
## 9 16170 9.69 16.2
## 10 18419 9.82 18.4
## # ... with 72 more rowsggplot(data = galaxy, mapping = aes(x = k_speed)) +
geom_point(y = 0.5) +
scale_y_continuous(limits = c(0,1), breaks = NULL) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
legend.key = element_blank(),
plot.title = element_text(hjust = 0.5),
strip.background = element_blank())
The log speed data seem to show some distinct clusters.
Single Normal
We will start with density estimation with a single normal distribution as a starting point.
Model
\[\begin{align*} y_{i} | (\mu, \sigma^{2}) &\sim N(\mu, \sigma^{2})\\ \\ p(\mu,\sigma^{2}) &= p(\mu|\sigma^{2})p(\sigma^{2})\\ &= N(\mu | m, s^{2}) Inverse-Gamma(\sigma^{2} | \alpha, \beta)\\ &= \left[ \frac{1}{\sqrt{2\pi s^{2}}} \exp \left( - \frac{(\mu - m)^{2}}{2 \times s^{2}} \right) \right] \left[ \frac{\beta^{\alpha}}{\Gamma(\alpha)} (\sigma^{2})^{-\alpha - 1} \exp \left( - \frac{\beta}{\sigma^{2}} \right) \right]\\ &\text{where}\\ m &= 0\\ s^{2} &= 10^{3}\\ \alpha &= 10^{-3}\\ \beta &= 10^{-3}\\ \end{align*} \]
Implementation
normal_stan_code <- biostan::read_stan_file("./bayesianideas_density_normal.stan")## Warning: replacing previous import 'rstan::loo' by 'rstanarm::loo' when loading 'biostan'biostan::print_stan_code(normal_stan_code, section = NULL)## data {
## // Hyperparameters
## real alpha;
## real beta;
## real m;
## real s_squared;
##
## // Define variables in data
## // Number of observations (an integer)
## int<lower=0> n;
## // Outcome (a real vector of length n)
## real y[n];
## }
##
## transformed data {
## real s;
## s = sqrt(s_squared);
## }
##
## parameters {
## // Define parameters to estimate
## // Population mean (a real number)
## real mu;
## // Population variance (a positive real number)
## real<lower=0> sigma_squared;
## }
##
## transformed parameters {
## // Population standard deviation (a positive real number)
## real<lower=0> sigma;
## // Standard deviation (derived from variance)
## sigma = sqrt(sigma_squared);
## }
##
## model {
## // Prior part of Bayesian inference
## // Mean
## mu ~ normal(m, s);
## // sigma^2 has inverse gamma (alpha = 1, beta = 1) prior
## sigma_squared ~ inv_gamma(alpha, beta);
##
## // Likelihood part of Bayesian inference
## // Outcome model N(mu, sigma^2) (use SD rather than Var)
## y ~ normal(mu, sigma);
## }
## Fitting
normal_stan_fit <-
rstan::stan(model_code = normal_stan_code,
data = list(alpha = 10^(-3), beta = 10^(-3),
m = 0, s_squared = 10^(3),
n = nrow(galaxy),
y = galaxy$k_speed))normal_stan_fit## Inference for Stan model: 1d920a28545f72dd74d51f633b0d128b.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 20.83 0.01 0.51 19.84 20.49 20.82 21.17 21.85 3687 1
## sigma_squared 21.40 0.06 3.38 15.81 19.01 21.05 23.39 29.01 3190 1
## sigma 4.61 0.01 0.36 3.98 4.36 4.59 4.84 5.39 3274 1
## lp__ -166.28 0.02 0.99 -168.91 -166.68 -165.98 -165.56 -165.31 1664 1
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:35:27 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).The sampling process seems ok.
normal_stan_draws <- tidybayes::tidy_draws(normal_stan_fit)
## Evaluation grid
eval_grid <- seq(from = floor(min(galaxy$k_speed)), to = ceiling(max(galaxy$k_speed)), length.out = 10^2)
eval_df <- data_frame(x = eval_grid)
## Evaluate each density function over a grid
normal_stan_draws <- normal_stan_draws %>%
## Density function at each parameter draw
mutate(density_fun = pmap(list(mu, sigma),
.f = function(mu, sigma) {
partial(dnorm, mean = mu, sd = sigma)
}),
data = list(eval_df)) %>%
## Grid evaluation
mutate(data = pmap(list(density_fun, data),
.f = function(density_fun, data) {
data %>%
mutate(density = density_fun(x = x))
})) %>%
select(.chain, .iteration, .draw, density_fun, data)Plot density function draws.
normal_stan_draws_unnest <- normal_stan_draws %>%
select(-density_fun) %>%
unnest()
normal_stan_draws_unnest_mean <-
normal_stan_draws_unnest %>%
group_by(x) %>%
summarise(density = mean(density))
normal_stan_draws_unnest %>%
ggplot(mapping = aes(x = x, y = density, group = interaction(.chain, .iteration, .draw))) +
geom_line(size = 0.1, alpha = 1/20) +
geom_line(data = normal_stan_draws_unnest_mean,
color = "red", alpha = 0.5, size = 2,
mapping = aes(group = NULL)) +
geom_point(data = galaxy,
mapping = aes(x = k_speed, group = NULL),
y = 0) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
legend.key = element_blank(),
plot.title = element_text(hjust = 0.5),
strip.background = element_blank())
Apparently, a single normal distribution model is not a sufficient description of the density.
Finite Mixture of Normals
Now we examine if modeling the distribution as a mixture of several underlying cluster-specific normals better fit the data. Here we assume some fixed number of clusters \(H\).
Model
Let \(z_{i} \in \left\{ 1, \dots, H \right\}\) be the cluster membership latent variable and \(\mathbf{z}_{i} = (I(z_{i}=1),\dots,I(z_{i}=H))^{T}\) be the vector indicator version. We have the following model.
\[\begin{align*} y_{i} | (z_{i}, \mu_{1},\dots,\mu_{H}, \sigma^{2}_{1},\dots,\sigma^{2}_{H}) &\sim N(\mu_{z_{i}}, \sigma^{2}_{z_{i}})\\ \mathbf{z}_{i} | \boldsymbol{\pi} &\sim Multinomial(1, \boldsymbol{\pi})\\ \\ p(\mu_{h},\sigma^{2}_{h}) &= p(\mu_{h}|\sigma^{2}_{h})p(\sigma^{2}_{h})\\ &= N(\mu_{h} | m, s^{2}) Inverse-Gamma(\sigma^{2}_{h} | \alpha, \beta)\\ &= \left[ \frac{1}{\sqrt{2\pi s^{2}}} \exp \left( - \frac{(\mu_{h} - m)^{2}}{2 \times s^{2}} \right) \right] \left[ \frac{\beta^{\alpha}}{\Gamma(\alpha)} (\sigma^{2}_{h})^{-\alpha - 1} \exp \left( - \frac{\beta}{\sigma^{2}_{h}} \right) \right]\\ &\text{where}\\ m &= 0\\ s^{2} &= 10^{3}\\ \alpha &= 10^{-3}\\ \beta &= 10^{-3}\\ \\ p(\boldsymbol{\pi}) &\sim Dirichlet \left( \frac{\alpha}{H}, \dots, \frac{\alpha}{H} \right)\\ \end{align*} \]
Summing out the latent categorical variable \(z_{i}\) results in the following (conditioning on parameters suppressed for simplicity) marginal density. Note the cluster membership latent variable \(z_i\) is not measured, thus, it is a discrete parameter. Stan’s Hamiltonian Monte Carlo (HMC) cannot deal with discrete parameters, this marginalization step is required fro Stan. JAGS seems to allow a discrete parameter and accepts the original model above.
\[\begin{align*} p(y) &= \sum^{H}_{z=1} p(y|z) p(z)\\ &= \sum^{H}_{z=1} \pi_{z} p(y|z)\\ &= \sum^{H}_{h=1} \pi_h N(y | \mu_{h}, \sigma^{2}_{h})\\ \end{align*} \]
Another layer of complexity is the label switching issue. That is, the cluster IDs \(\left\{ 1, \dots, H \right\}\) do not specify which cluster corresponds to which ID when all the cluster-specific prior for the normal distribution \(p(\mu_{h},\sigma^{2}_{h})\) is the same. The solution is somehow make cluster IDs distinguishable. Ideally, this should be done on the substantive basis to set a different prior for each hypothesized latent cluster. In the one dimensional case we have, a common solution seems to constrain \(\mu_{1} \le \mu_{2} \le \dots \le \mu_{H}\).
Implementation
normal_fixed_mixture_stan_code <- biostan::read_stan_file("./bayesianideas_density_normal_fixed_mixture.stan")
biostan::print_stan_code(normal_fixed_mixture_stan_code, section = NULL)## data {
## // Hyperparameters
## real alpha;
## real beta;
## real m;
## real s_squared;
## real<lower=0> dirichlet_alpha;
##
## // Define variables in data
## // Number of observations (an integer)
## int<lower=0> n;
## // Outcome (a real vector of length n)
## real y[n];
## // Number of latent clusters
## int<lower=1> H;
## }
##
## transformed data {
## real s;
## s = sqrt(s_squared);
## }
##
## parameters {
## // Define parameters to estimate
## // Population mean (a real number)
## ordered[H] mu;
## // Population variance (a positive real number)
## real<lower=0> sigma_squared[H];
## // Cluster probability
## simplex[H] Pi;
## }
##
## transformed parameters {
## // Population standard deviation (a positive real number)
## real<lower=0> sigma[H];
## // Standard deviation (derived from variance)
## sigma = sqrt(sigma_squared);
## }
##
## model {
## // Temporary vector for loop use. Need to come first before priors.
## real contributions[H];
##
## // Prior part of Bayesian inference
## // All vectorized
## // Mean
## mu ~ normal(m, s);
## // sigma^2 has inverse gamma (alpha = 1, beta = 1) prior
## sigma_squared ~ inv_gamma(alpha, beta);
## // cluster probability vector
## Pi ~ dirichlet(rep_vector(dirichlet_alpha / H, H));
##
## // Likelihood part of Bayesian inference
## // Outcome model N(mu, sigma^2) (use SD rather than Var)
## for (i in 1:n) {
## // Loop over individuals
## // z[i] in {1,...,H} gives the cluster membership.
## /* y[i] ~ normal(mu[z[i]], sigma[z[i]]); */
##
## for (h in 1:H) {
## // Loop over clusters within each individual
## // Log likelihood contributions log(Pi[h] * N(y[i] | mu[h],sigma[h]))
## contributions[h] = log(Pi[h]) + normal_lpdf(y[i] | mu[h], sigma[h]);
## }
##
## // log(sum(exp(contribution element)))
## target += log_sum_exp(contributions);
##
## }
## }
## Please note the ordering constrain on the \(\mu_{h}\).
Fitting
6 latent cluster model
Let us try a 6 latent cluster model with otherwise similar vague priors.
normal_fixed_mixture_stan_fit6 <-
rstan::stan(model_code = normal_fixed_mixture_stan_code,
data = list(alpha = 10^(-3), beta = 10^(-3),
m = 0, s_squared = 10^(3),
n = nrow(galaxy),
y = galaxy$k_speed,
dirichlet_alpha = 1,
H = 6))## Warning: There were 2097 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Warning: There were 877 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problemsnormal_fixed_mixture_stan_fit6## Inference for Stan model: 81a2617d8fd6a1e89b16e5b25e6f2bf0.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75%
## mu[1] -2.422000e+01 7.41 2.036000e+01 -71.75 -3.752000e+01 -2.077000e+01 -8.980000e+00
## mu[2] -3.710000e+00 5.47 1.421000e+01 -39.38 -1.206000e+01 4.300000e-01 9.300000e+00
## mu[3] 1.097000e+01 1.54 2.730000e+00 9.32 9.630000e+00 9.780000e+00 1.014000e+01
## mu[4] 1.997000e+01 1.72 2.990000e+00 10.83 2.076000e+01 2.125000e+01 2.149000e+01
## mu[5] 2.715000e+01 2.71 4.840000e+00 21.03 2.181000e+01 2.691000e+01 3.213000e+01
## mu[6] 3.810000e+01 4.57 1.259000e+01 25.54 3.224000e+01 3.338000e+01 3.945000e+01
## sigma_squared[1] 1.916584e+305 NaN Inf 2572.96 1.207992e+62 9.890222e+129 1.529336e+206
## sigma_squared[2] 1.113038e+305 NaN Inf 0.12 1.020000e+00 7.272890e+74 3.395365e+172
## sigma_squared[3] 1.333803e+303 NaN Inf 0.09 1.800000e-01 3.100000e-01 1.630000e+00
## sigma_squared[4] 7.661851e+303 NaN Inf 3.38 4.380000e+00 5.150000e+00 7.590000e+00
## sigma_squared[5] 1.269605e+305 NaN Inf 0.54 4.320000e+00 1.084000e+01 5.117357e+83
## sigma_squared[6] 3.624563e+304 NaN Inf 0.42 1.980000e+00 1.908000e+01 4.112080e+21
## Pi[1] 0.000000e+00 0.00 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00
## Pi[2] 2.000000e-02 0.03 4.000000e-02 0.00 0.000000e+00 0.000000e+00 3.000000e-02
## Pi[3] 7.000000e-02 0.03 5.000000e-02 0.00 3.000000e-02 7.000000e-02 1.000000e-01
## Pi[4] 6.400000e-01 0.26 3.700000e-01 0.00 4.900000e-01 8.400000e-01 8.800000e-01
## Pi[5] 2.400000e-01 0.26 3.600000e-01 0.00 0.000000e+00 3.000000e-02 3.500000e-01
## Pi[6] 3.000000e-02 0.01 3.000000e-02 0.00 0.000000e+00 3.000000e-02 5.000000e-02
## sigma[1] 2.999131e+151 NaN 4.368142e+152 50.72 1.099085e+31 9.908190e+64 1.233232e+103
## sigma[2] 1.368041e+151 NaN 3.333833e+152 0.35 1.010000e+00 2.695962e+37 1.837954e+86
## sigma[3] 9.915762e+149 NaN 3.651237e+151 0.30 4.300000e-01 5.500000e-01 1.280000e+00
## sigma[4] 2.632644e+150 NaN 8.750334e+151 1.84 2.090000e+00 2.270000e+00 2.760000e+00
## sigma[5] 1.865682e+151 NaN 3.558709e+152 0.73 2.080000e+00 3.290000e+00 7.118843e+41
## sigma[6] 7.664756e+150 NaN 1.902523e+152 0.65 1.410000e+00 4.370000e+00 6.412413e+10
## lp__ -2.079400e+02 0.18 3.420000e+00 -215.56 -2.100200e+02 -2.076400e+02 -2.055200e+02
## 97.5% n_eff Rhat
## mu[1] 6.840000e+00 8 1.22
## mu[2] 9.970000e+00 7 1.23
## mu[3] 1.948000e+01 3 1.79
## mu[4] 2.188000e+01 3 1.78
## mu[5] 3.455000e+01 3 1.80
## mu[6] 7.815000e+01 8 1.30
## sigma_squared[1] 4.883207e+299 NaN NaN
## sigma_squared[2] 1.344096e+291 NaN NaN
## sigma_squared[3] 8.024482e+269 NaN NaN
## sigma_squared[4] 6.289566e+254 NaN NaN
## sigma_squared[5] 1.782747e+279 NaN NaN
## sigma_squared[6] 1.459042e+276 NaN NaN
## Pi[1] 1.000000e-02 2124 1.00
## Pi[2] 1.300000e-01 2 2.72
## Pi[3] 1.500000e-01 3 1.79
## Pi[4] 9.300000e-01 2 9.90
## Pi[5] 9.100000e-01 2 12.10
## Pi[6] 1.100000e-01 7 1.37
## sigma[1] 6.980213e+149 NaN 1.00
## sigma[2] 3.659495e+145 NaN 1.00
## sigma[3] 8.935240e+134 NaN 1.00
## sigma[4] 2.503350e+127 NaN 1.00
## sigma[5] 4.125863e+139 NaN 1.00
## sigma[6] 1.207856e+138 NaN 1.00
## lp__ -2.022000e+02 352 1.03
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:38:20 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all <- function(fit, H) {
print(traceplot(fit, inc_warmup = TRUE, pars = sprintf("mu[%s]", seq_len(H))))
print(traceplot(fit, inc_warmup = TRUE, pars = sprintf("sigma[%s]", seq_len(H))))
print(traceplot(fit, inc_warmup = TRUE, pars = sprintf("Pi[%s]", seq_len(H))))
}
traceplot_all(normal_fixed_mixture_stan_fit6, 6)
This model broke down badly.
1 latent cluster model
Let us sanity check with just one cluster.
normal_fixed_mixture_stan_fit1 <-
rstan::stan(model_code = normal_fixed_mixture_stan_code,
data = list(alpha = 10^(-3), beta = 10^(-3),
m = 0, s_squared = 10^(3),
n = nrow(galaxy),
y = galaxy$k_speed,
dirichlet_alpha = 1,
H = 1))normal_fixed_mixture_stan_fit1## Inference for Stan model: 81a2617d8fd6a1e89b16e5b25e6f2bf0.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu[1] 20.84 0.01 0.51 19.85 20.48 20.84 21.19 21.85 3514 1
## sigma_squared[1] 21.32 0.06 3.36 15.73 18.89 21.05 23.45 28.44 2709 1
## Pi[1] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 NaN NaN
## sigma[1] 4.60 0.01 0.36 3.97 4.35 4.59 4.84 5.33 2708 1
## lp__ -241.64 0.02 0.97 -244.32 -242.04 -241.35 -240.94 -240.66 1798 1
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:38:32 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit1, 1)
This gave a similar result to the first single normal fit, except that the additional \(\pi\) parameter, which can only be 1 in this case.
2 latent cluster model
This model has two latent clusters.
normal_fixed_mixture_stan_fit2 <-
rstan::stan(model_code = normal_fixed_mixture_stan_code,
data = list(alpha = 10^(-3), beta = 10^(-3),
m = 0, s_squared = 10^(3),
n = nrow(galaxy),
y = galaxy$k_speed,
dirichlet_alpha = 1,
H = 2))normal_fixed_mixture_stan_fit2## Inference for Stan model: 81a2617d8fd6a1e89b16e5b25e6f2bf0.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu[1] 14.91 3.72 5.31 9.40 9.70 13.19 20.82 21.61 2 8.20
## mu[2] 21.86 0.25 0.69 20.97 21.43 21.75 22.11 23.61 8 1.18
## sigma_squared[1] 20.33 23.53 36.37 0.10 0.23 1.89 12.32 118.28 2 2.44
## sigma_squared[2] 28.67 25.61 41.62 2.77 6.79 10.20 21.64 146.14 3 2.01
## Pi[1] 0.29 0.19 0.27 0.04 0.09 0.15 0.50 0.83 2 4.95
## Pi[2] 0.71 0.19 0.27 0.17 0.50 0.85 0.91 0.96 2 4.95
## sigma[1] 2.89 2.39 3.46 0.31 0.48 1.37 3.32 10.88 2 4.37
## sigma[2] 4.41 2.03 3.03 1.66 2.61 3.19 4.64 12.09 2 3.05
## lp__ -223.09 1.08 2.18 -227.94 -224.55 -222.75 -221.37 -219.86 4 1.40
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:38:39 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit2, 2)
This model also exhibits poor mixing. Particularly, \(\sigma^{2}_{1}\) has a extremely broad posterior distribution.
3 latent cluster model
normal_fixed_mixture_stan_fit3 <-
rstan::stan(model_code = normal_fixed_mixture_stan_code,
data = list(alpha = 10^(-3), beta = 10^(-3),
m = 0, s_squared = 10^(3),
n = nrow(galaxy),
y = galaxy$k_speed,
dirichlet_alpha = 1,
H = 3))## Warning: There were 769 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Warning: There were 357 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problemsnormal_fixed_mixture_stan_fit3## Inference for Stan model: 81a2617d8fd6a1e89b16e5b25e6f2bf0.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu[1] 5.210000e+00 8.99 1.544000e+01 -40.92 9.12 9.71 10.74 1.974000e+01 3 1.83
## mu[2] 1.798000e+01 3.42 4.900000e+00 9.47 11.38 20.93 21.38 2.182000e+01 2 6.24
## mu[3] 2.665000e+01 3.60 5.470000e+00 21.02 21.58 22.68 32.60 3.469000e+01 2 2.73
## sigma_squared[1] 1.473928e+303 NaN Inf 0.10 0.23 1.78 194.90 6.083956e+257 NaN NaN
## sigma_squared[2] 1.270297e+303 NaN Inf 0.12 2.02 4.89 8.04 1.898486e+248 NaN NaN
## sigma_squared[3] 1.335000e+01 1.62 5.217000e+01 0.45 2.93 5.18 10.99 6.987000e+01 1043 1.01
## Pi[1] 1.000000e-01 0.06 1.000000e-01 0.00 0.03 0.08 0.15 3.400000e-01 2 2.39
## Pi[2] 4.500000e-01 0.29 4.100000e-01 0.00 0.03 0.38 0.87 9.200000e-01 2 11.56
## Pi[3] 4.400000e-01 0.28 4.000000e-01 0.01 0.04 0.43 0.86 9.400000e-01 2 8.85
## sigma[1] 1.031206e+150 NaN 3.838272e+151 0.32 0.48 1.33 13.96 7.798720e+128 NaN 1.00
## sigma[2] 1.003322e+150 NaN 3.563156e+151 0.34 1.42 2.21 2.84 1.377584e+124 NaN 1.00
## sigma[3] 2.850000e+00 0.18 2.280000e+00 0.67 1.71 2.28 3.32 8.360000e+00 163 1.04
## lp__ -2.155700e+02 6.38 9.260000e+00 -229.52 -224.61 -217.83 -206.46 -2.032100e+02 2 4.36
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:39:11 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit3, 3)
Adding another latent cluster does not help.
2 latent cluster model with increased Dirichlet alpha
normal_fixed_mixture_stan_fit2a <-
rstan::stan(model_code = normal_fixed_mixture_stan_code,
data = list(alpha = 10^(-3), beta = 10^(-3),
m = 0, s_squared = 10^(3),
n = nrow(galaxy),
y = galaxy$k_speed,
dirichlet_alpha = 2,
H = 2))normal_fixed_mixture_stan_fit2a## Inference for Stan model: 81a2617d8fd6a1e89b16e5b25e6f2bf0.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu[1] 18.35 3.53 4.99 9.47 17.85 21.12 21.35 21.77 2 19.76
## mu[2] 22.22 0.12 1.00 21.02 21.59 21.99 22.54 24.86 71 1.03
## sigma_squared[1] 2.86 1.06 1.68 0.12 1.80 3.23 4.00 5.67 3 2.14
## sigma_squared[2] 68.16 24.73 49.47 8.33 25.90 64.85 92.59 185.34 4 1.39
## Pi[1] 0.57 0.20 0.28 0.06 0.33 0.70 0.77 0.86 2 4.33
## Pi[2] 0.43 0.20 0.28 0.14 0.23 0.30 0.67 0.94 2 4.33
## sigma[1] 1.56 0.43 0.64 0.34 1.34 1.80 2.00 2.38 2 2.95
## sigma[2] 7.66 1.87 3.09 2.89 5.08 8.05 9.62 13.61 3 1.86
## lp__ -225.52 0.99 2.17 -229.99 -226.77 -225.49 -224.35 -221.36 5 1.31
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:39:19 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit2a, 2)
Increasing the mass parameter for the Dirichlet prior on the cluster probabilities does not seem to help.
2 latent cluster model with more precise prior on means
normal_fixed_mixture_stan_fit2b <-
rstan::stan(model_code = normal_fixed_mixture_stan_code,
data = list(alpha = 10^(-3), beta = 10^(-3),
m = 0, s_squared = 1,
n = nrow(galaxy),
y = galaxy$k_speed,
dirichlet_alpha = 1,
H = 2))## Warning: There were 2511 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Warning: Examine the pairs() plot to diagnose sampling problemsnormal_fixed_mixture_stan_fit2b## Inference for Stan model: 81a2617d8fd6a1e89b16e5b25e6f2bf0.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75%
## mu[1] -7.000000e-02 0.03 1.000000e+00 -2.02 -7.600000e-01 -8.000000e-02 6.300000e-01
## mu[2] 4.730000e+00 0.04 1.170000e+00 2.42 3.930000e+00 4.750000e+00 5.530000e+00
## sigma_squared[1] 2.085790e+305 NaN Inf 18673.09 2.247663e+58 2.045329e+130 8.156128e+210
## sigma_squared[2] 2.897300e+02 2.21 6.128000e+01 184.97 2.459800e+02 2.838300e+02 3.264500e+02
## Pi[1] 1.000000e-02 0.00 1.000000e-02 0.00 0.000000e+00 0.000000e+00 1.000000e-02
## Pi[2] 9.900000e-01 0.00 1.000000e-02 0.97 9.900000e-01 1.000000e+00 1.000000e+00
## sigma[1] 2.320744e+151 NaN 4.561715e+152 136.65 1.499209e+29 1.430113e+65 2.843413e+105
## sigma[2] 1.693000e+01 0.06 1.780000e+00 13.60 1.568000e+01 1.685000e+01 1.807000e+01
## lp__ -3.626000e+02 0.06 1.560000e+00 -366.41 -3.634600e+02 -3.622700e+02 -3.614500e+02
## 97.5% n_eff Rhat
## mu[1] 1.830000e+00 816 1.00
## mu[2] 7.000000e+00 948 1.00
## sigma_squared[1] 4.192219e+293 NaN NaN
## sigma_squared[2] 4.261800e+02 767 1.00
## Pi[1] 3.000000e-02 1807 1.00
## Pi[2] 1.000000e+00 1807 1.00
## sigma[1] 6.474488e+146 NaN 1.00
## sigma[2] 2.064000e+01 782 1.00
## lp__ -3.605200e+02 693 1.01
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:39:28 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit2b, 2)
2 latent cluster model with more precise prior on means and variance
normal_fixed_mixture_stan_fit2c <-
rstan::stan(model_code = normal_fixed_mixture_stan_code,
data = list(alpha = 10^(1), beta = 10^(1),
m = 0, s_squared = 1,
n = nrow(galaxy),
y = galaxy$k_speed,
dirichlet_alpha = 1,
H = 2))normal_fixed_mixture_stan_fit2c## Inference for Stan model: 81a2617d8fd6a1e89b16e5b25e6f2bf0.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu[1] 0.01 0.02 1.01 -1.90 -0.68 0.00 0.68 2.02 2408 1
## mu[2] 6.43 0.02 1.28 4.01 5.57 6.39 7.29 8.99 2738 1
## sigma_squared[1] 1.12 0.01 0.41 0.58 0.84 1.03 1.32 2.15 1841 1
## sigma_squared[2] 188.63 0.79 39.63 120.63 160.70 185.91 213.14 272.45 2489 1
## Pi[1] 0.01 0.00 0.01 0.00 0.00 0.00 0.01 0.03 4566 1
## Pi[2] 0.99 0.00 0.01 0.97 0.99 1.00 1.00 1.00 4566 1
## sigma[1] 1.04 0.00 0.18 0.76 0.92 1.02 1.15 1.47 2053 1
## sigma[2] 13.66 0.03 1.43 10.98 12.68 13.63 14.60 16.51 2466 1
## lp__ -426.78 0.04 1.71 -431.02 -427.68 -426.44 -425.52 -424.48 1537 1
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:39:37 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit2c, 2)
Additionally making the prior for variance more precise did not help.
2 latent cluster model with more chains and iterations
normal_fixed_mixture_stan_fit2d <-
rstan::stan(model_code = normal_fixed_mixture_stan_code,
data = list(alpha = 10^(-3), beta = 10^(-3),
m = 0, s_squared = 10^(3),
n = nrow(galaxy),
y = galaxy$k_speed,
dirichlet_alpha = 1,
H = 2),
chains = 12,
iter = 10000)## Warning: There were 3061 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Warning: Examine the pairs() plot to diagnose sampling problemsnormal_fixed_mixture_stan_fit2d## Inference for Stan model: 81a2617d8fd6a1e89b16e5b25e6f2bf0.
## 12 chains, each with iter=10000; warmup=5000; thin=1;
## post-warmup draws per chain=5000, total post-warmup draws=60000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu[1] 1.825000e+01 1.62 4.100000e+00 9.50 17.47 20.09 21.12 2.166000e+01 6 4.16
## mu[2] 2.329000e+01 2.00 6.750000e+00 20.89 21.34 21.66 22.22 4.507000e+01 11 1.49
## sigma_squared[1] 3.635000e+01 15.26 4.342000e+01 0.13 3.04 11.07 65.74 1.379300e+02 8 1.96
## sigma_squared[2] 1.246745e+304 NaN Inf 2.56 3.95 10.17 77.00 5.256578e+188 NaN NaN
## Pi[1] 4.500000e-01 0.12 3.100000e-01 0.06 0.19 0.32 0.74 1.000000e+00 6 4.43
## Pi[2] 5.500000e-01 0.12 3.100000e-01 0.00 0.26 0.68 0.81 9.400000e-01 6 4.43
## sigma[1] 4.760000e+00 1.44 3.700000e+00 0.36 1.74 3.32 8.11 1.174000e+01 7 3.27
## sigma[2] 1.940740e+150 NaN 1.116418e+152 1.60 1.99 3.19 8.78 2.292520e+94 NaN 1.00
## lp__ -2.255600e+02 2.34 5.980000e+00 -244.56 -225.75 -223.96 -222.34 -2.203300e+02 7 3.54
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:39:52 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit2d, 2)
2 latent cluster model with tau (precision) ~ Gamma(a,b) parametrization
normal_fixed_mixture_tau_stan_code <- biostan::read_stan_file("./bayesianideas_density_normal_fixed_mixture_tau.stan")
biostan::print_stan_code(normal_fixed_mixture_tau_stan_code, section = NULL)## data {
## // Hyperparameters
## real alpha;
## real beta;
## real m;
## real<lower=0> s_squared;
## real<lower=0> dirichlet_alpha;
##
## // Define variables in data
## // Number of observations (an integer)
## int<lower=0> n;
## // Outcome (a real vector of length n)
## real y[n];
## // Number of latent clusters
## int<lower=1> H;
## }
##
## transformed data {
## real s;
## s = sqrt(s_squared);
## }
##
## parameters {
## // Define parameters to estimate
## // Population mean (a real number)
## ordered[H] mu;
## // Population variance (a positive real number)
## real<lower=0> tau[H];
## // Cluster probability
## simplex[H] Pi;
## }
##
## transformed parameters {
## // Population standard deviation (a positive real number)
## real<lower=0> sigma[H];
## // Standard deviation (derived from variance)
## for (h in 1:H) {
## sigma[h] = sqrt(1 / tau[h]);
## }
## }
##
## model {
## // Temporary vector for loop use. Need to come first before priors.
## real contributions[H];
##
## // Prior part of Bayesian inference
## // All vectorized
## // Mean
## mu ~ normal(m, s);
## // tau = 1/sigma^2 has gamma prior
## tau ~ gamma(alpha, beta);
## // cluster probability vector
## Pi ~ dirichlet(rep_vector(dirichlet_alpha / H, H));
##
## // Likelihood part of Bayesian inference
## // Outcome model N(mu, sigma^2) (use SD rather than Var)
## for (i in 1:n) {
## // Loop over individuals
## // z[i] in {1,...,H} gives the cluster membership.
## /* y[i] ~ normal(mu[z[i]], sigma[z[i]]); */
##
## for (h in 1:H) {
## // Loop over clusters within each individual
## // Log likelihood contributions log(Pi[h] * N(y[i] | mu[h],sigma[h]))
## contributions[h] = log(Pi[h]) + normal_lpdf(y[i] | mu[h], sigma[h]);
## }
##
## // log(sum(exp(contribution element)))
## target += log_sum_exp(contributions);
##
## }
## }
## normal_fixed_mixture_stan_fit2e <-
rstan::stan(model_code = normal_fixed_mixture_tau_stan_code,
data = list(alpha = 10^(-3), beta = 10^(-3),
m = 0, s_squared = 10^(3),
n = nrow(galaxy),
y = galaxy$k_speed,
dirichlet_alpha = 2,
H = 2))normal_fixed_mixture_stan_fit2e## Inference for Stan model: 5785463fcca3d68ec4c4ac44df3fa279.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu[1] 19.35 0.78 1.77 15.44 18.13 19.67 20.93 21.59 5 1.31
## mu[2] 21.61 0.32 0.75 20.81 21.22 21.44 21.72 23.75 6 1.25
## tau[1] 0.08 0.09 0.13 0.01 0.01 0.02 0.06 0.39 2 3.51
## tau[2] 0.21 0.08 0.13 0.01 0.10 0.24 0.30 0.42 2 2.34
## Pi[1] 0.38 0.14 0.21 0.14 0.23 0.29 0.50 0.83 2 2.95
## Pi[2] 0.62 0.14 0.21 0.17 0.50 0.71 0.77 0.86 2 2.95
## sigma[1] 7.01 2.13 3.32 1.60 4.60 7.82 9.12 12.64 2 2.31
## sigma[2] 3.72 2.22 3.30 1.54 1.82 2.02 3.36 11.86 2 3.23
## lp__ -224.71 0.75 1.94 -229.37 -225.80 -224.44 -223.21 -221.96 7 1.21
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:41:01 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit2e, 2)
2 latent cluster model on original scale
normal_fixed_mixture_stan_fit2f <-
rstan::stan(model_code = normal_fixed_mixture_stan_code,
data = list(alpha = 10^(-3), beta = 10^(-3),
m = 0, s_squared = 10^(3),
n = nrow(galaxy),
y = galaxy$speed,
dirichlet_alpha = 1,
H = 2))## Warning: There were 2294 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Warning: Examine the pairs() plot to diagnose sampling problemsnormal_fixed_mixture_stan_fit2f## Inference for Stan model: 81a2617d8fd6a1e89b16e5b25e6f2bf0.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75%
## mu[1] -1.718000e+01 0.91 2.604000e+01 -71.29 -34.43 -16.19 8.000000e-01
## mu[2] 1.866000e+01 0.67 2.637000e+01 -30.24 1.32 17.65 3.541000e+01
## sigma_squared[1] 9.674271e+304 NaN Inf 2253892.88 447903069.27 598067102.04 4.866345e+118
## sigma_squared[2] 9.458461e+304 NaN Inf 322136077.94 449097279.12 619141633.42 2.836931e+129
## Pi[1] 5.000000e-01 0.35 4.900000e-01 0.00 0.00 0.50 1.000000e+00
## Pi[2] 5.000000e-01 0.35 4.900000e-01 0.00 0.00 0.50 1.000000e+00
## sigma[1] 1.496723e+151 NaN 3.107134e+152 1501.11 21163.72 24455.41 2.205482e+59
## sigma[2] 1.762656e+151 NaN 3.070790e+152 17948.15 21191.92 24882.55 5.300077e+64
## lp__ -9.361600e+02 0.08 1.770000e+00 -940.75 -936.99 -935.75 -9.348900e+02
## 97.5% n_eff Rhat
## mu[1] 3.056000e+01 816 1.01
## mu[2] 7.396000e+01 1566 1.00
## sigma_squared[1] 1.195252e+284 NaN NaN
## sigma_squared[2] 1.406913e+288 NaN NaN
## Pi[1] 1.000000e+00 2 63.56
## Pi[2] 1.000000e+00 2 63.56
## sigma[1] 8.750030e+141 NaN 1.00
## sigma[2] 1.183391e+144 NaN 1.00
## lp__ -9.339500e+02 501 1.01
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:41:11 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit2f, 2)
2 latent cluster model on original scale with smaller steps
normal_fixed_mixture_stan_fit2g <-
rstan::stan(model_code = normal_fixed_mixture_stan_code,
data = list(alpha = 10^(-3), beta = 10^(-3),
m = 0, s_squared = 10^(3),
n = nrow(galaxy),
y = galaxy$speed,
dirichlet_alpha = 1,
H = 2),
control = list(adapt_delta = 0.99999999999))## Warning: There were 881 divergent transitions after warmup. Increasing adapt_delta above 0.99999999999 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problemsnormal_fixed_mixture_stan_fit2g## Inference for Stan model: 81a2617d8fd6a1e89b16e5b25e6f2bf0.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75%
## mu[1] -1.668000e+01 0.82 2.618000e+01 -69.23 -3.362000e+01 -1.639000e+01 1.060000e+00
## mu[2] 2.068000e+01 0.62 2.663000e+01 -29.03 2.400000e+00 1.971000e+01 3.738000e+01
## sigma_squared[1] 3.456769e+305 NaN Inf 282.17 3.045594e+52 2.427263e+125 9.218386e+205
## sigma_squared[2] 4.622979e+08 2000917.10 7.392890e+07 339411787.76 4.083811e+08 4.547803e+08 5.072759e+08
## Pi[1] 1.000000e-02 0.00 1.000000e-02 0.00 0.000000e+00 0.000000e+00 1.000000e-02
## Pi[2] 9.900000e-01 0.00 1.000000e-02 0.97 9.900000e-01 1.000000e+00 1.000000e+00
## sigma[1] 5.600726e+151 NaN 5.853424e+152 16.80 1.743782e+26 4.925787e+62 9.594199e+102
## sigma[2] 2.143411e+04 45.85 1.696330e+03 18423.13 2.020844e+04 2.132558e+04 2.252279e+04
## lp__ -9.359200e+02 0.07 1.650000e+00 -940.19 -9.366500e+02 -9.355700e+02 -9.347500e+02
## 97.5% n_eff Rhat
## mu[1] 3.509000e+01 1014 1.00
## mu[2] 7.620000e+01 1838 1.00
## sigma_squared[1] 1.349108e+302 NaN NaN
## sigma_squared[2] 6.302219e+08 1365 1.00
## Pi[1] 3.000000e-02 824 1.00
## Pi[2] 1.000000e+00 824 1.00
## sigma[1] 1.161489e+151 NaN 1.00
## sigma[2] 2.510422e+04 1369 1.00
## lp__ -9.338500e+02 618 1.01
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:41:35 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit2g, 2)
2 latent cluster model with log-normal prior on sigma
normal_fixed_mixture_lognormal_stan_code <- biostan::read_stan_file("./bayesianideas_density_normal_fixed_mixture_lognormal.stan")
biostan::print_stan_code(normal_fixed_mixture_lognormal_stan_code, section = NULL)## data {
## // Hyperparameters
## real mean_log_normal;
## real<lower=0> sd_log_normal;
## real m;
## real<lower=0> s_squared;
## real<lower=0> dirichlet_alpha;
##
## // Define variables in data
## // Number of observations (an integer)
## int<lower=0> n;
## // Outcome (a real vector of length n)
## real y[n];
## // Number of latent clusters
## int<lower=1> H;
## }
##
## transformed data {
## real s;
## s = sqrt(s_squared);
## }
##
## parameters {
## // Define parameters to estimate
## // Population mean (a real number)
## ordered[H] mu;
## // Population variance (a positive real number)
## real<lower=0> sigma[H];
## // Cluster probability
## simplex[H] Pi;
## }
##
## model {
## // Temporary vector for loop use. Need to come first before priors.
## real contributions[H];
##
## // Prior part of Bayesian inference
## // All vectorized
## // Mean
## mu ~ normal(m, s);
## //
## sigma ~ lognormal(mean_log_normal, sd_log_normal);
## // cluster probability vector
## Pi ~ dirichlet(rep_vector(dirichlet_alpha / H, H));
##
## // Likelihood part of Bayesian inference
## // Outcome model N(mu, sigma^2) (use SD rather than Var)
## for (i in 1:n) {
## // Loop over individuals
## // z[i] in {1,...,H} gives the cluster membership.
## /* y[i] ~ normal(mu[z[i]], sigma[z[i]]); */
##
## for (h in 1:H) {
## // Loop over clusters within each individual
## // Log likelihood contributions log(Pi[h] * N(y[i] | mu[h],sigma[h]))
## contributions[h] = log(Pi[h]) + normal_lpdf(y[i] | mu[h], sigma[h]);
## }
##
## // log(sum(exp(contribution element)))
## target += log_sum_exp(contributions);
##
## }
## }
## normal_fixed_mixture_stan_fit2h <-
rstan::stan(model_code = normal_fixed_mixture_lognormal_stan_code,
data = list(mean_log_normal = 0, sd_log_normal = 2,
m = 0, s_squared = 10^(3),
n = nrow(galaxy),
y = galaxy$k_speed,
dirichlet_alpha = 1,
H = 2),
control = list(adapt_delta = 0.99999999999))normal_fixed_mixture_stan_fit2h## Inference for Stan model: 7c962c9cc4a3894467631c06c6c28588.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu[1] 14.14 3.15 4.63 9.37 9.71 10.44 18.87 20.92 2 3.73
## mu[2] 21.62 0.19 0.42 20.89 21.32 21.58 21.90 22.51 5 1.30
## sigma[1] 4.62 2.92 4.27 0.32 0.49 3.33 8.44 11.94 2 3.69
## sigma[2] 2.57 0.45 0.68 1.60 1.93 2.60 3.18 3.67 2 2.93
## Pi[1] 0.17 0.06 0.10 0.04 0.09 0.15 0.25 0.38 3 1.86
## Pi[2] 0.83 0.06 0.10 0.62 0.75 0.85 0.91 0.96 3 1.86
## lp__ -223.00 0.63 1.84 -227.32 -224.03 -222.78 -221.74 -220.09 9 1.16
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:42:21 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit2h, 2)
2 latent cluster model with more precise prior on means and variance (density grid)
normal_fixed_mixture_grid_stan_code <- biostan::read_stan_file("./bayesianideas_density_normal_fixed_mixture_grid.stan")
biostan::print_stan_code(normal_fixed_mixture_grid_stan_code, section = NULL)## data {
## // Hyperparameters
## real alpha;
## real beta;
## real m;
## real s_squared;
## real<lower=0> dirichlet_alpha;
##
## // Define variables in data
## // Number of observations (an integer)
## int<lower=0> n;
## // Outcome (a real vector of length n)
## real y[n];
## // Number of latent clusters
## int<lower=1> H;
##
## // Grid evaluation
## real grid_max;
## real grid_min;
## int<lower=1> grid_length;
## }
##
## transformed data {
## real s;
## real grid_step;
##
## s = sqrt(s_squared);
## grid_step = (grid_max - grid_min) / (grid_length - 1);
## }
##
## parameters {
## // Define parameters to estimate
## // Population mean (a real number)
## ordered[H] mu;
## // Population variance (a positive real number)
## real<lower=0> sigma_squared[H];
## // Cluster probability
## simplex[H] Pi;
## }
##
## transformed parameters {
## // Population standard deviation (a positive real number)
## real<lower=0> sigma[H];
## // Standard deviation (derived from variance)
## sigma = sqrt(sigma_squared);
## }
##
## model {
## // Temporary vector for loop use. Need to come first before priors.
## real contributions[H];
##
## // Prior part of Bayesian inference
## // All vectorized
## // Mean
## mu ~ normal(m, s);
## // sigma^2 has inverse gamma (alpha = 1, beta = 1) prior
## sigma_squared ~ inv_gamma(alpha, beta);
## // cluster probability vector
## Pi ~ dirichlet(rep_vector(dirichlet_alpha / H, H));
##
## // Likelihood part of Bayesian inference
## // Outcome model N(mu, sigma^2) (use SD rather than Var)
## for (i in 1:n) {
## // Loop over individuals
## // z[i] in {1,...,H} gives the cluster membership.
## /* y[i] ~ normal(mu[z[i]], sigma[z[i]]); */
##
## for (h in 1:H) {
## // Loop over clusters within each individual
## // Log likelihood contributions log(Pi[h] * N(y[i] | mu[h],sigma[h]))
## contributions[h] = log(Pi[h]) + normal_lpdf(y[i] | mu[h], sigma[h]);
## }
##
## // log(sum(exp(contribution element)))
## target += log_sum_exp(contributions);
##
## }
## }
##
## generated quantities {
##
## real log_f[grid_length];
##
## for (g in 1:grid_length) {
## // Definiting here avoids reporting of these intermediates.
## real contributions[H];
## real grid_value;
##
## grid_value = grid_min + grid_step * (g - 1);
## for (h in 1:H) {
## contributions[h] = log(Pi[h]) + normal_lpdf(grid_value | mu[h], sigma[h]);
## }
##
## log_f[g] = log_sum_exp(contributions);
## }
##
## }
## grid_max <- 40
grid_min <- -20
grid_length <- 100normal_fixed_mixture_stan_fit2i <-
rstan::stan(model_code = normal_fixed_mixture_grid_stan_code,
data = list(alpha = 10^(1), beta = 10^(1),
m = 0, s_squared = 1,
n = nrow(galaxy),
y = galaxy$k_speed,
dirichlet_alpha = 1,
H = 2,
grid_max = grid_max, grid_min = grid_min, grid_length = grid_length))normal_fixed_mixture_stan_fit2i## Inference for Stan model: 963ec46299674d2582fc7e9c128240d0.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu[1] 0.62 0.77 1.42 -1.84 -0.43 0.48 1.67 3.37 3 1.53
## mu[2] 5.55 1.17 1.98 1.82 3.89 5.90 7.02 8.76 3 1.76
## sigma_squared[1] 74.10 89.67 128.62 0.61 0.89 1.18 42.81 356.89 2 5.86
## sigma_squared[2] 140.67 57.85 87.73 0.69 60.36 167.08 200.37 269.86 2 2.68
## Pi[1] 0.25 0.30 0.43 0.00 0.00 0.01 0.29 1.00 2 53.94
## Pi[2] 0.75 0.30 0.43 0.00 0.71 0.99 1.00 1.00 2 53.94
## sigma[1] 5.04 4.91 6.98 0.78 0.95 1.09 4.54 18.89 2 10.78
## sigma[2] 10.46 3.87 5.59 0.83 7.11 12.93 14.16 16.43 2 4.74
## log_f[1] -5.30 0.28 0.58 -6.62 -5.65 -5.26 -4.81 -4.52 4 1.37
## log_f[2] -5.22 0.27 0.55 -6.48 -5.55 -5.18 -4.75 -4.48 4 1.37
## log_f[3] -5.15 0.25 0.52 -6.35 -5.46 -5.10 -4.70 -4.45 4 1.37
## log_f[4] -5.07 0.24 0.50 -6.22 -5.37 -5.03 -4.65 -4.41 4 1.36
## log_f[5] -5.00 0.23 0.47 -6.09 -5.28 -4.96 -4.59 -4.38 4 1.36
## log_f[6] -4.93 0.22 0.45 -5.96 -5.19 -4.89 -4.54 -4.34 4 1.36
## log_f[7] -4.86 0.20 0.42 -5.83 -5.11 -4.82 -4.50 -4.31 4 1.36
## log_f[8] -4.79 0.19 0.40 -5.71 -5.02 -4.75 -4.45 -4.28 4 1.35
## log_f[9] -4.73 0.18 0.38 -5.59 -4.94 -4.69 -4.40 -4.24 4 1.35
## log_f[10] -4.66 0.17 0.35 -5.48 -4.86 -4.62 -4.35 -4.21 4 1.35
## log_f[11] -4.60 0.16 0.33 -5.37 -4.79 -4.56 -4.31 -4.18 4 1.34
## log_f[12] -4.54 0.15 0.31 -5.26 -4.71 -4.51 -4.27 -4.15 5 1.34
## log_f[13] -4.48 0.13 0.29 -5.15 -4.64 -4.45 -4.23 -4.12 5 1.33
## log_f[14] -4.42 0.12 0.27 -5.05 -4.57 -4.39 -4.19 -4.09 5 1.32
## log_f[15] -4.37 0.11 0.25 -4.95 -4.51 -4.34 -4.16 -4.06 5 1.31
## log_f[16] -4.32 0.10 0.23 -4.86 -4.44 -4.29 -4.13 -4.03 5 1.30
## log_f[17] -4.27 0.09 0.21 -4.76 -4.38 -4.24 -4.10 -4.00 5 1.29
## log_f[18] -4.22 0.08 0.19 -4.67 -4.32 -4.19 -4.07 -3.97 5 1.28
## log_f[19] -4.17 0.07 0.17 -4.59 -4.26 -4.14 -4.04 -3.94 5 1.26
## log_f[20] -4.13 0.07 0.16 -4.50 -4.20 -4.10 -4.01 -3.92 6 1.24
## log_f[21] -4.08 0.06 0.14 -4.42 -4.15 -4.06 -3.98 -3.90 6 1.22
## log_f[22] -4.04 0.05 0.13 -4.35 -4.10 -4.02 -3.95 -3.87 7 1.19
## log_f[23] -4.00 0.04 0.11 -4.27 -4.05 -3.98 -3.93 -3.85 8 1.16
## log_f[24] -3.96 0.03 0.10 -4.20 -4.00 -3.95 -3.90 -3.82 10 1.13
## log_f[25] -3.93 0.02 0.09 -4.13 -3.96 -3.92 -3.87 -3.80 15 1.09
## log_f[26] -3.89 0.01 0.08 -4.07 -3.93 -3.88 -3.84 -3.77 31 1.05
## log_f[27] -3.86 0.00 0.07 -4.01 -3.89 -3.85 -3.82 -3.75 463 1.02
## log_f[28] -3.83 0.00 0.06 -3.96 -3.86 -3.82 -3.79 -3.72 2941 1.00
## log_f[29] -3.80 0.00 0.06 -3.93 -3.83 -3.79 -3.76 -3.69 481 1.02
## log_f[30] -3.76 0.02 0.07 -3.91 -3.80 -3.76 -3.72 -3.63 18 1.07
## log_f[31] -3.73 0.03 0.08 -3.89 -3.77 -3.73 -3.68 -3.55 9 1.14
## log_f[32] -3.69 0.04 0.10 -3.88 -3.75 -3.70 -3.65 -3.47 6 1.20
## log_f[33] -3.66 0.05 0.11 -3.87 -3.73 -3.67 -3.61 -3.42 5 1.26
## log_f[34] -3.64 0.05 0.12 -3.86 -3.71 -3.64 -3.58 -3.38 5 1.30
## log_f[35] -3.63 0.05 0.11 -3.84 -3.70 -3.62 -3.56 -3.38 5 1.30
## log_f[36] -3.62 0.05 0.11 -3.83 -3.68 -3.61 -3.55 -3.40 5 1.27
## log_f[37] -3.61 0.04 0.10 -3.82 -3.67 -3.60 -3.54 -3.42 6 1.21
## log_f[38] -3.60 0.04 0.10 -3.82 -3.66 -3.59 -3.53 -3.42 8 1.16
## log_f[39] -3.60 0.04 0.11 -3.82 -3.67 -3.58 -3.52 -3.41 8 1.15
## log_f[40] -3.59 0.05 0.12 -3.84 -3.67 -3.58 -3.51 -3.40 6 1.21
## log_f[41] -3.59 0.06 0.12 -3.86 -3.68 -3.58 -3.50 -3.39 5 1.31
## log_f[42] -3.59 0.07 0.13 -3.87 -3.68 -3.58 -3.50 -3.38 4 1.41
## log_f[43] -3.60 0.07 0.14 -3.88 -3.69 -3.58 -3.49 -3.37 4 1.49
## log_f[44] -3.60 0.08 0.14 -3.89 -3.70 -3.58 -3.49 -3.35 3 1.53
## log_f[45] -3.60 0.08 0.15 -3.90 -3.71 -3.58 -3.49 -3.35 3 1.55
## log_f[46] -3.60 0.08 0.15 -3.90 -3.72 -3.58 -3.49 -3.34 3 1.56
## log_f[47] -3.61 0.09 0.16 -3.91 -3.72 -3.59 -3.49 -3.34 3 1.57
## log_f[48] -3.61 0.09 0.16 -3.93 -3.74 -3.59 -3.49 -3.34 3 1.58
## log_f[49] -3.62 0.09 0.16 -3.94 -3.75 -3.60 -3.50 -3.34 3 1.59
## log_f[50] -3.63 0.09 0.17 -3.95 -3.76 -3.61 -3.51 -3.34 3 1.60
## log_f[51] -3.65 0.09 0.17 -3.97 -3.78 -3.62 -3.52 -3.35 3 1.61
## log_f[52] -3.66 0.10 0.17 -3.98 -3.79 -3.64 -3.54 -3.36 3 1.61
## log_f[53] -3.68 0.10 0.17 -4.00 -3.81 -3.65 -3.55 -3.37 3 1.62
## log_f[54] -3.69 0.10 0.17 -4.02 -3.83 -3.67 -3.57 -3.39 3 1.63
## log_f[55] -3.71 0.10 0.17 -4.03 -3.85 -3.69 -3.59 -3.41 3 1.64
## log_f[56] -3.74 0.10 0.17 -4.05 -3.88 -3.71 -3.61 -3.43 3 1.65
## log_f[57] -3.76 0.10 0.17 -4.07 -3.90 -3.74 -3.63 -3.46 3 1.66
## log_f[58] -3.78 0.10 0.17 -4.10 -3.92 -3.76 -3.66 -3.48 3 1.67
## log_f[59] -3.81 0.10 0.17 -4.12 -3.95 -3.79 -3.69 -3.51 3 1.68
## log_f[60] -3.84 0.10 0.17 -4.14 -3.98 -3.82 -3.72 -3.55 3 1.69
## log_f[61] -3.87 0.09 0.16 -4.17 -4.01 -3.85 -3.75 -3.58 3 1.70
## log_f[62] -3.90 0.09 0.16 -4.19 -4.04 -3.88 -3.79 -3.62 3 1.71
## log_f[63] -3.94 0.09 0.16 -4.22 -4.07 -3.92 -3.82 -3.66 3 1.72
## log_f[64] -3.98 0.09 0.16 -4.25 -4.11 -3.95 -3.86 -3.71 3 1.73
## log_f[65] -4.01 0.09 0.15 -4.29 -4.14 -3.99 -3.90 -3.75 3 1.74
## log_f[66] -4.06 0.09 0.15 -4.32 -4.18 -4.03 -3.95 -3.80 3 1.74
## log_f[67] -4.10 0.08 0.14 -4.35 -4.22 -4.08 -3.99 -3.85 3 1.74
## log_f[68] -4.14 0.08 0.14 -4.39 -4.26 -4.12 -4.04 -3.90 3 1.73
## log_f[69] -4.19 0.08 0.14 -4.43 -4.31 -4.17 -4.09 -3.95 3 1.72
## log_f[70] -4.24 0.08 0.13 -4.47 -4.35 -4.22 -4.14 -4.01 3 1.70
## log_f[71] -4.28 0.07 0.13 -4.52 -4.39 -4.27 -4.19 -4.06 3 1.66
## log_f[72] -4.34 0.07 0.12 -4.56 -4.44 -4.32 -4.25 -4.12 3 1.62
## log_f[73] -4.39 0.07 0.12 -4.61 -4.49 -4.38 -4.31 -4.18 3 1.56
## log_f[74] -4.44 0.06 0.12 -4.66 -4.54 -4.43 -4.36 -4.24 4 1.49
## log_f[75] -4.50 0.06 0.11 -4.71 -4.59 -4.49 -4.42 -4.30 4 1.41
## log_f[76] -4.56 0.05 0.11 -4.77 -4.64 -4.55 -4.48 -4.36 4 1.33
## log_f[77] -4.62 0.05 0.11 -4.83 -4.70 -4.62 -4.54 -4.42 5 1.26
## log_f[78] -4.68 0.04 0.11 -4.90 -4.76 -4.68 -4.60 -4.48 7 1.19
## log_f[79] -4.75 0.04 0.11 -4.97 -4.82 -4.75 -4.67 -4.54 9 1.13
## log_f[80] -4.82 0.03 0.12 -5.05 -4.89 -4.81 -4.73 -4.60 15 1.08
## log_f[81] -4.88 0.02 0.12 -5.14 -4.96 -4.88 -4.80 -4.66 31 1.05
## log_f[82] -4.95 0.01 0.13 -5.22 -5.04 -4.95 -4.86 -4.72 118 1.02
## log_f[83] -5.03 0.00 0.14 -5.32 -5.11 -5.02 -4.93 -4.79 2015 1.01
## log_f[84] -5.10 0.00 0.14 -5.42 -5.19 -5.09 -5.00 -4.85 4278 1.00
## log_f[85] -5.18 0.00 0.16 -5.52 -5.27 -5.16 -5.06 -4.91 4221 1.00
## log_f[86] -5.25 0.00 0.17 -5.63 -5.35 -5.24 -5.13 -4.97 3933 1.00
## log_f[87] -5.33 0.00 0.18 -5.73 -5.44 -5.31 -5.20 -5.03 2694 1.01
## log_f[88] -5.41 0.01 0.19 -5.85 -5.53 -5.40 -5.27 -5.09 518 1.01
## log_f[89] -5.50 0.02 0.21 -5.97 -5.63 -5.47 -5.34 -5.15 140 1.02
## log_f[90] -5.58 0.03 0.23 -6.09 -5.72 -5.56 -5.42 -5.21 72 1.03
## log_f[91] -5.67 0.04 0.24 -6.21 -5.82 -5.64 -5.49 -5.27 42 1.04
## log_f[92] -5.76 0.05 0.26 -6.34 -5.92 -5.73 -5.57 -5.33 30 1.05
## log_f[93] -5.85 0.06 0.28 -6.47 -6.02 -5.82 -5.64 -5.39 24 1.06
## log_f[94] -5.94 0.07 0.30 -6.60 -6.12 -5.91 -5.72 -5.45 20 1.07
## log_f[95] -6.03 0.08 0.32 -6.74 -6.23 -6.00 -5.81 -5.51 17 1.08
## log_f[96] -6.13 0.09 0.34 -6.87 -6.34 -6.09 -5.89 -5.57 15 1.09
## log_f[97] -6.23 0.10 0.36 -7.02 -6.45 -6.19 -5.97 -5.63 13 1.10
## log_f[98] -6.33 0.11 0.38 -7.16 -6.56 -6.29 -6.05 -5.70 12 1.10
## log_f[99] -6.43 0.12 0.40 -7.32 -6.68 -6.39 -6.14 -5.76 11 1.11
## log_f[100] -6.53 0.13 0.42 -7.47 -6.79 -6.49 -6.23 -5.83 11 1.12
## lp__ -429.78 3.77 5.55 -441.20 -435.17 -427.14 -425.78 -424.57 2 3.54
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 05:43:09 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit2i, 2)
tidybayes::tidy_draws(normal_fixed_mixture_stan_fit2i) %>%
select(.chain, .iteration, .draw, starts_with("log_f")) %>%
gather(key = key, value = value, starts_with("log_f")) %>%
mutate(key = gsub("log_f|\\[|\\]", "", key) %>% as.integer(),
x = factor(key, labels = seq(from = grid_min, to = grid_max, length.out = grid_length)) %>%
as.character() %>%
as.numeric(),
value = exp(value)) %>%
ggplot(mapping = aes(x = x, y = value, group = interaction(.chain, .iteration, .draw))) +
geom_line(size = 0.1, alpha = 1/20) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
legend.key = element_blank(),
plot.title = element_text(hjust = 0.5),
strip.background = element_blank())
2 latent cluster model with more precise prior on variance, concentrated probability (density grid)
normal_fixed_mixture_stan_fit2j <-
rstan::stan(model_code = normal_fixed_mixture_grid_stan_code,
data = list(alpha = 10^(1), beta = 10^(1),
m = 0, s_squared = 10,
n = nrow(galaxy),
y = galaxy$k_speed,
dirichlet_alpha = 2,
H = 2,
grid_max = grid_max, grid_min = grid_min, grid_length = grid_length))normal_fixed_mixture_stan_fit2j## Inference for Stan model: 963ec46299674d2582fc7e9c128240d0.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu[1] 11.30 2.11 3.14 8.95 9.45 9.74 10.94 19.02 2 3.26
## mu[2] 21.53 0.12 0.37 20.81 21.29 21.53 21.77 22.24 9 1.14
## sigma_squared[1] 10.65 12.03 17.73 0.52 0.75 0.94 6.16 53.64 2 3.58
## sigma_squared[2] 7.01 1.72 2.67 1.99 5.54 7.72 8.84 11.11 2 2.36
## Pi[1] 0.15 0.06 0.11 0.04 0.08 0.11 0.16 0.43 3 1.79
## Pi[2] 0.85 0.06 0.11 0.57 0.84 0.89 0.92 0.96 3 1.79
## sigma[1] 2.26 1.64 2.35 0.72 0.87 0.97 2.25 7.32 2 5.84
## sigma[2] 2.59 0.38 0.57 1.41 2.35 2.78 2.97 3.33 2 2.69
## log_f[1] -86.44 26.59 39.78 -139.39 -114.29 -100.02 -59.71 -15.33 2 2.93
## log_f[2] -84.00 25.83 38.64 -135.42 -111.06 -97.19 -58.00 -14.95 2 2.93
## log_f[3] -81.60 25.07 37.51 -131.50 -107.86 -94.40 -56.32 -14.58 2 2.93
## log_f[4] -79.23 24.33 36.39 -127.64 -104.71 -91.65 -54.66 -14.21 2 2.93
## log_f[5] -76.90 23.60 35.30 -123.83 -101.62 -88.93 -53.03 -13.85 2 2.93
## log_f[6] -74.60 22.87 34.22 -120.09 -98.57 -86.27 -51.43 -13.48 2 2.93
## log_f[7] -72.34 22.16 33.16 -116.40 -95.59 -83.67 -49.85 -13.12 2 2.93
## log_f[8] -70.11 21.46 32.11 -112.78 -92.65 -81.08 -48.29 -12.76 2 2.93
## log_f[9] -67.92 20.77 31.08 -109.21 -89.75 -78.53 -46.77 -12.40 2 2.93
## log_f[10] -65.77 20.10 30.06 -105.70 -86.89 -76.02 -45.27 -12.06 2 2.93
## log_f[11] -63.66 19.43 29.07 -102.25 -84.07 -73.58 -43.79 -11.71 2 2.93
## log_f[12] -61.57 18.77 28.08 -98.89 -81.29 -71.17 -42.34 -11.38 2 2.93
## log_f[13] -59.53 18.13 27.12 -95.59 -78.56 -68.79 -40.92 -11.08 2 2.93
## log_f[14] -57.52 17.49 26.17 -92.35 -75.90 -66.45 -39.52 -10.78 2 2.93
## log_f[15] -55.55 16.87 25.24 -89.14 -73.26 -64.17 -38.15 -10.51 2 2.93
## log_f[16] -53.61 16.25 24.32 -85.97 -70.70 -61.91 -36.80 -10.23 2 2.93
## log_f[17] -51.71 15.65 23.42 -82.85 -68.14 -59.70 -35.48 -9.96 2 2.92
## log_f[18] -49.85 15.06 22.54 -79.80 -65.66 -57.54 -34.18 -9.70 2 2.92
## log_f[19] -48.02 14.48 21.67 -76.84 -63.22 -55.41 -32.91 -9.44 2 2.92
## log_f[20] -46.23 13.91 20.82 -73.93 -60.82 -53.33 -31.67 -9.19 2 2.92
## log_f[21] -44.48 13.35 19.98 -71.05 -58.49 -51.28 -30.45 -8.94 2 2.92
## log_f[22] -42.76 12.80 19.16 -68.29 -56.19 -49.27 -29.26 -8.70 2 2.92
## log_f[23] -41.07 12.27 18.36 -65.56 -53.94 -47.32 -28.10 -8.50 2 2.92
## log_f[24] -39.42 11.74 17.57 -62.88 -51.74 -45.40 -26.96 -8.25 2 2.92
## log_f[25] -37.81 11.22 16.80 -60.26 -49.60 -43.53 -25.84 -8.03 2 2.92
## log_f[26] -36.23 10.72 16.05 -57.72 -47.49 -41.69 -24.76 -7.79 2 2.91
## log_f[27] -34.69 10.23 15.31 -55.20 -45.42 -39.89 -23.69 -7.57 2 2.91
## log_f[28] -33.19 9.74 14.58 -52.71 -43.39 -38.13 -22.66 -7.38 2 2.91
## log_f[29] -31.71 9.27 13.88 -50.30 -41.42 -36.42 -21.65 -7.19 2 2.91
## log_f[30] -30.28 8.80 13.18 -47.96 -39.49 -34.73 -20.66 -6.97 2 2.91
## log_f[31] -28.87 8.35 12.50 -45.62 -37.61 -33.09 -19.70 -6.78 2 2.91
## log_f[32] -27.50 7.91 11.84 -43.38 -35.78 -31.49 -18.77 -6.60 2 2.90
## log_f[33] -26.17 7.47 11.19 -41.21 -33.98 -29.94 -17.86 -6.42 2 2.90
## log_f[34] -24.86 7.04 10.55 -38.99 -32.25 -28.41 -16.98 -6.24 2 2.90
## log_f[35] -23.57 6.62 9.92 -36.89 -30.54 -26.91 -16.12 -6.08 2 2.89
## log_f[36] -22.29 6.20 9.29 -34.80 -28.87 -25.41 -15.17 -5.92 2 2.88
## log_f[37] -21.00 5.77 8.66 -32.75 -27.14 -23.92 -13.69 -5.77 2 2.85
## log_f[38] -19.69 5.33 8.01 -30.65 -25.35 -22.43 -12.36 -5.63 2 2.81
## log_f[39] -18.31 4.85 7.34 -28.41 -23.50 -20.79 -11.27 -5.49 2 2.70
## log_f[40] -16.82 4.33 6.63 -26.15 -21.66 -18.96 -9.82 -5.37 2 2.53
## log_f[41] -15.19 3.75 5.86 -23.63 -19.53 -16.92 -9.09 -5.26 2 2.28
## log_f[42] -13.39 3.10 5.02 -21.14 -17.25 -14.63 -8.40 -5.14 3 2.01
## log_f[43] -11.47 2.40 4.08 -18.35 -14.58 -12.05 -7.39 -5.04 3 1.75
## log_f[44] -9.55 1.68 3.09 -15.49 -11.78 -9.64 -6.62 -4.92 3 1.53
## log_f[45] -7.79 1.01 2.13 -12.41 -9.23 -7.57 -5.92 -4.78 4 1.33
## log_f[46] -6.30 0.42 1.33 -9.47 -7.09 -6.02 -5.27 -4.45 10 1.13
## log_f[47] -5.14 0.02 0.82 -7.05 -5.57 -5.06 -4.59 -3.75 2794 1.00
## log_f[48] -4.31 0.30 0.68 -5.66 -4.82 -4.27 -3.77 -3.15 5 1.26
## log_f[49] -3.83 0.44 0.73 -5.33 -4.40 -3.64 -3.26 -2.79 3 1.86
## log_f[50] -3.68 0.45 0.72 -5.17 -4.23 -3.47 -3.15 -2.65 3 2.01
## log_f[51] -3.87 0.33 0.62 -5.06 -4.37 -3.78 -3.40 -2.86 4 1.47
## log_f[52] -4.40 0.01 0.59 -5.68 -4.74 -4.39 -4.01 -3.30 1603 1.01
## log_f[53] -5.19 0.31 0.87 -7.13 -5.76 -5.04 -4.51 -3.94 8 1.17
## log_f[54] -5.91 0.66 1.17 -8.07 -6.77 -6.06 -4.80 -4.02 3 1.61
## log_f[55] -6.06 0.74 1.21 -8.04 -6.93 -6.30 -4.98 -3.94 3 1.95
## log_f[56] -5.74 0.64 1.05 -7.56 -6.49 -5.93 -4.86 -3.86 3 1.85
## log_f[57] -5.31 0.48 0.84 -6.85 -5.91 -5.40 -4.65 -3.79 3 1.65
## log_f[58] -4.89 0.33 0.64 -6.14 -5.33 -4.91 -4.43 -3.71 4 1.42
## log_f[59] -4.50 0.18 0.48 -5.49 -4.81 -4.48 -4.17 -3.63 7 1.19
## log_f[60] -4.14 0.04 0.36 -4.91 -4.37 -4.12 -3.90 -3.49 78 1.03
## log_f[61] -3.80 0.01 0.30 -4.44 -3.99 -3.79 -3.59 -3.25 1032 1.01
## log_f[62] -3.48 0.09 0.28 -4.08 -3.65 -3.45 -3.28 -3.01 11 1.12
## log_f[63] -3.16 0.10 0.25 -3.72 -3.30 -3.13 -2.98 -2.77 7 1.20
## log_f[64] -2.86 0.06 0.20 -3.33 -2.96 -2.83 -2.72 -2.55 10 1.13
## log_f[65] -2.58 0.00 0.15 -2.91 -2.66 -2.57 -2.49 -2.33 3596 1.00
## log_f[66] -2.35 0.05 0.13 -2.60 -2.43 -2.36 -2.28 -2.04 8 1.15
## log_f[67] -2.17 0.09 0.16 -2.42 -2.28 -2.21 -2.10 -1.79 3 1.67
## log_f[68] -2.05 0.12 0.18 -2.31 -2.18 -2.11 -1.98 -1.63 3 2.11
## log_f[69] -2.00 0.12 0.19 -2.25 -2.13 -2.05 -1.92 -1.58 3 2.12
## log_f[70] -2.01 0.10 0.17 -2.25 -2.12 -2.05 -1.94 -1.63 3 1.68
## log_f[71] -2.08 0.05 0.14 -2.31 -2.17 -2.09 -2.02 -1.79 8 1.16
## log_f[72] -2.22 0.01 0.14 -2.54 -2.28 -2.21 -2.14 -2.00 645 1.02
## log_f[73] -2.41 0.11 0.21 -3.03 -2.47 -2.36 -2.28 -2.15 4 1.40
## log_f[74] -2.66 0.20 0.33 -3.53 -2.73 -2.54 -2.45 -2.30 3 1.90
## log_f[75] -2.95 0.29 0.45 -4.08 -3.05 -2.77 -2.66 -2.48 2 2.31
## log_f[76] -3.27 0.36 0.56 -4.64 -3.43 -3.04 -2.90 -2.69 2 2.53
## log_f[77] -3.62 0.41 0.64 -5.16 -3.81 -3.37 -3.19 -2.93 2 2.48
## log_f[78] -3.98 0.43 0.68 -5.64 -4.22 -3.73 -3.52 -3.20 2 2.22
## log_f[79] -4.36 0.41 0.69 -6.06 -4.60 -4.14 -3.89 -3.51 3 1.88
## log_f[80] -4.76 0.37 0.67 -6.43 -5.06 -4.60 -4.29 -3.85 3 1.56
## log_f[81] -5.18 0.29 0.65 -6.76 -5.52 -5.07 -4.73 -4.20 5 1.30
## log_f[82] -5.64 0.18 0.63 -7.06 -6.01 -5.56 -5.19 -4.59 12 1.12
## log_f[83] -6.12 0.02 0.65 -7.48 -6.55 -6.08 -5.66 -4.94 892 1.02
## log_f[84] -6.64 0.02 0.72 -8.08 -7.11 -6.62 -6.15 -5.31 1748 1.00
## log_f[85] -7.20 0.15 0.83 -8.85 -7.74 -7.19 -6.63 -5.62 30 1.05
## log_f[86] -7.79 0.33 0.99 -9.76 -8.45 -7.78 -7.14 -5.88 9 1.14
## log_f[87] -8.42 0.50 1.18 -10.70 -9.21 -8.40 -7.63 -6.13 6 1.24
## log_f[88] -9.08 0.68 1.40 -11.73 -10.03 -9.11 -8.12 -6.35 4 1.35
## log_f[89] -9.78 0.86 1.65 -12.81 -10.93 -9.88 -8.59 -6.60 4 1.45
## log_f[90] -10.51 1.05 1.92 -13.96 -11.86 -10.70 -9.09 -6.84 3 1.54
## log_f[91] -11.29 1.26 2.21 -15.16 -12.85 -11.56 -9.65 -7.08 3 1.63
## log_f[92] -12.09 1.47 2.52 -16.42 -13.89 -12.47 -10.20 -7.32 3 1.71
## log_f[93] -12.94 1.69 2.85 -17.73 -14.98 -13.41 -10.83 -7.57 3 1.78
## log_f[94] -13.82 1.93 3.20 -19.12 -16.11 -14.40 -11.40 -7.82 3 1.85
## log_f[95] -14.73 2.17 3.56 -20.55 -17.28 -15.43 -12.02 -8.09 3 1.91
## log_f[96] -15.68 2.43 3.94 -22.04 -18.51 -16.51 -12.66 -8.35 3 1.96
## log_f[97] -16.67 2.70 4.34 -23.57 -19.79 -17.63 -13.40 -8.61 3 2.01
## log_f[98] -17.69 2.98 4.76 -25.16 -21.11 -18.78 -14.12 -8.88 3 2.06
## log_f[99] -18.75 3.27 5.19 -26.83 -22.49 -19.98 -14.80 -9.17 3 2.10
## log_f[100] -19.85 3.57 5.64 -28.55 -23.90 -21.22 -15.41 -9.46 3 2.14
## lp__ -293.85 8.55 12.19 -317.32 -301.47 -287.27 -285.94 -284.78 2 7.77
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 06:06:36 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit2j, 2)
tidybayes::tidy_draws(normal_fixed_mixture_stan_fit2j) %>%
select(.chain, .iteration, .draw, starts_with("log_f")) %>%
gather(key = key, value = value, starts_with("log_f")) %>%
mutate(key = gsub("log_f|\\[|\\]", "", key) %>% as.integer(),
x = factor(key, labels = seq(from = grid_min, to = grid_max, length.out = grid_length)) %>%
as.character() %>%
as.numeric(),
value = exp(value)) %>%
ggplot(mapping = aes(x = x, y = value, color = factor(.chain),
group = interaction(.chain, .iteration, .draw))) +
geom_line(size = 0.1, alpha = 1/20) +
facet_wrap(~ .chain) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
legend.key = element_blank(),
plot.title = element_text(hjust = 0.5),
strip.background = element_blank())
5 latent cluster model with more precise prior on variance (density grid)
normal_fixed_mixture_stan_fit5 <-
rstan::stan(model_code = normal_fixed_mixture_grid_stan_code,
data = list(alpha = 10^(1), beta = 10^(1),
m = 0, s_squared = 10,
n = nrow(galaxy),
y = galaxy$k_speed,
dirichlet_alpha = 1,
H = 5,
grid_max = grid_max, grid_min = grid_min, grid_length = grid_length))## Warning: There were 11 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Warning: Examine the pairs() plot to diagnose sampling problemsnormal_fixed_mixture_stan_fit5## Inference for Stan model: 963ec46299674d2582fc7e9c128240d0.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu[1] -2.49 0.10 2.45 -7.30 -4.11 -2.46 -0.86 2.06 603 1.02
## mu[2] 0.39 0.46 2.29 -3.83 -1.25 0.33 1.89 5.20 25 1.06
## mu[3] 4.39 2.18 3.63 -1.50 1.59 3.63 8.83 10.05 3 1.86
## mu[4] 12.42 3.46 4.92 8.94 9.44 9.75 11.78 21.59 2 10.51
## mu[5] 21.56 0.05 0.40 20.76 21.31 21.57 21.82 22.31 62 1.04
## sigma_squared[1] 1.11 0.01 0.39 0.58 0.83 1.03 1.29 2.07 3207 1.00
## sigma_squared[2] 1.10 0.01 0.39 0.59 0.83 1.03 1.28 2.10 2993 1.00
## sigma_squared[3] 1.05 0.05 0.37 0.55 0.80 0.98 1.22 1.99 58 1.03
## sigma_squared[4] 3.93 3.81 6.36 0.52 0.75 0.94 3.16 19.77 3 1.90
## sigma_squared[5] 6.70 2.07 3.11 0.84 5.13 7.75 8.80 10.90 2 2.84
## Pi[1] 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.02 4684 1.00
## Pi[2] 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.02 4257 1.00
## Pi[3] 0.02 0.03 0.04 0.00 0.00 0.00 0.03 0.13 2 2.48
## Pi[4] 0.22 0.16 0.26 0.04 0.07 0.10 0.17 0.91 3 2.16
## Pi[5] 0.75 0.19 0.29 0.00 0.81 0.90 0.92 0.96 2 2.42
## sigma[1] 1.04 0.00 0.17 0.76 0.91 1.02 1.14 1.44 3563 1.00
## sigma[2] 1.03 0.00 0.17 0.77 0.91 1.01 1.13 1.45 3411 1.00
## sigma[3] 1.01 0.02 0.17 0.74 0.89 0.99 1.11 1.41 47 1.04
## sigma[4] 1.58 0.80 1.20 0.72 0.87 0.97 1.75 4.45 2 3.10
## sigma[5] 2.48 0.50 0.74 0.92 2.26 2.78 2.97 3.30 2 3.31
## log_f[1] -95.78 7.77 22.72 -134.65 -110.62 -98.56 -83.91 -44.15 9 1.16
## log_f[2] -92.26 7.26 22.12 -129.97 -106.89 -94.96 -80.32 -42.34 9 1.14
## log_f[3] -88.70 6.73 21.55 -125.24 -103.03 -91.18 -76.73 -40.58 10 1.13
## log_f[4] -85.12 6.17 21.03 -121.17 -99.31 -87.58 -73.14 -39.06 12 1.11
## log_f[5] -81.50 5.59 20.53 -116.67 -95.66 -83.91 -69.34 -37.00 13 1.10
## log_f[6] -77.85 4.98 20.06 -112.61 -91.84 -80.39 -65.72 -34.36 16 1.09
## log_f[7] -74.16 4.33 19.61 -108.39 -88.10 -76.65 -61.31 -32.32 20 1.07
## log_f[8] -70.43 3.65 19.16 -104.17 -84.25 -72.80 -57.24 -30.68 27 1.06
## log_f[9] -66.68 2.94 18.72 -99.75 -80.44 -68.71 -53.37 -28.44 41 1.05
## log_f[10] -62.91 2.16 18.30 -95.91 -76.07 -64.80 -49.36 -25.85 72 1.04
## log_f[11] -59.13 1.44 17.89 -92.00 -72.27 -60.45 -45.66 -23.83 154 1.03
## log_f[12] -55.36 0.78 17.47 -87.89 -68.34 -56.19 -42.13 -21.86 498 1.02
## log_f[13] -51.61 0.57 17.02 -83.71 -64.32 -52.08 -38.56 -19.88 884 1.01
## log_f[14] -47.91 0.40 16.55 -80.05 -60.06 -47.63 -35.24 -17.77 1713 1.01
## log_f[15] -44.26 0.34 16.01 -75.52 -56.10 -43.59 -32.05 -15.53 2188 1.00
## log_f[16] -40.70 0.32 15.42 -71.23 -51.65 -39.55 -29.09 -13.66 2327 1.00
## log_f[17] -37.24 0.30 14.77 -67.65 -47.46 -35.73 -26.20 -11.94 2433 1.00
## log_f[18] -33.91 0.28 14.06 -63.58 -43.24 -32.20 -23.37 -10.45 2450 1.00
## log_f[19] -30.75 0.27 13.29 -59.61 -39.34 -28.90 -20.81 -9.31 2395 1.00
## log_f[20] -27.77 0.26 12.48 -55.75 -35.61 -25.82 -18.37 -8.33 2311 1.00
## log_f[21] -25.01 0.25 11.62 -51.71 -32.07 -22.95 -16.14 -7.68 2210 1.00
## log_f[22] -22.47 0.24 10.73 -48.13 -28.79 -20.41 -14.14 -7.24 2084 1.00
## log_f[23] -20.16 0.22 9.82 -44.57 -25.92 -18.12 -12.54 -6.74 1947 1.01
## log_f[24] -18.10 0.21 8.90 -40.45 -23.21 -16.12 -11.28 -6.46 1822 1.01
## log_f[25] -16.30 0.19 8.00 -36.91 -20.63 -14.45 -10.27 -6.17 1725 1.01
## log_f[26] -14.74 0.18 7.13 -33.22 -18.36 -12.97 -9.41 -5.90 1653 1.01
## log_f[27] -13.43 0.16 6.31 -29.47 -16.50 -11.89 -8.78 -5.59 1509 1.01
## log_f[28] -12.34 0.15 5.57 -26.87 -14.94 -10.95 -8.28 -5.41 1413 1.01
## log_f[29] -11.47 0.13 4.95 -24.26 -13.77 -10.25 -7.90 -5.30 1444 1.01
## log_f[30] -10.79 0.12 4.43 -22.32 -12.89 -9.69 -7.60 -5.26 1428 1.01
## log_f[31] -10.28 0.11 4.01 -20.66 -12.16 -9.30 -7.41 -5.24 1446 1.01
## log_f[32] -9.94 0.10 3.71 -19.25 -11.67 -9.08 -7.34 -5.14 1505 1.01
## log_f[33] -9.73 0.09 3.50 -18.37 -11.44 -8.96 -7.22 -5.10 1456 1.01
## log_f[34] -9.65 0.11 3.38 -18.02 -11.34 -8.94 -7.19 -5.11 935 1.01
## log_f[35] -9.69 0.10 3.38 -17.79 -11.47 -8.95 -7.24 -5.12 1057 1.01
## log_f[36] -9.86 0.09 3.49 -18.34 -11.69 -9.06 -7.28 -5.20 1677 1.01
## log_f[37] -10.16 0.08 3.67 -19.05 -12.15 -9.30 -7.45 -5.29 1871 1.01
## log_f[38] -10.57 0.08 3.90 -20.20 -12.68 -9.67 -7.63 -5.36 2177 1.01
## log_f[39] -11.07 0.08 4.12 -21.07 -13.39 -10.20 -7.86 -5.40 2651 1.00
## log_f[40] -11.55 0.08 4.17 -20.62 -14.25 -10.73 -8.26 -5.56 2896 1.00
## log_f[41] -11.90 0.07 3.98 -20.09 -14.75 -11.37 -8.72 -5.66 3058 1.00
## log_f[42] -11.87 0.06 3.54 -18.94 -14.39 -11.67 -9.20 -5.83 2993 1.00
## log_f[43] -11.25 0.06 2.92 -17.19 -13.17 -11.07 -9.22 -5.94 2613 1.00
## log_f[44] -10.00 0.05 2.30 -14.82 -11.50 -9.87 -8.34 -5.94 2539 1.01
## log_f[45] -8.33 0.03 1.77 -12.21 -9.42 -8.18 -7.05 -5.32 3156 1.00
## log_f[46] -6.64 0.02 1.30 -9.66 -7.42 -6.52 -5.71 -4.46 3683 1.00
## log_f[47] -5.22 0.01 0.91 -7.25 -5.77 -5.13 -4.57 -3.70 3850 1.00
## log_f[48] -4.18 0.01 0.61 -5.54 -4.57 -4.13 -3.75 -3.16 3610 1.00
## log_f[49] -3.59 0.01 0.44 -4.52 -3.86 -3.55 -3.28 -2.80 3161 1.00
## log_f[50] -3.43 0.01 0.41 -4.32 -3.70 -3.41 -3.14 -2.70 3088 1.00
## log_f[51] -3.72 0.01 0.48 -4.78 -4.03 -3.69 -3.38 -2.89 3744 1.00
## log_f[52] -4.43 0.01 0.65 -5.87 -4.84 -4.39 -3.98 -3.31 4603 1.00
## log_f[53] -5.45 0.02 0.83 -7.25 -5.99 -5.40 -4.86 -3.99 1452 1.01
## log_f[54] -6.31 0.20 0.89 -8.13 -6.92 -6.28 -5.67 -4.68 20 1.07
## log_f[55] -6.44 0.28 0.84 -8.09 -7.00 -6.45 -5.87 -4.83 9 1.14
## log_f[56] -6.05 0.26 0.76 -7.56 -6.54 -6.05 -5.54 -4.60 9 1.14
## log_f[57] -5.53 0.20 0.65 -6.83 -5.95 -5.52 -5.08 -4.33 11 1.12
## log_f[58] -5.03 0.14 0.54 -6.15 -5.37 -5.00 -4.65 -4.04 15 1.09
## log_f[59] -4.57 0.08 0.45 -5.50 -4.85 -4.55 -4.26 -3.75 28 1.05
## log_f[60] -4.14 0.03 0.37 -4.95 -4.38 -4.12 -3.89 -3.47 197 1.02
## log_f[61] -3.76 0.01 0.31 -4.44 -3.96 -3.74 -3.54 -3.21 2990 1.00
## log_f[62] -3.41 0.00 0.26 -4.01 -3.57 -3.39 -3.24 -2.96 2850 1.00
## log_f[63] -3.10 0.01 0.22 -3.60 -3.23 -3.08 -2.95 -2.74 1217 1.01
## log_f[64] -2.82 0.00 0.17 -3.20 -2.92 -2.81 -2.70 -2.53 1290 1.01
## log_f[65] -2.57 0.00 0.13 -2.86 -2.65 -2.56 -2.48 -2.34 3724 1.00
## log_f[66] -2.36 0.03 0.12 -2.58 -2.43 -2.36 -2.29 -2.11 20 1.07
## log_f[67] -2.20 0.06 0.14 -2.41 -2.28 -2.22 -2.14 -1.83 6 1.24
## log_f[68] -2.09 0.07 0.15 -2.30 -2.18 -2.12 -2.04 -1.67 4 1.35
## log_f[69] -2.04 0.07 0.15 -2.25 -2.13 -2.07 -1.99 -1.63 4 1.35
## log_f[70] -2.04 0.05 0.13 -2.26 -2.13 -2.06 -1.99 -1.71 7 1.19
## log_f[71] -2.11 0.00 0.12 -2.35 -2.18 -2.11 -2.05 -1.88 2624 1.01
## log_f[72] -2.23 0.04 0.14 -2.61 -2.29 -2.21 -2.15 -2.03 15 1.09
## log_f[73] -2.40 0.09 0.20 -2.96 -2.45 -2.35 -2.28 -2.15 5 1.30
## log_f[74] -2.60 0.13 0.26 -3.33 -2.65 -2.53 -2.44 -2.30 4 1.42
## log_f[75] -2.84 0.16 0.31 -3.69 -2.90 -2.76 -2.65 -2.49 4 1.45
## log_f[76] -3.11 0.17 0.34 -4.04 -3.20 -3.02 -2.89 -2.70 4 1.42
## log_f[77] -3.42 0.17 0.36 -4.39 -3.54 -3.34 -3.18 -2.95 4 1.35
## log_f[78] -3.76 0.16 0.38 -4.73 -3.92 -3.70 -3.51 -3.22 5 1.27
## log_f[79] -4.14 0.15 0.39 -5.10 -4.34 -4.10 -3.87 -3.51 7 1.19
## log_f[80] -4.57 0.12 0.41 -5.49 -4.80 -4.53 -4.29 -3.86 12 1.11
## log_f[81] -5.03 0.08 0.45 -5.99 -5.31 -5.00 -4.73 -4.24 30 1.05
## log_f[82] -5.53 0.02 0.50 -6.60 -5.85 -5.51 -5.19 -4.65 464 1.02
## log_f[83] -6.08 0.01 0.57 -7.32 -6.44 -6.05 -5.68 -5.06 3712 1.00
## log_f[84] -6.66 0.01 0.66 -8.09 -7.09 -6.63 -6.21 -5.50 3780 1.00
## log_f[85] -7.29 0.02 0.77 -8.93 -7.79 -7.26 -6.76 -5.92 1857 1.01
## log_f[86] -7.96 0.07 0.89 -9.82 -8.54 -7.93 -7.36 -6.32 179 1.02
## log_f[87] -8.68 0.13 1.02 -10.78 -9.34 -8.65 -7.97 -6.73 63 1.03
## log_f[88] -9.43 0.19 1.17 -11.81 -10.18 -9.42 -8.64 -7.18 37 1.04
## log_f[89] -10.23 0.26 1.33 -12.89 -11.09 -10.23 -9.35 -7.60 26 1.06
## log_f[90] -11.06 0.33 1.50 -14.01 -12.04 -11.08 -10.09 -8.02 21 1.07
## log_f[91] -11.94 0.40 1.68 -15.21 -13.04 -11.98 -10.88 -8.51 18 1.08
## log_f[92] -12.87 0.47 1.87 -16.44 -14.08 -12.91 -11.71 -8.96 16 1.08
## log_f[93] -13.83 0.55 2.08 -17.75 -15.17 -13.88 -12.56 -9.41 14 1.09
## log_f[94] -14.84 0.63 2.29 -19.13 -16.31 -14.89 -13.46 -9.90 13 1.10
## log_f[95] -15.88 0.71 2.51 -20.54 -17.50 -15.96 -14.40 -10.42 12 1.11
## log_f[96] -16.97 0.80 2.74 -22.04 -18.75 -17.07 -15.35 -10.94 12 1.11
## log_f[97] -18.10 0.89 2.99 -23.59 -20.04 -18.23 -16.35 -11.49 11 1.12
## log_f[98] -19.28 0.99 3.24 -25.17 -21.38 -19.42 -17.39 -12.03 11 1.12
## log_f[99] -20.49 1.09 3.50 -26.85 -22.75 -20.66 -18.49 -12.64 10 1.13
## log_f[100] -21.75 1.19 3.77 -28.56 -24.18 -21.94 -19.61 -13.24 10 1.13
## lp__ -327.20 7.21 10.55 -348.59 -335.42 -322.35 -319.81 -316.75 2 3.74
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 30 06:08:08 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot_all(normal_fixed_mixture_stan_fit5, 5)
tidybayes::tidy_draws(normal_fixed_mixture_stan_fit5) %>%
select(.chain, .iteration, .draw, starts_with("log_f")) %>%
gather(key = key, value = value, starts_with("log_f")) %>%
mutate(key = gsub("log_f|\\[|\\]", "", key) %>% as.integer(),
x = factor(key, labels = seq(from = grid_min, to = grid_max, length.out = grid_length)) %>%
as.character() %>%
as.numeric(),
value = exp(value)) %>%
ggplot(mapping = aes(x = x, y = value, color = factor(.chain),
group = interaction(.chain, .iteration, .draw))) +
geom_line(size = 0.1, alpha = 1/20) +
facet_wrap(~ .chain) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
legend.key = element_blank(),
plot.title = element_text(hjust = 0.5),
strip.background = element_blank())
- Top Page: http://rpubs.com/kaz_yos/
- Github: https://github.com/kaz-yos
print(sessionInfo())## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.14.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2.2 DPpackage_1.1-7.4 rstan_2.18.1 StanHeaders_2.18.0 forcats_0.3.0
## [6] stringr_1.3.1 dplyr_0.7.7 purrr_0.2.5 readr_1.1.1 tidyr_0.8.2
## [11] tibble_1.4.2 ggplot2_3.1.0 tidyverse_1.2.1 doRNG_1.7.1 rngtools_1.3.1
## [16] pkgmaker_0.27 registry_0.5 doParallel_1.0.14 iterators_1.0.10 foreach_1.4.4
## [21] knitr_1.20
##
## loaded via a namespace (and not attached):
## [1] minqa_1.2.4 colorspace_1.3-2 ggridges_0.5.1 rsconnect_0.8.8
## [5] rprojroot_1.3-2 ggstance_0.3.1 markdown_0.8 base64enc_0.1-3
## [9] rstudioapi_0.8 svUnit_0.7-12 DT_0.4 fansi_0.4.0
## [13] lubridate_1.7.4 xml2_1.2.0 codetools_0.2-15 splines_3.5.1
## [17] shinythemes_1.1.1 bayesplot_1.6.0 jsonlite_1.5 nloptr_1.2.1
## [21] broom_0.5.0 shiny_1.1.0 compiler_3.5.1 httr_1.3.1
## [25] backports_1.1.2 assertthat_0.2.0 Matrix_1.2-14 lazyeval_0.2.1
## [29] cli_1.0.1 later_0.7.5 htmltools_0.3.6 prettyunits_1.0.2
## [33] tools_3.5.1 igraph_1.2.2 coda_0.19-2 gtable_0.2.0
## [37] glue_1.3.0 reshape2_1.4.3 Rcpp_0.12.19 cellranger_1.1.0
## [41] nlme_3.1-137 crosstalk_1.0.0 ps_1.2.0 lme4_1.1-18-1
## [45] rvest_0.3.2 mime_0.6 miniUI_0.1.1.1 tidybayes_1.0.3
## [49] gtools_3.8.1 MASS_7.3-51 zoo_1.8-4 scales_1.0.0
## [53] rstanarm_2.18.1 colourpicker_1.0 hms_0.4.2 promises_1.0.1
## [57] inline_0.3.15 shinystan_2.5.0 yaml_2.2.0 gridExtra_2.3
## [61] loo_2.0.0 stringi_1.2.4 dygraphs_1.1.1.6 pkgbuild_1.0.2
## [65] bibtex_0.4.2 rlang_0.3.0.1 pkgconfig_2.0.2 matrixStats_0.54.0
## [69] evaluate_0.12 lattice_0.20-35 splines2_0.2.8 bindr_0.1.1
## [73] rstantools_1.5.1 htmlwidgets_1.3 labeling_0.3 processx_3.2.0
## [77] tidyselect_0.2.5 biostan_0.1.0 plyr_1.8.4 magrittr_1.5
## [81] R6_2.3.0 pillar_1.3.0 haven_1.1.2 withr_2.1.2
## [85] xts_0.11-1 survival_2.43-1 modelr_0.1.2 crayon_1.3.4
## [89] arrayhelpers_1.0-20160527 utf8_1.1.4 rmarkdown_1.10 grid_3.5.1
## [93] readxl_1.1.0 callr_3.0.0 threejs_0.3.1 digest_0.6.18
## [97] xtable_1.8-3 httpuv_1.4.5 stats4_3.5.1 munsell_0.5.0
## [101] shinyjs_1.0## Record execution time and multicore use
end_time <- Sys.time()
diff_time <- difftime(end_time, start_time, units = "auto")
cat("Started ", as.character(start_time), "\n",
"Finished ", as.character(end_time), "\n",
"Time difference of ", diff_time, " ", attr(diff_time, "units"), "\n",
"Used ", foreach::getDoParWorkers(), " cores\n",
"Used ", foreach::getDoParName(), " as backend\n",
sep = "")## Started 2018-11-30 06:02:39
## Finished 2018-11-30 06:09:28
## Time difference of 6.819778 mins
## Used 12 cores
## Used doParallelMC as backend