bmixt {bmixture} | R Documentation |
Sampling algorithm for mixture of t-distributions
Description
This function consists of several sampling algorithms for Bayesian estimation for finite mixture of t-distributions.
Usage
bmixt( data, k = "unknown", iter = 1000, burnin = iter / 2, lambda = 1, df = 1,
k.start = NULL, mu.start = NULL, sig.start = NULL, pi.start = NULL,
k.max = 30, trace = TRUE )
Arguments
data |
vector of data with size |
k |
number of components of mixture distribution. Defult is |
iter |
number of iteration for the sampling algorithm. |
burnin |
number of burn-in iteration for the sampling algorithm. |
lambda |
For the case |
df |
Degrees of freedom (> 0, maybe non-integer). df = Inf is allowed. |
k.start |
For the case |
mu.start |
Initial value for parameter of mixture distribution. |
sig.start |
Initial value for parameter of mixture distribution. |
pi.start |
Initial value for parameter of mixture distribution. |
k.max |
For the case |
trace |
Logical: if TRUE (default), tracing information is printed. |
Details
Sampling from finite mixture of t-distribution, with density:
Pr(x|k, \underline{\pi}, \underline{\mu}, \underline{\sigma}) = \sum_{i=1}^{k} \pi_{i} t_p(x|\mu_{i}, \sigma^2_{i}),
where k
is the number of components of mixture distribution (as a defult we assume is unknown
).
The prior distributions are defined as below
P(K=k) \propto \frac{\lambda^k}{k!}, \ \ \ k=1,...,k_{max},
\pi_{i} | k \sim Dirichlet( 1,..., 1 ),
\mu_{i} | k \sim N( \epsilon, \kappa ),
\sigma_i | k \sim IG( g, h ),
where IG
denotes an inverted gamma distribution. For more details see Stephens, M. (2000), doi: 10.1214/aos/1016120364.
Value
An object with S3
class "bmixt"
is returned:
all_k |
a vector which includes the waiting times for all iterations. It is needed for monitoring the convergence of the BD-MCMC algorithm. |
all_weights |
a vector which includes the waiting times for all iterations. It is needed for monitoring the convergence of the BD-MCMC algorithm. |
pi_sample |
a vector which includes the MCMC samples after burn-in from parameter |
mu_sample |
a vector which includes the MCMC samples after burn-in from parameter |
sig_sample |
a vector which includes the MCMC samples after burn-in from parameter |
data |
original data. |
Author(s)
Reza Mohammadi a.mohammadi@uva.nl
References
Stephens, M. (2000) Bayesian analysis of mixture models with an unknown number of components-an alternative to reversible jump methods. Annals of statistics, 28(1):40-74, doi: 10.1214/aos/1016120364
Richardson, S. and Green, P. J. (1997) On Bayesian analysis of mixtures with an unknown number of components. Journal of the Royal Statistical Society: series B, 59(4):731-792, doi: 10.1111/1467-9868.00095
Green, P. J. (1995) Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika, 82(4):711-732, doi: 10.1093/biomet/82.4.711
Cappe, O., Christian P. R., and Tobias, R. (2003) Reversible jump, birth and death and more general continuous time Markov chain Monte Carlo samplers. Journal of the Royal Statistical Society: Series B, 65(3):679-700
Mohammadi, A., Salehi-Rad, M. R., and Wit, E. C. (2013) Using mixture of Gamma distributions for Bayesian analysis in an M/G/1 queue with optional second service. Computational Statistics, 28(2):683-700, doi: 10.1007/s00180-012-0323-3
Mohammadi, A., and Salehi-Rad, M. R. (2012) Bayesian inference and prediction in an M/G/1 with optional second service. Communications in Statistics-Simulation and Computation, 41(3):419-435, doi: 10.1080/03610918.2011.588358
Wade, S. and Ghahramani, Z. (2018) Bayesian Cluster Analysis: Point Estimation and Credible Balls (with Discussion). Bayesian Analysis, 13(2):559-626, doi: 10.1214/17-BA1073
See Also
Examples
## Not run:
set.seed( 20 )
# simulating data from mixture of Normal with 3 components
n = 2000
weight = c( 0.3, 0.5, 0.2 )
mean = c( 0 , 10 , 3 )
sd = c( 1 , 1 , 1 )
data = rmixnorm( n = n, weight = weight, mean = mean, sd = sd )
# plot for simulation data
hist( data, prob = TRUE, nclass = 30, col = "gray" )
x = seq( -20, 20, 0.05 )
densmixnorm = dmixnorm( x, weight, mean, sd )
lines( x, densmixnorm, lwd = 2 )
# Runing bdmcmc algorithm for the above simulation data set
bmixt.obj = bmixt( data, k = 3, iter = 5000 )
summary( bmixt.obj )
## End(Not run)