example

library(dsample)

Multi-Modes Detection

Please run demo(mix2) and demo(mix3).

6.1 Space Shuttle Challenger disaster

Data are taken from Dalal, Fowlkes, and Hoadley (1989). Details are described in Dezfuli et al. (2009) on pages 144–146.

expr <- str2expression("
  lp <- 0
  for(i in 1:len) lp <- lp + 
    y[i] * log(exp(alpha + beta*temp[i])/(1+exp(alpha + beta*temp[i])))
  for(i in 1:len) lp <- lp + 
    (1-y[i])*log(1/(1+exp(alpha + beta*temp[i])))
  lp <- lp + alpha - exp(alpha)/b
  lp <- exp(lp)
")

sets <- list(
  alpha=runif(n=nd, min=10, max=20), 
  beta=runif(n=nd, min=-0.3, max=-0.15)
)

smp <- dsample(expr=expr, rpmat=sets, nk=1e3, n=1e3)
op <- summary(smp)
op$means
#>      alpha       beta 
#> 15.0410844 -0.2324951
op$stdevs
#>      alpha       beta 
#> 1.31718727 0.02120767

6.2 Beetle

Data are taken from Prentice (1976). Details are described in OpenBUGS Examples Vol 2. Beetles.

expr <- str2expression("
  sigma <- exp(log.sigma)
  m1 <- exp(log.m1)
  
  lp <- 0
  for(i in 1:len) lp <- lp + 
    yi[i]*m1*log((exp((wi[i]-mu)/sigma)/(1+exp((wi[i]-mu)/sigma))))
  for(i in 1:len) lp <- lp + 
    (ni[i]-yi[i])*log(( 1- (exp((wi[i]-mu)/sigma)/(1+exp((wi[i]-mu)/sigma)))^m1 ))
  lp <- lp + (a-1)*log.m1 - 2*(e+1)*log.sigma
  lp <- lp - 0.5*((mu-c1)/d)^2
  lp <- lp - m1/b - 1/(f*sigma^2)
  lp <- exp(lp)
")

sets <- list(
  mu=runif(nd, min=1.75, max=1.85), 
  log.sigma=runif(nd, min=-5, max=-3), 
  log.m1=runif(nd, min=-2, max=0.1)
)

smp <- dsample(expr=expr, rpmat=sets, nk=1e3, n=1e3)
op <- summary(smp)
op$means
#>        mu log.sigma    log.m1 
#>  1.814511 -4.100180 -1.203697
op$stdevs
#>         mu  log.sigma     log.m1 
#> 0.01447934 0.25005453 0.34182293

6.3 Dugongs

Data are taken from Ratkowsky (1986). Details are described in OpenBUGS Examples Vol 2.Dugongs.

expr <- str2expression("
  lp <- (len/2 + k - 1)*log(tau)
  for(i in 1:len) lp <- lp - 
    tau*0.5*(y.length[i] - alpha+beta*gamma^x.age[i])^2
  lp <- lp - tau*k - tau.alpha*alpha^2*0.5 - tau.beta*beta^2*0.5
  lp <- exp(lp)
")

sets <- list(
  alpha=runif(nd, min=2, max=3), 
  beta=runif(nd, min=0.5, max=1.5), 
  gamma=runif(nd, min=0.5, max=1.5), 
  tau=runif(nd, min=0.2, max=200)
)

smp <- dsample(expr=expr, rpmat=sets, nk=1e3, n=1e3)
op <- summary(smp)
op$means
#>       alpha        beta       gamma         tau 
#>   2.6458769   0.9649206   0.8455442 114.6309473
op$stdevs
#>       alpha        beta       gamma         tau 
#>  0.17992098  0.15649336  0.09861807 40.79853512

6.4 British Coal Mining Data

Data are taken from Diggle and Marron (1988).

expr <- str2expression("
  ll <- 0
  ll <- ll + (cum.x.until.k[kappa]-0.5)*log(theta) + 
        (cum.x.after.k[kappa]-0.5)*log(lambda) - 
        kappa*theta -  (len-kappa)*lambda
  lp <- ll  + 1.5*log(alpha) + 1.5*log(beta) - 
        (theta+1)*alpha - (lambda+1)*beta
  lp <- exp(lp)
")

sets <- list(
  kappa=sample(x=30:50, size=nd, replace=TRUE),
  theta=runif(nd, min=2.2, max=4),
  lambda=runif(nd, min=0.6, max=1.4),
  alpha=runif(nd, min=0, max=2),
  beta=runif(nd, min=0, max=4)
)

smp <- dsample(expr=expr, rpmat=sets, nk=1e3, n=1e3)
op <- summary(smp)
op$means
#>      kappa      theta     lambda      alpha       beta 
#> 40.0470000  3.0769057  0.9175435  0.6127736  1.3342198
op$stdevs
#>     kappa     theta    lambda     alpha      beta 
#> 2.7659431 0.3010445 0.1250694 0.4005116 0.8059321

6.5 Nuclear Pump Data

Data are taken from Gaver and O’Muircheartaigh (1987). Details are described in OpenBUGS Examples Vol 2..

expr <- str2expression("
  ll <- 0
  for(i in 1:len){
    sum.cmd <- gsub(' ', '', paste('ll <- ll +(failure[', i,']+alpha-1)*log(lambda', i,')'))
    eval(parse(text=sum.cmd))
  }
  for(i in 1:len){
    sum.cmd <- gsub(' ', '', paste('ll <- ll - (time[', i,']+bb)*lambda', i))
    eval(parse(text=sum.cmd))
  }
  
  lp <- ll + (10*alpha+gg-1)*log(bb) - delta*bb
  lp <- exp(lp)
")

sets <- list(
  bb=runif(nd, 0, 4),
  lambda1=runif(nd, 0, 0.2),
  lambda2=runif(nd, 0, 0.4),
  lambda3=runif(nd, 0, 0.25),
  lambda4=runif(nd, 0, 0.25),
  lambda5=runif(nd, 0, 2),
  lambda6=runif(nd, 0, 1.5),
  lambda7=runif(nd, 0, 2),
  lambda8=runif(nd, 0, 2),
  lambda9=runif(nd, 0, 4),
  lambda10=runif(nd, 0, 3.5)
)

smp <- dsample(expr=expr, rpmat=sets, nk=5e4, n=3e3)
op <- summary(smp)
op$means
#>         bb    lambda1    lambda2    lambda3    lambda4    lambda5    lambda6 
#> 1.18648090 0.05849313 0.10509804 0.08872813 0.11502873 0.64308619 0.60354310 
#>    lambda7    lambda8    lambda9   lambda10 
#> 0.66471406 0.60911674 1.56126854 1.92201563
op$stdevs
#>         bb    lambda1    lambda2    lambda3    lambda4    lambda5    lambda6 
#> 0.69196508 0.03641480 0.08606147 0.03959884 0.03555977 0.34804315 0.20101222 
#>    lambda7    lambda8    lambda9   lambda10 
#> 0.44382590 0.47142139 0.75903501 0.50815639

References

Dalal, Siddhartha R., Edward B. Fowlkes, and Bruce Hoadley. 1989. “Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure.” Journal of the American Statistical Association 84 (408): 945–57. http://www.jstor.org/stable/2290069.
Dezfuli, Homayoon, Dana Kelly, Curtis L. Smith, Kurt G. Vedros, and William J. Galyean. 2009. “Bayesian Inference for NASA Probabilistic Risk and Reliability Analysis.” In.
Diggle, Peter, and J. S. Marron. 1988. “Equivalence of Smoothing Parameter Selectors in Density and Intensity Estimation.” Journal of the American Statistical Association 83 (403): 793–800. http://www.jstor.org/stable/2289308.
Gaver, Donald P., and I. G. O’Muircheartaigh. 1987. “Robust Empirical Bayes Analyses of Event Rates.” Technometrics 29 (1): 1–15. https://doi.org/10.1080/00401706.1987.10488178.
Liang, Faming, Chuanhai Liu, and Raymond J Carroll. 2007. “Stochastic Approximation in Monte Carlo Computation.” Journal of the American Statistical Association 102 (477): 305–20. https://doi.org/10.1198/016214506000001202.
Lunn, David J., Andrew Thomas, Nicky Best, and David Spiegelhalter. 2000. “WinBUGS a Bayesian Modelling Framework: Concepts, Structure, and Extensibility.” Statistics and Computing 10 (4): 325–37. https://doi.org/http://dx.doi.org/10.1023/A:1008929526011.
Prentice, Ross L. 1976. “A Generalization of the Probit and Logit Methods for Dose Response Curves.” Biometrics 32 (4): 761–68. http://www.jstor.org/stable/2529262.
Ratkowsky, D. A. 1986. “Statistical Properties of Alternative Parameterizations of the von Bertalanffy Growth Curve.” Canadian Journal of Fisheries and Aquatic Sciences 43 (4): 742–47. https://doi.org/10.1139/f86-091.
Wang, Liqun, and James Fu. 2007. “A practical sampling approach for a Bayesian mixture model with unknown number of components.” Statistical Papers 48 (4): 631–53. https://doi.org/10.1007/s00362-007-0361-4.
Wang, Liqun, and Chel Hee Lee. 2014. “Discretization-Based Direct Random Sample Generation.” Computational Statistics & Data Analysis 71: 1001–10. https://doi.org/https://doi.org/10.1016/j.csda.2013.06.011.
West, Mike. 1993. “Approximating Posterior Distributions by Mixture.” Journal of the Royal Statistical Society. Series B (Methodological) 55 (2): 409–22. http://www.jstor.org/stable/2346202.