bmixgamma function

Sampling algorithm for mixture of gamma distributions

Sampling algorithm for mixture of gamma distributions

This function consists of several sampling algorithms for Bayesian estimation for a mixture of Gamma distributions.

bmixgamma( data, k = "unknown", iter = 1000, burnin = iter / 2, lambda = 1, mu = NULL, nu = NULL, kesi = NULL, tau = NULL, k.start = NULL, alpha.start = NULL, beta.start = NULL, pi.start = NULL, k.max = 30, trace = TRUE )

Arguments

  • data: vector of data with size n.
  • k: number of components of mixture distribution. It can take an integer values.
  • iter: number of iteration for the sampling algorithm.
  • burnin: number of burn-in iteration for the sampling algorithm.
  • lambda: For the case k = "unknown", it is the parameter of the prior distribution of number of components k.
  • mu: parameter of alpha in mixture distribution.
  • nu: parameter of alpha in mixture distribution.
  • kesi: parameter of beta in mixture distribution.
  • tau: parameter of beta in mixture distribution.
  • k.start: For the case k = "unknown", initial value for number of components of mixture distribution.
  • alpha.start: Initial value for parameter of mixture distribution.
  • beta.start: Initial value for parameter of mixture distribution.
  • pi.start: Initial value for parameter of mixture distribution.
  • k.max: For the case k = "unknown", maximum value for the number of components of mixture distribution.
  • trace: Logical: if TRUE (default), tracing information is printed.

Details

Sampling from finite mixture of Gamma distribution, with density:

Pr(xk,π,α,β)=i=1kπiGamma(xαi,βi), Pr(x|k, \underline{\pi}, \underline{\alpha}, \underline{\beta}) = \sum_{i=1}^{k} \pi_{i} Gamma(x|\alpha_{i}, \beta_{i}),

where k is the number of components of mixture distribution (as a defult we assume is unknown) and

Gamma(xαi,βi)=(βi)αiΓ(αi)xαi1eβix. Gamma(x|\alpha_{i}, \beta_{i})=\frac{(\beta_{i})^{\alpha_{i}}}{\Gamma(\alpha_{i})} x^{\alpha_{i}-1} e^{-\beta_{i}x}.

The prior distributions are defined as below

P(K=k)λkk!,   k=1,...,kmax, P(K=k) \propto \frac{\lambda^k}{k!}, \ \ \ k=1,...,k_{max}, πikDirichlet(1,...,1), \pi_{i} | k \sim Dirichlet( 1,..., 1 ), αikGamma(ν,υ), \alpha_{i} | k \sim Gamma(\nu, \upsilon), βikGamma(η,τ), \beta_i | k \sim Gamma(\eta, \tau),

for more details see Mohammadi et al. (2013), tools:::Rd_expr_doi("10.1007/s00180-012-0323-3") .

Returns

An object with S3 class "bmixgamma" 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 pi of mixture distribution.

  • alpha_sample: a vector which includes the MCMC samples after burn-in from parameter alpha of mixture distribution.

  • beta_sample: a vector which includes the MCMC samples after burn-in from parameter beta of mixture distribution.

  • data: original data.

References

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, tools:::Rd_expr_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, tools:::Rd_expr_doi("10.1080/03610918.2011.588358")

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, tools:::Rd_expr_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, tools:::Rd_expr_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, tools:::Rd_expr_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

Wade, S. and Ghahramani, Z. (2018) Bayesian Cluster Analysis: Point Estimation and Credible Balls (with Discussion). Bayesian Analysis, 13(2):559-626, tools:::Rd_expr_doi("10.1214/17-BA1073")

Author(s)

Reza Mohammadi a.mohammadi@uva.nl

See Also

bmixnorm, bmixt, bmixgamma

Examples

## Not run: set.seed( 70 ) # simulating data from mixture of gamma with two components n = 1000 # number of observations weight = c( 0.6, 0.4 ) alpha = c( 12 , 1 ) beta = c( 3 , 2 ) data = rmixgamma( n = n, weight = weight, alpha = alpha, beta = beta ) # plot for simulation data hist( data, prob = TRUE, nclass = 50, col = "gray" ) x = seq( 0, 10, 0.05 ) truth = dmixgamma( x, weight, alpha, beta ) lines( x, truth, lwd = 2 ) # Runing bdmcmc algorithm for the above simulation data set bmixgamma.obj = bmixgamma( data, iter = 1000 ) summary( bmixgamma.obj ) plot( bmixgamma.obj ) ## End(Not run)