Package 'LaplacesDemon'

Description

Provides a complete environment for Bayesian inference using a variety of different samplers (see ?LaplacesDemon for an overview).

Details

The DESCRIPTION file:

Index of help topics:

ABB                     Approximate Bayesian Bootstrap
AcceptanceRate          Acceptance Rate
BMK.Diagnostic          BMK Convergence Diagnostic
BayesFactor             Bayes Factor
BayesTheorem            Bayes' Theorem
BayesianBootstrap       The Bayesian Bootstrap
BigData                 Big Data
Blocks                  Blocks
CSF                     Cumulative Sample Function
CenterScale             Centering and Scaling
Combine                 Combine Demonoid Objects
Consort                 Consort with Laplace's Demon
Cov2Prec                Precision
ESS                     Effective Sample Size due to Autocorrelation
GIV                     Generate Initial Values
GaussHermiteQuadRule    Math Utility Functions
Gelfand.Diagnostic      Gelfand's Convergence Diagnostic
Gelman.Diagnostic       Gelman and Rubin's MCMC Convergence Diagnostic
Geweke.Diagnostic       Geweke's Convergence Diagnostic
Hangartner.Diagnostic   Hangartner's Convergence Diagnostic
Heidelberger.Diagnostic
                        Heidelberger and Welch's MCMC Convergence
                        Diagnostic
IAT                     Integrated Autocorrelation Time
Importance              Variable Importance
IterativeQuadrature     Iterative Quadrature
Juxtapose               Juxtapose MCMC Algorithm Inefficiency
KLD                     Kullback-Leibler Divergence (KLD)
KS.Diagnostic           Kolmogorov-Smirnov Convergence Diagnostic
LML                     Logarithm of the Marginal Likelihood
LPL.interval            Lowest Posterior Loss Interval
LaplaceApproximation    Laplace Approximation
LaplacesDemon           Laplace's Demon
LaplacesDemon-package   Complete Environment for Bayesian Inference
LaplacesDemon.RAM       LaplacesDemon RAM Estimate
Levene.Test             Levene's Test
LossMatrix              Loss Matrix
MCSE                    Monte Carlo Standard Error
MISS                    Multiple Imputation Sequential Sampling
MinnesotaPrior          Minnesota Prior
Mode                    The Mode(s) of a Vector
Model.Spec.Time         Model Specification Time
PMC                     Population Monte Carlo
PMC.RAM                 PMC RAM Estimate
PosteriorChecks         Posterior Checks
Raftery.Diagnostic      Raftery and Lewis's diagnostic
RejectionSampling       Rejection Sampling
SIR                     Sampling Importance Resampling
SensitivityAnalysis     Sensitivity Analysis
Stick                   Truncated Stick-Breaking
Thin                    Thin
Validate                Holdout Validation
VariationalBayes        Variational Bayes
WAIC                    Widely Applicable Information Criterion
as.covar                Proposal Covariance
as.indicator.matrix     Matrix Utility Functions
as.initial.values       Initial Values
as.parm.names           Parameter Names
as.ppc                  As Posterior Predictive Check
burnin                  Burn-in
caterpillar.plot        Caterpillar Plot
cloglog                 The log-log and complementary log-log functions
cond.plot               Conditional Plots
dStick                  Truncated Stick-Breaking Prior Distribution
dalaplace               Asymmetric Laplace Distribution: Univariate
dallaplace              Asymmetric Log-Laplace Distribution
daml                    Asymmetric Multivariate Laplace Distribution
dbern                   Bernoulli Distribution
dcat                    Categorical Distribution
dcrmrf                  Continuous Relaxation of a Markov Random Field
                        Distribution
ddirichlet              Dirichlet Distribution
de.Finetti.Game         de Finetti's Game
deburn                  De-Burn
delicit                 Prior Elicitation
demonchoice             Demon Choice Data Set
demonfx                 Demon FX Data Set
demonsessions           Demon Sessions Data Set
demonsnacks             Demon Snacks Data Set
demontexas              Demon Space-Time Data Set
dgpd                    Generalized Pareto Distribution
dgpois                  Generalized Poisson Distribution
dhalfcauchy             Half-Cauchy Distribution
dhalfnorm               Half-Normal Distribution
dhalft                  Half-t Distribution
dhs                     Horseshoe Distribution
dhuangwand              Huang-Wand Distribution
dhyperg                 Hyperprior-g Prior and Zellner's g-Prior
dinvbeta                Inverse Beta Distribution
dinvchisq               (Scaled) Inverse Chi-Squared Distribution
dinvgamma               Inverse Gamma Distribution
dinvgaussian            Inverse Gaussian Distribution
dinvmatrixgamma         Inverse Matrix Gamma Distribution
dinvwishart             Inverse Wishart Distribution
dinvwishartc            Inverse Wishart Distribution: Cholesky
                        Parameterization
dlaplace                Laplace Distribution: Univariate Symmetric
dlaplacem               Mixture of Laplace Distributions
dlaplacep               Laplace Distribution: Precision
                        Parameterization
dlasso                  LASSO Distribution
dllaplace               Log-Laplace Distribution: Univariate Symmetric
dlnormp                 Log-Normal Distribution: Precision
                        Parameterization
dmatrixgamma            Matrix Gamma Distribution
dmatrixnorm             Matrix Normal Distribution
dmvc                    Multivariate Cauchy Distribution
dmvcc                   Multivariate Cauchy Distribution: Cholesky
                        Parameterization
dmvcp                   Multivariate Cauchy Distribution: Precision
                        Parameterization
dmvcpc                  Multivariate Cauchy Distribution:
                        Precision-Cholesky Parameterization
dmvl                    Multivariate Laplace Distribution
dmvlc                   Multivariate Laplace Distribution: Cholesky
                        Parameterization
dmvn                    Multivariate Normal Distribution
dmvnc                   Multivariate Normal Distribution: Cholesky
                        Parameterization
dmvnp                   Multivariate Normal Distribution: Precision
                        Parameterization
dmvnpc                  Multivariate Normal Distribution:
                        Precision-Cholesky Parameterization
dmvpe                   Multivariate Power Exponential Distribution
dmvpec                  Multivariate Power Exponential Distribution:
                        Cholesky Parameterization
dmvpolya                Multivariate Polya Distribution
dmvt                    Multivariate t Distribution
dmvtc                   Multivariate t Distribution: Cholesky
                        Parameterization
dmvtp                   Multivariate t Distribution: Precision
                        Parameterization
dmvtpc                  Multivariate t Distribution: Precision-Cholesky
                        Parameterization
dnorminvwishart         Normal-Inverse-Wishart Distribution
dnormlaplace            Normal-Laplace Distribution: Univariate
                        Asymmetric
dnormm                  Mixture of Normal Distributions
dnormp                  Normal Distribution: Precision Parameterization
dnormv                  Normal Distribution: Variance Parameterization
dnormwishart            Normal-Wishart Distribution
dpareto                 Pareto Distribution
dpe                     Power Exponential Distribution: Univariate
                        Symmetric
dsdlaplace              Skew Discrete Laplace Distribution: Univariate
dsiw                    Scaled Inverse Wishart Distribution
dslaplace               Skew-Laplace Distribution: Univariate
dst                     Student t Distribution: Univariate
dstp                    Student t Distribution: Precision
                        Parameterization
dtrunc                  Truncated Distributions
dwishart                Wishart Distribution
dwishartc               Wishart Distribution: Cholesky Parameterization
dyangberger             Yang-Berger Distribution
interval                Constrain to Interval
is.appeased             Appeased
is.bayesfactor          Logical Check of Classes
is.bayesian             Logical Check of a Bayesian Model
is.constant             Logical Check of a Constant
is.constrained          Logical Check of Constraints
is.data                 Logical Check of Data
is.model                Logical Check of a Model
is.proper               Logical Check of Propriety
is.stationary           Logical Check of Stationarity
joint.density.plot      Joint Density Plot
joint.pr.plot           Joint Probability Region Plot
logit                   The logit and inverse-logit functions
p.interval              Probability Interval
plot.bmk                Plot Hellinger Distances
plot.demonoid           Plot samples from the output of Laplace's Demon
plot.demonoid.ppc       Plots of Posterior Predictive Checks
plot.importance         Plot Variable Importance
plot.iterquad           Plot the output of 'IterativeQuadrature'
plot.iterquad.ppc       Plots of Posterior Predictive Checks
plot.juxtapose          Plot MCMC Juxtaposition
plot.laplace            Plot the output of 'LaplaceApproximation'
plot.laplace.ppc        Plots of Posterior Predictive Checks
plot.miss               Plot samples from the output of MISS
plot.pmc                Plot samples from the output of PMC
plot.pmc.ppc            Plots of Posterior Predictive Checks
plot.vb                 Plot the output of 'VariationalBayes'
plot.vb.ppc             Plots of Posterior Predictive Checks
plotMatrix              Plot a Numerical Matrix
plotSamples             Plot Samples
predict.demonoid        Posterior Predictive Checks
predict.iterquad        Posterior Predictive Checks
predict.laplace         Posterior Predictive Checks
predict.pmc             Posterior Predictive Checks
predict.vb              Posterior Predictive Checks
print.demonoid          Print an object of class 'demonoid' to the
                        screen.
print.heidelberger      Print an object of class 'heidelberger' to the
                        screen.
print.iterquad          Print an object of class 'iterquad' to the
                        screen.
print.laplace           Print an object of class 'laplace' to the
                        screen.
print.miss              Print an object of class 'miss' to the screen.
print.pmc               Print an object of class 'pmc' to the screen.
print.raftery           Print an object of class 'raftery' to the
                        screen.
print.vb                Print an object of class 'vb' to the screen.
server_Listening        Server Listening
summary.demonoid.ppc    Posterior Predictive Check Summary
summary.iterquad.ppc    Posterior Predictive Check Summary
summary.laplace.ppc     Posterior Predictive Check Summary
summary.miss            MISS Summary
summary.pmc.ppc         Posterior Predictive Check Summary
summary.vb.ppc          Posterior Predictive Check Summary

The goal of LaplacesDemon, often referred to as LD, is to provide a complete and self-contained Bayesian environment within R. For example, this package includes dozens of MCMC algorithms, Laplace Approximation, iterative quadrature, variational Bayes, parallelization, big data, PMC, over 100 examples in the “Examples” vignette, dozens of additional probability distributions, numerous MCMC diagnostics, Bayes factors, posterior predictive checks, a variety of plots, elicitation, parameter and variable importance, Bayesian forms of test statistics (such as Durbin-Watson, Jarque-Bera, etc.), validation, and numerous additional utility functions, such as functions for multimodality, matrices, or timing your model specification. Other vignettes include an introduction to Bayesian inference, as well as a tutorial.

No further development of this package is currently being done as the original maintainer has stopped working on the package. Contributions to this package are welcome at https://github.com/LaplacesDemonR/LaplacesDemon.

The main function in this package is the LaplacesDemon function, and the best place to start is probably with the LaplacesDemon Tutorial vignette.

Author(s)

Byron Hall [aut], Martina Hall [aut], Statisticat, LLC [aut], Eric Brown [ctb], Richard Hermanson [ctb], Emmanuel Charpentier [ctb], Daniel Heck [ctb], Stephane Laurent [ctb], Quentin F. Gronau [ctb], Henrik Singmann [cre]

Maintainer: Henrik Singmann <[email protected]>

Description

This function performs multiple imputation (MI) with the Approximate Bayesian Bootstrap (ABB) of Rubin and Schenker (1986).

Usage

ABB(X, K=1)

Arguments

Details

The Approximate Bayesian Bootstrap (ABB) is a modified form of the BayesianBootstrap (Rubin, 1981) that is used for multiple imputation (MI). Imputation is a family of statistical methods for replacing missing values with estimates. Introduced by Rubin and Schenker (1986) and Rubin (1987), MI is a family of imputation methods that includes multiple estimates, and therefore includes variability of the estimates.

The data, $\textbf{X}$, are assumed to be independent and identically distributed (IID), contain both observed and missing values, and its missing values are assumed to be ignorable (meaning enough information is available in the data that the missingness mechanism can be ignored, if the information is used properly) and Missing Completely At Random (MCAR). When ABB is used in conjunction with a propensity score (described below), missing values may be Missing At Random (MAR).

ABB does not add auxiliary information, but performs imputation with two sampling (with replacement) steps. First, $\textbf{X}^\star_{obs}$ is sampled from $\textbf{X}_{obs}$. Then, $\textbf{X}^\star_{mis}$ is sampled from $\textbf{X}^\star_{obs}$. The result is a sample of the posterior predictive distribution of $(\textbf{X}_{mis}|\textbf{X}_{obs})$. The first sampling step is also known as hotdeck imputation, and the second sampling step changes the variance. Since auxiliary information is not included, ABB is appropriate for missing values that are ignorable and MCAR.

Auxiliary information may be included in the process of imputation by introducing a propensity score (Rosenbaum and Rubin, 1983; Rosenbaum and Rubin, 1984), which is an estimate of the probability of missingness. The propensity score is often the result of a binary logit model, where missingness is predicted as a function of other variables. The propensity scores are discretized into quantile-based groups, usually quintiles. Each quintile must have both observed and missing values. ABB is applied to each quintile. This is called within-class imputation. It is assumed that the missing mechanism depends only on the variables used to estimate the propensity score.

With $K=1$, ABB may be used in MCMC, such as in LaplacesDemon, more commonly along with a propensity score for missingness. MI is performed, despite $K=1$, because imputation occurs at each MCMC iteration. The practical advantage of this form of imputation is the ease with which it may be implemented. For example, full-likelihood imputation should perform better, but requires a chain to be updated for each missing value.

An example of a limitation of ABB with propensity scores is to consider imputing missing values of income from age in a context where age and income have a positive relationship, and where the highest incomes are missing systematically. ABB with propensity scores should impute these highest missing incomes given the highest observed ages, but is unable to infer beyond the observed data.

ABB has been extended (Parzen et al., 2005) to reduce bias, by introducing a correction factor that is applied to the MI variance estimate. This correction may be applied to output from ABB.

Value

This function returns a list with $K$ components, one for each set of imputations. Each component contains a vector of imputations equal in length to the number of missing values in the data.

ABB does not currently return the mean of the imputations, or the between-imputation variance or within-imputation variance.

Author(s)

Statisticat, LLC [email protected]

References

Parzen, M., Lipsitz, S.R., and Fitzmaurice, G.M. (2005). "A Note on Reducing the Bias of the Approximate Bayesian Bootstrap Imputation Variance Estimator". Biometrika, 92, 4, p. 971–974.

Rosenbaum, P.R. and Rubin, D.B. (1983). "The Central Role of the Propensity Score in Observational Studies for Causal Effects". Biometrika, 70, p. 41–55.

Rosenbaum, P.R. and Rubin, D.B. (1984). "Reducing Bias in Observational Studies Using Subclassification in the Propensity Score". Journal of the American Statistical Association, 79, p. 516–524.

Rubin, D.B. (1981). "The Bayesian Bootstrap". Annals of Statistics, 9, p. 130–134.

Rubin, D.B. (1987). "Multiple Imputation for Nonresponse in Surveys". John Wiley and Sons: New York, NY.

Rubin, D.B. and Schenker, N. (1986). "Multiple Imputation for Interval Estimation from Simple Random Samples with Ignorable Nonresponse". Journal of the American Statistical Association, 81, p. 366–374.

See Also

BayesianBootstrap, LaplacesDemon, and MISS.

Examples

library(LaplacesDemon)

### Create Data
J <- 10 #Number of variables
m <- 20 #Number of missings
N <- 50 #Number of records
mu <- runif(J, 0, 100)
sigma <- runif(J, 0, 100)
X <- matrix(0, N, J)
for (j in 1:J) X[,j] <- rnorm(N, mu[j], sigma[j])

### Create Missing Values
M1 <- rep(0, N*J)
M2 <- sample(N*J, m)
M1[M2] <- 1
M <- matrix(M1, N, J)
X <- ifelse(M == 1, NA, X)

### Approximate Bayesian Bootstrap
imp <- ABB(X, K=1)

### Replace Missing Values in X (when K=1)
X.imp <- X
X.imp[which(is.na(X.imp))] <- unlist(imp)
X.imp

Description

The Acceptance.Rate function calculates the acceptance rate per chain from a matrix of posterior MCMC samples.

Usage

AcceptanceRate(x)

Arguments

Details

The acceptance rate of an MCMC algorithm is the percentage of iterations in which the proposals were accepted.

Optimal Acceptance Rates
The optimal acceptance rate varies with the number of parameters and by algorithm. Algorithms with componentwise Gaussian proposals have an optimal acceptance rate of 0.44, regardless of the number of parameters. Algorithms that update with multivariate Gaussian proposals tend to have an optimal acceptance rate that ranges from 0.44 for one parameter (one IID Gaussian target distribution) to 0.234 for an infinite number of parameters (IID Gaussian target distributions), and 0.234 is approached quickly as the number of parameters increases. The AHMC, HMC, and THMC algorithms have an optimal acceptance rate of 0.67, except with the algorithm specification L=1, where the optimal acceptance rate is 0.574. The target acceptance rate is specified in HMCDA and NUTS, and the recommended rate is 0.65 and 0.60 respectively. Some algorithms have an acceptance rate of 1, such as AGG, ESS, GG, GS (MISS only), SGLD, or Slice.

Global and Local Acceptance Rates
LaplacesDemon reports the global acceptance rate for the un-thinned chains. However, componentwise algorithms make a proposal per parameter, and therefore have a local acceptance rate for each parameter. Since only the global acceptance rate is reported, the AcceptanceRate function may be used to calculate the local acceptance rates from a matrix of un-thinned posterior samples.

Thinning
Thinned samples tend to have higher local acceptance rates than un-thinned samples. With enough samples and enough thinning, local acceptance rates approach 1. Local acceptance rates do not need to approach the optimal acceptance rates above. Conversely, local acceptance rates do not need to approach 1, because too much information may possibly be discarded by thinning. For more information on thinning, see the Thin function.

Diagnostics
The AcceptanceRate function may be used to calculate local acceptance rates on a matrix of thinned or un-thinned samples. Any chain with a local acceptance rate that is an outlier may be studied for reasons that may cause the outlier. A local acceptance rate outlier does not violate theory and is often acceptable, but may indicate a potential problem. Only some of the many potential problems include: identifiability, model misspecification, multicollinearity, multimodality, choice of prior distributions, or becoming trapped in a low-probability space. The solution to local acceptance rate outliers tends to be either changing the MCMC algorithm or re-specifying the model or priors. For example, an MCMC algorithm that makes multivariate Gaussian proposals for a large number of parameters may have low global and local acceptance rates when far from the target distributions.

Value

The AcceptanceRate function returns a vector of acceptance rates, one for each chain.

Author(s)

Statisticat, LLC. [email protected]

See Also

LaplacesDemon, MISS, PosteriorChecks, and Thin.

Examples

library(LaplacesDemon)
AcceptanceRate(matrix(rnorm(5000),1000,5))

Description

This function returns the most recent covariance matrix or a list of blocking covariance matrices from an object of class demonoid, the most recent covariance matrix from iterquad, laplace, or vb, the most recent covariance matrix from the chain with the lowest deviance in an object of class demonoid.hpc, and a number of covariance matrices of an object of class pmc equal to the number of mixture components. The returned covariance matrix or matrices are intended to be the initial proposal covariance matrix or matrices for future updates. A variance vector from an object of class demonoid or demonoid.hpc is converted to a covariance matrix.

Usage

as.covar(x)

Arguments

Details

Unless it is known beforehand how many iterations are required for iterative quadrature, Laplace Approximation, or Variational Bayes to converge, MCMC to appear converged, or the normalized perplexity to stabilize in PMC, multiple updates are necessary. An additional update, however, should not begin with the same proposal covariance matrix or matrices as the original update, because it will have to repeat the work already accomplished. For this reason, the as.covar function may be used at the end of an update to change the previous initial values to the latest values.

The as.covar function is most helpful with objects of class pmc that have multiple mixture components. For more information, see PMC.

Value

The returned value is a matrix (or array in the case of PMC with multiple mixture components) of the latest observed or proposal covariance, which may now be used as an initial proposal covariance matrix or matrices for a future update.

Author(s)

Statisticat, LLC [email protected]

See Also

IterativeQuadrature, LaplaceApproximation, LaplacesDemon, LaplacesDemon.hpc, PMC, and VariationalBayes.

Description

This function returns the most recent posterior samples from an object of class demonoid or demonoid.hpc, the posterior means of an object of class iterquad, the posterior modes of an object of class laplace or vb, the posterior means of an object of class pmc with one mixture component, or the latest means of the importance sampling distribution of an object of class pmc with multiple mixture components. The returned values are intended to be the initial values for future updates.

Usage

as.initial.values(x)

Arguments

Details

Unless it is known beforehand how many iterations are required for IterativeQuadrature, LaplaceApproximation, or VariationalBayes to converge, MCMC in LaplacesDemon to appear converged, or the normalized perplexity to stabilize in PMC, multiple updates are necessary. An additional update, however, should not begin with the same initial values as the original update, because it will have to repeat the work already accomplished. For this reason, the as.initial.values function may be used at the end of an update to change the previous initial values to the latest values.

When using LaplacesDemon.hpc, as.initial.values should be used when the output is of class demonoid.hpc, before the Combine function is used to combine the multiple chains for use with Consort and other functions, because the Combine function returns an object of class demonoid, and the number of chains will become unknown. The Consort function may suggest using as.initial.values, but when applied to an object of class demonoid, it will return the latest values as if there were only one chain.

Value

The returned value is a vector (or matrix in the case of an object of class demonoid.hpc, or pmc with multiple mixture components) of the latest values, which may now be used as initial values for a future update.

Author(s)

Statisticat, LLC. [email protected]

See Also

Combine, IterativeQuadrature, LaplaceApproximation, LaplacesDemon, LaplacesDemon.hpc, PMC, and VariationalBayes.

Description

This function creates a vector of parameter names from a list of parameters, and the list may contain any combination of scalars, vectors, matrices, upper-triangular matrices, and arrays.

Usage

as.parm.names(x, uppertri=NULL)

Arguments

Details

Each model function for IterativeQuadrature, LaplaceApproximation, LaplacesDemon, PMC, or VariationalBayes requires a vector of parameters (specified at first as Initial.Values) and a list of data. One component in the list of data must be named parm.names. Each element of parm.names is a name associated with the corresponding parameter in Initial.Values.

The parm.names vector is easy to program explicitly for a simple model, but can require considerably more programming effort for more complicated models. The as.parm.names function is a utility function designed to minimize programming by the user.

For example, a simple model may only require parm.names <- c("alpha", "beta[1]", "beta[2]", "sigma"). A more complicated model may contain hundreds of parameters that are a combination of scalars, vectors, matrices, upper-triangular matrices, and arrays, and is the reason for the as.parm.names function. The code for the above is as.parm.names(list(alpha=0, beta=rep(0,2), sigma=0)).

In the case of an upper-triangular matrix, simply pass the full matrix to as.parm.names and indicate that only the upper-triangular will be used via the uppertri argument. For example, as.parm.names(list(beta=rep(0,J),U=diag(K)), uppertri=c(0,1)) creates parameter names for a vector of $\beta$ parameters of length $J$ and an upper-triangular matrix $\textbf{U}$ of dimension $K$.

Numerous examples may be found in the accompanying “Examples” vignette.

Value

This function returns a vector of parameter names.

Author(s)

Statisticat, LLC. [email protected]

See Also

IterativeQuadrature LaplaceApproximation, LaplacesDemon, PMC, and VariationalBayes.

Examples

library(LaplacesDemon)
N <- 100
J <- 5
y <- rnorm(N,0,1)
X <- matrix(runif(N*J,-2,2),N,J)
S <- diag(J)
T <- diag(2)
mon.names <- c("LP","sigma")
parm.names <- as.parm.names(list(log.sigma=0, beta=rep(0,J), S=diag(J),
     T=diag(2)), uppertri=c(0,0,0,1))
MyData <- list(J=J, N=N, S=S, T=T, X=X, mon.names=mon.names,
     parm.names=parm.names, y=y)
MyData

Description

This function converts an object of class demonoid.val to an object of class demonoid.ppc.

Usage

as.ppc(x, set=3)

Arguments

Details

After using the Validate function for holdout validation, it is often suggested to perform posterior predictive checks. The as.ppc function converts the output object of Validate, which is an object of class demonoid.val, to an object of class demonoid.ppc. The returned object is the same as if it were created with the predict.demonoid function, rather than the Validate function.

After this conversion, the user may use posterior predictive checks, as usual, with the summary.demonoid.ppc function.

Value

The returned object is an object of class demonoid.ppc.

Author(s)

Statisticat, LLC. [email protected]

See Also

predict.demonoid, summary.demonoid.ppc, and Validate.

Description

This function calculates Bayes factors for two or more fitted objects of class demonoid, iterquad, laplace, pmc, or vb that were estimated respectively with the LaplacesDemon, IterativeQuadrature, LaplaceApproximation, PMC, or VariationalBayes functions, and indicates the strength of evidence in favor of the hypothesis (that each model, $\mathcal{M}_i$, is better than another model, $\mathcal{M}_j$).

Usage

BayesFactor(x)

Arguments

Details

Introduced by Harold Jeffreys, a 'Bayes factor' is a Bayesian alternative to frequentist hypothesis testing that is most often used for the comparison of multiple models by hypothesis testing, usually to determine which model better fits the data (Jeffreys, 1961). Bayes factors are notoriously difficult to compute, and the Bayes factor is only defined when the marginal density of $\textbf{y}$ under each model is proper (see is.proper). However, the Bayes factor is easy to approximate with the Laplace-Metropolis estimator (Lewis and Raftery, 1997) and other methods of approximating the logarithm of the marginal likelihood (for more information, see LML).

Hypothesis testing with Bayes factors is more robust than frequentist hypothesis testing, since the Bayesian form avoids model selection bias, evaluates evidence in favor of the null hypothesis, includes model uncertainty, and allows non-nested models to be compared (though of course the model must have the same dependent variable). Also, frequentist significance tests become biased in favor of rejecting the null hypothesis with sufficiently large sample size.

The Bayes factor for comparing two models may be approximated as the ratio of the marginal likelihood of the data in model 1 and model 2. Formally, the Bayes factor in this case is

$B = \frac{p(\textbf{y}|\mathcal{M}_1)}{p(\textbf{y}|\mathcal{M}_2)} = \frac{\int p(\textbf{y}|\Theta_1,\mathcal{M}_1)p(\Theta_1|\mathcal{M}_1)d\Theta_1}{\int p(\textbf{y}|\Theta_2,\mathcal{M}_2)p(\Theta_2|\mathcal{M}_2)d\Theta_2}$

where $p(\textbf{y}|\mathcal{M}_1)$ is the marginal likelihood of the data in model 1.

The IterativeQuadrature, LaplaceApproximation, LaplacesDemon, PMC, and VariationalBayes functions each return the LML, the approximate logarithm of the marginal likelihood of the data, in each fitted object of class iterquad, laplace, demonoid, pmc, or vb. The BayesFactor function calculates matrix B, a matrix of Bayes factors, where each element of matrix B is a comparison of two models. Each Bayes factor is calculated as the exponentiated difference of LML of model 1 ($\mathcal{M}_1$) and LML of model 2 ($\mathcal{M}_2$), and the hypothesis for each element of matrix B is that the model associated with the row is greater than the model associated with the column. For example, element B[3,2] is the Bayes factor that model 3 is greater than model 2. The 'Strength of Evidence' aids in the interpretation (Jeffreys, 1961).

A table for the interpretation of the strength of evidence for Bayes factors is available at https://web.archive.org/web/20150214194051/http://www.bayesian-inference.com/bayesfactors.

Each Bayes factor, B, is the posterior odds in favor of the hypothesis divided by the prior odds in favor of the hypothesis, where the hypothesis is usually $\mathcal{M}_1 > \mathcal{M}_2$. For example, when B[3,2]=2, the data favor $\mathcal{M}_3$ over $\mathcal{M}_2$ with 2:1 odds.

It is also popular to consider the natural logarithm of the Bayes factor. The scale of the logged Bayes factor is the same above and below one, which is more appropriate for visual comparisons. For example, when comparing two Bayes factors at 0.5 and 2, the logarithm of these Bayes factors is -0.69 and 0.69.

Gelman finds Bayes factors generally to be irrelevant, because they compute the relative probabilities of the models conditional on one of them being true. Gelman prefers approaches that measure the distance of the data to each of the approximate models (Gelman et al., 2004, p. 180), such as with posterior predictive checks (see the predict.iterquad function regarding iterative quadrature, predict.laplace function in the context of Laplace Approximation, predict.demonoid function in the context of MCMC, predict.pmc function in the context of PMC, or predict.vb function in the context of Variational Bayes). Kass et al. (1995) asserts this can be done without assuming one model is the true model.

Value

BayesFactor returns an object of class bayesfactor that is a list with the following components:

Author(s)

Statisticat, LLC.

References

Gelman, A., Carlin, J., Stern, H., and Rubin, D. (2004). "Bayesian Data Analysis, Texts in Statistical Science, 2nd ed.". Chapman and Hall, London.

Jeffreys, H. (1961). "Theory of Probability, Third Edition". Oxford University Press: Oxford, England.

Kass, R.E. and Raftery, A.E. (1995). "Bayes Factors". Journal of the American Statistical Association, 90(430), p. 773–795.

Lewis, S.M. and Raftery, A.E. (1997). "Estimating Bayes Factors via Posterior Simulation with the Laplace-Metropolis Estimator". Journal of the American Statistical Association, 92, p. 648–655.

See Also

is.bayesfactor, is.proper, IterativeQuadrature, LaplaceApproximation, LaplacesDemon, LML, PMC, predict.demonoid, predict.iterquad, predict.laplace, predict.pmc, predict.vb, and VariationalBayes.

Examples

# The following example fits a model as Fit1, then adds a predictor, and
# fits another model, Fit2. The two models are compared with Bayes
# factors.

library(LaplacesDemon)

##############################  Demon Data  ###############################
data(demonsnacks)
J <- 2
y <- log(demonsnacks$Calories)
X <- cbind(1, as.matrix(log(demonsnacks[,10]+1)))
X[,2] <- CenterScale(X[,2])

#########################  Data List Preparation  #########################
mon.names <- "LP"
parm.names <- as.parm.names(list(beta=rep(0,J), sigma=0))
pos.beta <- grep("beta", parm.names)
pos.sigma <- grep("sigma", parm.names)
PGF <- function(Data) {
     beta <- rnorm(Data$J)
     sigma <- runif(1)
     return(c(beta, sigma))
     }
MyData <- list(J=J, PGF=PGF, X=X, mon.names=mon.names,
     parm.names=parm.names, pos.beta=pos.beta, pos.sigma=pos.sigma, y=y)

##########################  Model Specification  ##########################
Model <- function(parm, Data)
     {
     ### Parameters
     beta <- parm[Data$pos.beta]
     sigma <- interval(parm[Data$pos.sigma], 1e-100, Inf)
     parm[Data$pos.sigma] <- sigma
     ### Log-Prior
     beta.prior <- sum(dnormv(beta, 0, 1000, log=TRUE))
     sigma.prior <- dhalfcauchy(sigma, 25, log=TRUE)
     ### Log-Likelihood
     mu <- tcrossprod(Data$X, t(beta))
     LL <- sum(dnorm(Data$y, mu, sigma, log=TRUE))
     ### Log-Posterior
     LP <- LL + beta.prior + sigma.prior
     Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LP,
          yhat=rnorm(length(mu), mu, sigma), parm=parm)
     return(Modelout)
     }

############################  Initial Values  #############################
Initial.Values <- GIV(Model, MyData, PGF=TRUE)

########################  Laplace Approximation  ##########################
Fit1 <- LaplaceApproximation(Model, Initial.Values, Data=MyData,
     Iterations=10000)
Fit1

##############################  Demon Data  ###############################
data(demonsnacks)
J <- 3
y <- log(demonsnacks$Calories)
X <- cbind(1, as.matrix(demonsnacks[,c(7,8)]))
X[,2] <- CenterScale(X[,2])
X[,3] <- CenterScale(X[,3])

#########################  Data List Preparation  #########################
mon.names <- c("sigma","mu[1]")
parm.names <- as.parm.names(list(beta=rep(0,J), sigma=0))
pos.beta <- grep("beta", parm.names)
pos.sigma <- grep("sigma", parm.names)
PGF <- function(Data) return(c(rnormv(Data$J,0,10), rhalfcauchy(1,5)))
MyData <- list(J=J, PGF=PGF, X=X, mon.names=mon.names,
     parm.names=parm.names, pos.beta=pos.beta, pos.sigma=pos.sigma, y=y)

############################  Initial Values  #############################
Initial.Values <- GIV(Model, MyData, PGF=TRUE)

########################  Laplace Approximation  ##########################
Fit2 <- LaplaceApproximation(Model, Initial.Values, Data=MyData,
     Iterations=10000)
Fit2

#############################  Bayes Factor  ##############################
Model.list <- list(M1=Fit1, M2=Fit2)
BayesFactor(Model.list)

Description

This function performs the Bayesian bootstrap of Rubin (1981), returning either bootstrapped weights or statistics.

Usage

BayesianBootstrap(X, n=1000, Method="weights", Status=NULL)

Arguments

Details

The term, ‘bootstrap’, comes from the German novel Adventures of Baron Munchausen by Rudolph Raspe, in which the hero saves himself from drowning by pulling on his own bootstraps. The idea of the statistical bootstrap is to evaluate properties of an estimator through the empirical, rather than theoretical, CDF.

Rubin (1981) introduced the Bayesian bootstrap. In contrast to the frequentist bootstrap which simulates the sampling distribution of a statistic estimating a parameter, the Bayesian bootstrap simulates the posterior distribution.

The data, $\textbf{X}$, are assumed to be independent and identically distributed (IID), and to be a representative sample of the larger (bootstrapped) population. Given that the data has $N$ rows in one bootstrap replication, the row weights are sampled from a Dirichlet distribution with all $N$ concentration parameters equal to $1$ (a uniform distribution over an open standard $N-1$ simplex). The distributions of a parameter inferred from considering many samples of weights are interpretable as posterior distributions on that parameter.

The Bayesian bootstrap is useful for estimating marginal posterior covariance and standard deviations for the posterior modes of LaplaceApproximation, especially when the model dimension (the number of parameters) is large enough that estimating the Hessian matrix of second partial derivatives is too computationally demanding.

Just as with the frequentist bootstrap, inappropriate use of the Bayesian bootstrap can lead to inappropriate inferences. The Bayesian bootstrap violates the likelihood principle, because the evaluation of a statistic of interest depends on data sets other than the observed data set. For more information on the likelihood principle, see https://web.archive.org/web/20150213002158/http://www.bayesian-inference.com/likelihood#likelihoodprinciple.

The BayesianBootstrap function has many uses, including creating test statistics on the population data given the observed data (supported here), imputation (with this variation: ABB), validation, and more.

Value

When Method="weights", this function returns a $N \times n$ matrix of weights, where the number of rows $N$ is equal to the number of rows in X.

For statistics, a matrix or array is returned, depending on the number of dimensions. The replicates are indexed by row in a matrix or in the first dimension of the array.

Author(s)

Bogumil Kaminski, [email protected] and Statisticat, LLC.

References

Rubin, D.B. (1981). "The Bayesian Bootstrap". The Annals of Statistics, 9(1), p. 130–134.

See Also

ABB, Hessian, LaplaceApproximation, and LaplacesDemon.

Examples

library(LaplacesDemon)

#Example 1: Samples
x <- 1:2
BB <- BayesianBootstrap(X=x, n=100, Method="weights"); BB

#Example 2: Mean, Univariate
x <- 1:2
BB <- BayesianBootstrap(X=x, n=100, Method=weighted.mean); BB

#Example 3: Mean, Multivariate
data(demonsnacks)
BB <- BayesianBootstrap(X=demonsnacks, n=100,
     Method=function(x,w) apply(x, 2, weighted.mean, w=w)); BB

#Example 4: Correlation
dye <- c(1.15, 1.70, 1.42, 1.38, 2.80, 4.70, 4.80, 1.41, 3.90)
efp <- c(1.38, 1.72, 1.59, 1.47, 1.66, 3.45, 3.87, 1.31, 3.75)
X <- matrix(c(dye,efp), length(dye), 2)
colnames(X) <- c("dye","efp")
BB <- BayesianBootstrap(X=X, n=100,
     Method=function(x,w) cov.wt(x, w, cor=TRUE)$cor); BB

#Example 5: Marginal Posterior Covariance
#The following example is commented out due to package build time.
#To run the following example, use the code from the examples in
#the LaplaceApproximation function for the data, model specification
#function, and initial values. Then perform the Laplace
#Approximation as below (with CovEst="Identity" and sir=FALSE) until
#convergence, set the latest initial values, then use the Bayesian
#bootstrap on the data, run the Laplace Approximation again to
#convergence, save the posterior modes, and repeat until S samples
#of the posterior modes are collected. Finally, calculate the
#parameter covariance or standard deviation.

#Fit <- LaplaceApproximation(Model, Initial.Values, Data=MyData,
#     Iterations=1000, Method="SPG",  CovEst="Identity", sir=FALSE)
#Initial.Values <- as.initial.values(Fit)
#S <- 100 #Number of bootstrapped sets of posterior modes (parameters)
#Z <- rbind(Fit$Summary1[,1]) #Bootstrapped parameters collected here
#N <- nrow(MyData$X) #Number of records
#MyData.B <- MyData
#for (s in 1:S) {
#     cat("\nIter:", s, "\n")
#     BB <- BayesianBootstrap(MyData$y, n=N)
#     z <- apply(BB, 2, function(x) sample.int(N, size=1, prob=x))
#     MyData.B$y <- MyData$y[z]
#     MyData.B$X <- MyData$X[z,]
#     Fit <- LaplaceApproximation(Model, Initial.Values, Data=MyData.B,
#          Iterations=1000, Method="SPG", CovEst="Identity", sir=FALSE)
#     Z <- rbind(Z, Fit$Summary1[,1])}
#cov(Z) #Bootstrapped marginal posterior covariance
#sqrt(diag(cov(Z))) #Bootstrapped marginal posterior standard deviations

Description

Bayes' theorem shows the relation between two conditional probabilities that are the reverse of each other. This theorem is named after Reverend Thomas Bayes (1702-1761), and is also referred to as Bayes' law or Bayes' rule (Bayes and Price, 1763). Bayes' theorem expresses the conditional probability, or ‘posterior probability’, of an event $A$ after $B$ is observed in terms of the 'prior probability' of $A$, prior probability of $B$, and the conditional probability of $B$ given $A$. Bayes' theorem is valid in all common interpretations of probability. This function provides one of several forms of calculations that are possible with Bayes' theorem.

Usage

BayesTheorem(PrA, PrBA)

Arguments

Details

Bayes' theorem provides an expression for the conditional probability of $A$ given $B$, which is equal to

$\Pr(A | B) = \frac{\Pr(B | A)\Pr(A)}{\Pr(B)}$

For example, suppose one asks the question: what is the probability of going to Hell, conditional on consorting (or given that a person consorts) with Laplace's Demon. By replacing $A$ with $Hell$ and $B$ with $Consort$, the question becomes

$\Pr(\mathrm{Hell} | \mathrm{Consort}) = \frac{\Pr(\mathrm{Consort} | \mathrm{Hell})\Pr(\mathrm{Hell})}{\Pr(\mathrm{Consort})}$

Note that a common fallacy is to assume that $\Pr(A | B) = \Pr(B | A)$, which is called the conditional probability fallacy.

Another way to state Bayes' theorem (and this is the form in the provided function) is

$\Pr(A_i | B) = \frac{\Pr(B | A_i)\Pr(A_i)}{\Pr(B | A_i)\Pr(A_i) +\dots+ \Pr(B | A_n)\Pr(A_n)}$

Let's examine our burning question, by replacing $A_i$ with Hell or Heaven, and replacing $B$ with Consort

• $\Pr(A_1) = \Pr(\mathrm{Hell})$

• $\Pr(A_2) = \Pr(\mathrm{Heaven})$

• $\Pr(B) = \Pr(\mathrm{Consort})$

• $\Pr(A_1 | B) = \Pr(\mathrm{Hell} | \mathrm{Consort})$

• $\Pr(A_2 | B) = \Pr(\mathrm{Heaven} | \mathrm{Consort})$

• $\Pr(B | A_1) = \Pr(\mathrm{Consort} | \mathrm{Hell})$

• $\Pr(B | A_2) = \Pr(\mathrm{Consort} | \mathrm{Heaven})$

Laplace's Demon was conjured and asked for some data. He was glad to oblige.

• 6 people consorted out of 9 who went to Hell.

• 5 people consorted out of 7 who went to Heaven.

• 75% of the population goes to Hell.

• 25% of the population goes to Heaven.

Now, Bayes' theorem is applied to the data. Four pieces are worked out as follows

• $\Pr(\mathrm{Consort} | \mathrm{Hell}) = 6/9 = 0.666$

• $\Pr(\mathrm{Consort} | \mathrm{Heaven}) = 5/7 = 0.714$

• $\Pr(\mathrm{Hell}) = 0.75$

• $\Pr(\mathrm{Heaven}) = 0.25$

Finally, the desired conditional probability $\Pr(\mathrm{Hell} | \mathrm{Consort})$ is calculated using Bayes' theorem

• $\Pr(\mathrm{Hell} | \mathrm{Consort}) = \frac{0.666(0.75)}{0.666(0.75) + 0.714(0.25)}$

• $\Pr(\mathrm{Hell} | \mathrm{Consort}) = 0.737$

The probability of someone consorting with Laplace's Demon and going to Hell is 73.7%, which is less than the prevalence of 75% in the population. According to these findings, consorting with Laplace's Demon does not increase the probability of going to Hell.

For an introduction to model-based Bayesian inference, see the accompanying vignette entitled “Bayesian Inference” or https://web.archive.org/web/20150206004608/http://www.bayesian-inference.com/bayesian.

Value

The BayesTheorem function returns the conditional probability of $A$ given $B$, known in Bayesian inference as the posterior. The returned object is of class bayestheorem.

Author(s)

Statisticat, LLC.

References

Bayes, T. and Price, R. (1763). "An Essay Towards Solving a Problem in the Doctrine of Chances". By the late Rev. Mr. Bayes, communicated by Mr. Price, in a letter to John Canton, M.A. and F.R.S. Philosophical Transactions of the Royal Statistical Society of London, 53, p. 370–418.

See Also

IterativeQuadrature, LaplaceApproximation, LaplacesDemon, PMC, and VariationalBayes.

Examples

# Pr(Hell|Consort) =
PrA <- c(0.75,0.25)
PrBA <- c(6/9, 5/7)
BayesTheorem(PrA, PrBA)

Description

This function enables Bayesian inference with data that is too large for computer memory (RAM) with the simplest method: reading in batches of data (where each batch is a section of rows), applying a function to the batch, and combining the results.

Usage

BigData(file, nrow, ncol, size=1, Method="add", CPUs=1, Type="PSOCK",
FUN, ...)

Arguments

Details

Big data is defined loosely here as data that is too large for computer memory (RAM). The BigData function uses the split-apply-combine strategy with a big data set. The unmanageable big data set is split into smaller, manageable pieces (batches), a function is applied to each batch, and results are combined.

Each iteration, the BigData function opens a connection to a big data set and keeps the connection open while the scan function reads in each batch of data (elsewhere, batches are often referred to chunks). A user-specified function is applied to each batch of data, the results are combined together, the connection is closed, and the results are returned.

As an introductory example, suppose a statistician updates a linear regression model, but the design matrix $\textbf{X}$ is too large for computer memory. Suppose the design matrix has 100 million rows, and the statistician specifies size=1e6. The statistician combines dependent variable $\textbf{y}$ with design matrix $\textbf{X}$. Each iteration in IterativeQuadrature, LaplaceApproximation, LaplacesDemon, PMC, or VariationalBayes, the BigData function sequentially reads in one million rows of the combined data $\textbf{X}$, calculates expectation vector $\mu$, and finally returns the sum of the log-likelihood. The sum of the log-likelihood is added together for all batches, and returned.

There are many limitations with this function.

This function is not fast, in the sense that the entire big data set is processed in batches, each iteration. With iterative methods, this may perform well, albeit slowly.

There are many functions that cannot be performed on batches, though most models in the Examples vignette may easily be updated with big data.

Large matrices of samples are unaddressed, only the data.

Although many (but not all) models may be estimated, many additional functions in this package will not work when applied after the model has updated. Instead, a batch or random sample of data (see the read.matrix function for sampling from big data) should be used in the usual way, in the Data argument, and the Model function coded in the usual way without the BigData function.

Parallel processing may be performed when the user specifies CPUs to be greater than one, implying that the specified number of CPUs exists and is available. Parallelization may be performed on a multicore computer or a computer cluster. Either a Simple Network of Workstations (SNOW) or Message Passing Interface (MPI) is used. Each call to BigData establishes and closes the parallelization, which is costly, and unfortunately results in copious output to the console. With small data sets, parallel processing may be slower, due to computer network communication. With larger data sets, the user should experience a faster run-time.

There have been several alternative approaches suggested for big data.

Huang and Gelman (2005) propose that the user creates batches by sampling from big data, updating a separate Bayesian model on each batch, and combining the results into a consensus posterior. This many-mini-model approach may be faster when feasible, because multiple models may be updated in parallel, say one per CPU. Such results will work with all functions in this package. With the many-mini-model approach, several methods are proposed for combining posterior samples from batch-level models, such as by using a normal approximation, updating from prior to posterior sequentially (the posterior from the last batch becomes the prior of the next batch), sample from the full posterior via importance sampling from the batched posteriors, and more.

Scott et al. (2013) propose a method that they call Consensus Monte Carlo, which consists of breaking the data down into chunks, calling each chunk a shard, and use a many-mini-model approach as well, but propose their own method of weighting the posteriors back together.

Balakrishnan and Madigan (2006) introduced a Sequential Monte Carlo (SMC) sampler, a refinement of an earlier proposal, that was designed for big data. It makes one pass through the massive data set, after an initial MCMC estimation on a small sample. Each particle is updated for each record, resulting in numerous evaluations per record.

Welling and Teh (2011) proposed a new class of MCMC sampler in which only a random sample of big data is used each iteration. The stochastic gradient Langevin dynamics (SGLD) algorithm is available in the LaplacesDemon function.

An important alternative to consider is using the ff package, where "ff" stands for fast access file. The ff package has been tested successfully with updating a model in LaplacesDemon. Once the big data set, say $\textbf{X}$, is an object of class ff_matrix, simply include it in the list of data as usual, and modify the Model specification function appropriately. For example, change mu <- tcrossprod(X, t(beta)) to mu <- tcrossprod(X[], t(beta)). The ff package is not included as a dependency in the LaplacesDemon package, so it must be installed and activated.

Value

The BigData function returns output that is the result of performing a user-specified function on batches of big data. Output is a matrix, and may have one or more column vectors.

Author(s)

Statisticat, LLC [email protected]

References

Balakrishnan, S. and Madigan, D. (2006). "A One-Pass Sequential Monte Carlo Method for Bayesian Analysis of Massive Datasets". Bayesian Analysis, 1(2), p. 345–362.

Huang, Z. and Gelman, A. (2005) "Sampling for Bayesian Computation with Large Datasets". SSRN eLibrary.

Scott, S.L., Blocker, A.W. and Bonassi, F.V. (2013). "Bayes and Big Data: The Consensus Monte Carlo Algorithm". In Bayes 250.

Welling, M. and Teh, Y.W. (2011). "Bayesian Learning via Stochastic Gradient Langevin Dynamics". Proceedings of the 28th International Conference on Machine Learning (ICML), p. 681–688.

See Also

IterativeQuadrature, LaplaceApproximation, LaplacesDemon, LaplacesDemon.RAM, PMC, PMC.RAM, read.matrix, and VariationalBayes.

Examples

### Below is an example of a linear regression model specification
### function in which BigData reads in a batch of 1,000 records of
### Data$N records from a data set that is too large to fully open
### in memory. The example simulates on 10,000 records, which is
### not big data; it's just a toy example. The data set is file X.csv,
### and the first column of matrix X is the dependent variable y. The
### user supplies a function to BigData along with parameters beta and
### sigma. When each batch of 1,000 records is read in,
### mu = XB is calculated, and then the LL is calculated as
### y ~ N(mu, sigma^2). These results are added together from all
### batches, and returned as LL.

library(LaplacesDemon)
N <- 10000
J <- 10 #Number of predictors, including the intercept
X <- matrix(1,N,J)
for (j in 2:J) {X[,j] <- rnorm(N,runif(1,-3,3),runif(1,0.1,1))}
beta.orig <- runif(J,-3,3)
e <- rnorm(N,0,0.1)
y <- as.vector(tcrossprod(beta.orig, X) + e)
mon.names <- c("LP","sigma")
parm.names <- as.parm.names(list(beta=rep(0,J), log.sigma=0))
PGF <- function(Data) return(c(rnormv(Data$J,0,0.01),
     log(rhalfcauchy(1,1))))
MyData <- list(J=J, PGF=PGF, N=N, mon.names=mon.names,
     parm.names=parm.names) #Notice that X and y are not included here
filename <- tempfile("X.csv")  
write.table(cbind(y,X), filename, sep=",", row.names=FALSE,
  col.names=FALSE)

Model <- function(parm, Data)
     {
     ### Parameters
     beta <- parm[1:Data$J]
     sigma <- exp(parm[Data$J+1])
     ### Log(Prior Densities)
     beta.prior <- sum(dnormv(beta, 0, 1000, log=TRUE))
     sigma.prior <- dhalfcauchy(sigma, 25, log=TRUE)
     ### Log-Likelihood
     LL <- BigData(file=filename, nrow=Data$N, ncol=Data$J+1, size=1000,
          Method="add", CPUs=1, Type="PSOCK",
          FUN=function(x, beta, sigma) sum(dnorm(x[,1], tcrossprod(x[,-1],
               t(beta)), sigma, log=TRUE)), beta, sigma)
     ### Log-Posterior
     LP <- LL + beta.prior + sigma.prior
     Modelout <- list(LP=LP, Dev=-2*LL, Monitor=c(LP,sigma),
          yhat=0,#rnorm(length(mu), mu, sigma),
          parm=parm)
     return(Modelout)
     }

### From here, the user may update the model as usual.

Description

The Blocks function returns a list of $N$ blocks of parameters, for use with some MCMC algorithms in the LaplacesDemon function. Blocks may be created either sequentially, or from a hierarchical clustering of the posterior correlation matrix.

Usage

Blocks(Initial.Values, N, PostCor=NULL)

Arguments

Details

Usually, there is more than one target distribution in MCMC, in which case it must be determined whether it is best to sample from target distributions individually, in groups, or all at once. Blockwise sampling (also called block updating) refers to splitting a multivariate vector into groups called blocks, and each block is sampled separately. A block may contain one or more parameters.

Parameters are usually grouped into blocks such that parameters within a block are as correlated as possible, and parameters between blocks are as independent as possible. This strategy retains as much of the parameter correlation as possible for blockwise sampling, as opposed to componentwise sampling where parameter correlation is ignored. The PosteriorChecks function can be used on the output of previous runs to find highly correlated parameters. See examples below.

Advantages of blockwise sampling are that a different MCMC algorithm may be used for each block (or parameter, for that matter), creating a more specialized approach (though different algorithms by block are not supported here), the acceptance of a newly proposed state is likely to be higher than sampling from all target distributions at once in high dimensions, and large proposal covariance matrices can be reduced in size, which is most helpful again in high dimensions.

Disadvantages of blockwise sampling are that correlations probably exist between parameters between blocks, and each block is updated while holding the other blocks constant, ignoring these correlations of parameters between blocks. Without simultaneously taking everything into account, the algorithm may converge slowly or never arrive at the proper solution. However, there are instances when it may be best when everything is not taken into account at once, such as in state-space models. Also, as the number of blocks increases, more computation is required, which slows the algorithm. In general, blockwise sampling allows a more specialized approach at the expense of accuracy, generalization, and speed. Blockwise sampling is offered in the following algorithms: Adaptive-Mixture Metropolis (AMM), Adaptive Metropolis-within-Gibbs (AMWG), Automated Factor Slice Sampler (AFSS), Elliptical Slice Sampler (ESS), Hit-And-Run Metropolis (HARM), Metropolis-within-Gibbs (MWG), Random-Walk Metropolis (RWM), Robust Adaptive Metropolis (RAM), Slice Sampler (Slice), and the Univariate Eigenvector Slice Sampler (UESS).

Large-dimensional models often require blockwise sampling. For example, with thousands of parameters, a componentwise algorithm must evaluate the model specification function once per parameter per iteration, resulting in an algorithm that may take longer than is acceptable to produce samples. Algorithms that require derivatives, such as the family of Hamiltonian Monte Carlo (HMC), require even more evaluations of the model specification function per iteration, and quickly become too costly in large dimensions. Finally, algorithms with multivariate proposals often have difficulty producing an accepted proposal in large-dimensional models. The most practical solution is to group parameters into $N$ blocks, and each iteration the algorithm evaluates the model specification function $N$ times, each with a reduced set of parameters.

The Blocks function performs either a sequential assignment of parameters to blocks when posterior correlation is not supplied, or uses hierarchical clustering to create blocks based on posterior correlation. If posterior correlation is supplied, then the user may specify a range of the number of blocks to consider, and the optimal number of blocks is considered to be the maximum of the mean silhouette width of each hierarchical clustering. Silhouette width is calculated as per the cluster package. Hierarchical clustering is performed on the distance matrix calculated from the dissimilarity matrix (1 - abs(PostCor)) of the posterior correlation matrix. With sequential assignment, the number of parameters per block is approximately equal. With hierarchical clustering, the number of parameters per block may vary widely. Creating blocks from hierarchical clustering performs well in practice, though there are many alternative methods the user may consider outside of this function, such as using factor analysis, model-based clustering, or other methods.

Aside from sequentially-assigned blocks, or blocks based on posterior correlation, it is also common to group parameters with similar uses, such as putting regression effects parameters into one block, and autocorrelation parameters into another block. Another popular way to group parameters into blocks is by time-period for some time-series models. These alternative blocking strategies are unsupported in the Blocks function, and best left to user discretion.

Some MCMC algorithms that accept blocked parameters also require blocked variance-covariance matrices. The Blocks function does not return these matrices, because it may not be necessary, or when it is, the user may prefer identity matrices, scaled identity matrices, or matrices with explicitly-defined elements.

If the user is looking for a place to begin with blockwise sampling, then the recommended, default approach (when blocked parameters by time-period are not desired in a time-series) is to begin with a trial run of the adaptive, unblocked HARM algorithm (since covariance matrices are not required) for the purposes of obtaining a posterior correlation matrix. Next, create blocks with the Blocks function based on the posterior correlation matrix obtained from the trial run. Finally, run the desired, blocked algorithm with the newly created blocks (and possibly user-specified covariance matrices), beginning where the trial run ended.

If hierarchical clustering is used, then it is important to note that hierarchical clustering has no idea that the user intends to perform blockwise sampling in MCMC. If hierarchical clustering returns numerous small blocks, then the user may consider combining some or all of those blocks. For example, if several 1-parameter blocks are returned, then blockwise sampling will equal componentwise sampling for those blocks, which will iterate slower. Conversely, if hierarchical clustering returns one or more big blocks, each with enough parameters that multivariate sampling will have difficulty getting an accepted proposal, or an accepted proposal that moves more than a small amount, then the user may consider subdividing these big blocks into smaller, more manageable blocks, though with the understanding that more posterior correlation is unaccounted for.

Value

The Blocks function returns an object of class blocks, which is a list. Each component of the list is a block of parameters, and parameters are indicated by their position in the initial values vector.

Author(s)

Statisticat, LLC. [email protected]

See Also

LaplacesDemon and PosteriorChecks.

Examples

library(LaplacesDemon)

### Create the default number of sequentially assigned blocks:
Initial.Values <- rep(0,1000)
MyBlocks <- Blocks(Initial.Values)
MyBlocks

### Or, a pre-specified number of sequentially assigned blocks:
#Initial.Values <- rep(0,1000)
#MyBlocks <- Blocks(Initial.Values, N=20)

### If scaled diagonal covariance matrices are desired:
#VarCov <- list()
#for (i in 1:length(MyBlocks))
#  VarCov[[i]] <- diag(length(MyBlocks[[i]]))*2.38^2/length(MyBlocks[[i]])

### Or, determine the number of blocks in the range of 2 to 50 from
### hierarchical clustering on the posterior correlation matrix of an
### object, say called Fit, output from LaplacesDemon:
#MyBlocks <- Blocks(Initial.Values, N=c(2,50),
#  PostCor=cor(Fit$Posterior1))
#lapply(MyBlocks, length) #See the number of parameters per block

### Or, create a pre-specified number of blocks from hierarchical
### clustering on the posterior correlation matrix of an object,
### say called Fit, output from LaplacesDemon:
#MyBlocks <- Blocks(Initial.Values, N=20, PostCor=cor(Fit$Posterior1))

### Posterior correlation from a previous trial run could be obtained
### with either method below (though cor() will be fastest because
### additional checks are not calculated for the parameters):
#rho <- cor(Fit$Posterior1)
#rho <- PosteriorChecks(Fit)$Posterior.Correlation

Description

Given a matrix of posterior samples from MCMC, the BMK.Diagnostic function calculates Hellinger distances between consecutive batches for each chain. This is useful for monitoring convergence of MCMC chains.

Usage

BMK.Diagnostic(X, batches=10)

Arguments

Details

Hellinger distance is used to quantify dissimilarity between two probability distributions. It is based on the Hellinger integral, introduced by Hellinger (1909). Traditionally, Hellinger distance is bound to the interval [0,1], though another popular form occurs in the interval [0,$\sqrt{2}$]. A higher value of Hellinger distance is associated with more dissimilarity between the distributions.

Convergence is assumed when Hellinger distances are below a threshold, indicating that posterior samples are similar between consecutive batches. If all Hellinger distances beyond a given batch of samples is below the threshold, then burnin is suggested to occur immediately before the first batch of satisfactory Hellinger distances.

As an aid to interpretation, consider a matrix of 1,000 posterior samples from three chains: beta[1], beta[2], and beta[3]. With 10 batches, the column names are: 100, 200, ..., 900. A Hellinger distance for the chain beta[1] at 100 is the Hellinger distance between two batches: samples 1-100, and samples 101:200.

A benefit to using BMK.Diagnostic is that the resulting Hellinger distances may easily be plotted with the plotMatrix function, allowing the user to see quickly which consecutive batches of which chains were dissimilar. This makes it easier to find problematic chains.

The BMK.Diagnostic is calculated automatically in the LaplacesDemon function, and is one of the criteria in the Consort function regarding the recommendation of when to stop updating the Markov chain Monte Carlo (MCMC) sampler in LaplacesDemon.

For more information on the related topics of burn-in and stationarity, see the burnin and is.stationary functions, and the accompanying vignettes.

Value

The BMK.Diagnostic function returns an object of class bmk that is a $J \times B$ matrix of Hellinger distances between consecutive batches for $J$ parameters of posterior samples. The number of columns, $B$ is equal to the number of batches minus one.

The BMK.Diagnostic function is similar to the bmkconverge function in package BMK.

References

Boone, E.L., Merrick, J.R. and Krachey, M.J. (2013). "A Hellinger Distance Approach to MCMC Diagnostics". Journal of Statistical Computation and Simulation, in press.

Hellinger, E. (1909). "Neue Begrundung der Theorie quadratischer Formen von unendlichvielen Veranderlichen" (in German). Journal fur die reine und angewandte Mathematik, 136, p. 210–271.

See Also

burnin, Consort, is.stationary, and LaplacesDemon.

Examples

library(LaplacesDemon)
N <- 1000 #Number of posterior samples
J <- 10 #Number of parameters
Theta <- matrix(runif(N*J),N,J)
colnames(Theta) <- paste("beta[", 1:J, "]", sep="")
for (i in 2:N) {Theta[i,1] <- Theta[i-1,1] + rnorm(1)}
HD <- BMK.Diagnostic(Theta, batches=10)
plot(HD, title="Hellinger distance between batches")

Description

The burnin function estimates the duration of burn-in in iterations for one or more Markov chains. “Burn-in” refers to the initial portion of a Markov chain that is not stationary and is still affected by its initial value.

Usage

burnin(x, method="BMK")

Arguments

Details

Burn-in is a colloquial term for the initial iterations in a Markov chain prior to its convergence to the target distribution. During burn-in, the chain is not considered to have “forgotten” its initial value.

Burn-in is not a theoretical part of MCMC, but its use is the norm because of the need to limit the number of posterior samples due to computer memory. If burn-in were retained rather than discarded, then more posterior samples would have to be retained. If a Markov chain starts anywhere close to the center of its target distribution, then burn-in iterations do not need to be discarded.

In the LaplacesDemon function, stationarity is estimated with the BMK.Diagnostic function on all thinned posterior samples of each chain, beginning at cumulative 10% intervals relative to the total number of samples, and the lowest number in which all chains are stationary is considered the burn-in.

The term, “burn-in”, originated in electronics regarding the initial testing of component failure at the factory to eliminate initial failures (Geyer, 2011). Although “burn-in' has been the standard term for decades, some are referring to these as “warm-up” iterations.

Value

The burnin function returns a vector equal in length to the number of MCMC chains in x, and each element indicates the maximum iteration in burn-in.

Author(s)

Statisticat, LLC. [email protected]

References

Geyer, C.J. (2011). "Introduction to Markov Chain Monte Carlo". In S Brooks, A Gelman, G Jones, and M Xiao-Li (eds.), "Handbook of Markov Chain Monte Carlo", p. 3–48. Chapman and Hall, Boca Raton, FL.

See Also

BMK.Diagnostic, deburn, Geweke.Diagnostic, KS.Diagnostic, and LaplacesDemon.

Examples

library(LaplacesDemon)
x <- rnorm(1000)
burnin(x)

Description

A caterpillar plot is a horizontal plot of 3 quantiles of selected distributions. This may be used to produce a caterpillar plot of posterior samples (parameters and monitored variables) from an object either of class demonoid, demonoid.hpc, iterquad, laplace, pmc, vb, or a matrix.

Usage

caterpillar.plot(x, Parms=NULL, Title=NULL)

Arguments

Details

Caterpillar plots are popular plots in Bayesian inference for summarizing the quantiles of posterior samples. A caterpillar plot is similar to a horizontal boxplot, though without quartiles, making it easier for the user to study more distributions in a single plot. The following quantiles are plotted as a line for each parameter: 0.025 and 0.975, with the exception of a generic matrix, where unimodal 95% HPD intervals are estimated (for more information, see p.interval). A vertical, gray line is included at zero. For all but class demonoid.hpc, the median appears as a black dot, and the quantile line is black. For class demonoid.hpc, the color of the median and quantile line differs by chain; the first chain is black and additional chains appear beneath.

Author(s)

Statisticat, LLC. [email protected]

See Also

IterativeQuadrature, LaplaceApproximation, LaplacesDemon, LaplacesDemon.hpc, PMC, p.interval, SIR, and VariationalBayes.

Examples

#An example is provided in the LaplacesDemon function.

Description

This function either centers and scales a continuous variable and provides options for binary variables, or returns an untransformed variable from a centered and scaled variable.

Usage

CenterScale(x, Binary="none", Inverse=FALSE, mu, sigma, Range, Min)

Arguments

Details

Gelman (2008) recommends centering and scaling continuous predictors to facilitate MCMC convergence and enable comparisons between coefficients of centered and scaled continuous predictors with coefficients of untransformed binary predictors. A continuous predictor is centered and scaled as follows: x.cs <- (x - mean(x)) / (2*sd(x)). This is an improvement over the usual practice of standardizing predictors, which is x.z <- (x - mean(x)) / sd(x), where coefficients cannot be validly compared between binary and continuous predictors.

In MCMC, such as in LaplacesDemon, a centered and scaled predictor often results in a higher effective sample size (ESS), and therefore the chain mixes better. Centering and scaling is a method of re-parameterization to improve mixing.

Griffin and Brown (2013) also assert that the user may not want to scale predictors that are measured on the same scale, since scaling in this case may increase noisy, low signals. In this case, centering (without scaling) is recommended. To center a predictor, subtract its mean.

Value

The CenterScale function returns a centered and scaled vector, or the untransformed vector.

References

Gelman, A. (2008). "Scaling Regression Inputs by Dividing by Two Standard Devations". Statistics in Medicine, 27, p. 2865–2873.

Griffin, J.E. and Brown, P.J. (2013) "Some Priors for Sparse Regression Modelling". Bayesian Analysis, 8(3), p. 691–702.

See Also

ESS, IterativeQuadrature, LaplaceApproximation, LaplacesDemon, and PMC.

Examples

### See the LaplacesDemon function for an example in use.
library(LaplacesDemon)
x <- rnorm(100,10,1)
x.cs <- CenterScale(x)
x.orig <- CenterScale(x.cs, Inverse=TRUE, mu=mean(x), sigma=sd(x))

Description

This function combines objects of class demonoid.

Usage

Combine(x, Data, Thinning=1)

Arguments

Details

The purpose of the Combine function is to enable a user to combine objects of class demonoid for one of three reasons. First, parallel chains from LaplacesDemon.hpc may be combined after convergence is assessed with Gelman.Diagnostic. Second, consecutive updates of single chains from LaplacesDemon or parallel chains from LaplacesDemon.hpc may be combined when the computer has insufficient random-access memory (RAM) for the user to update once with enough iterations. Third, consecutive single-chain or parallel-chain updates may be combined when it seems that the logarithm of the joint posterior distribution, LP, seems to be oscillating up and down, which is described in more detail below.

The most common use regards the combination of parallel chains output from LaplacesDemon.hpc. Typically, a user with parallel chains examines them graphically with the caterpillar.plot and plot (actually, plot.demonoid) functions, and assesses convergence with the Gelman.Diagnostic function. Thereafter, the parallel chain output in the object of class demonoid.hpc should be combined into a single object of class demonoid, before doing posterior predictive checks and making inferences. In this case, the Thinning argument usually is recommended to remain at its default.

It is also common with a high-dimensional model (a model with a large number of parameters) to need more posterior samples than allowed by the random-access memory (RAM) of the computer. In this case, it is best to use the LaplacesDemon.RAM function to estimate the amount of RAM that a given model will require with a given number of iterations, and then update LaplacesDemon almost as much as RAM allows, and save the output object of class demonoid. Then, the user is advised to continue onward with a consecutive update (after using as.initial.values and anything else appropriate to prepare for the consecutive update). Suppose a user desires to update a gigantic model with thousands of parameters, and with the aid of LaplacesDemon.RAM, estimates that they can safely update only 100,000 iterations, and that 150,000 iterations would exceed RAM and crash the computer. The patient user can update several consecutive models, each with retaining only 1,000 thinned posterior samples, and combine them later with the Combine function, by placing multiple objects into a list, as described below. In this way, it is possible for a user to update models that otherwise far exceed computer RAM.

Less commonly, multiple updates of single-chain objects should be combined into a single object of class demonoid. This is most useful in complicated models that are run for large numbers of iterations, where it may be suspected that stationarity has been achieved, but that thinning is insufficient, and the samples may be combined and thinned. If followed, then these suggestions may continue seemingly to infinity, and the unnormalized logarithm of the joint posterior density, LP, may seem to oscillate, sometimes improving and getting higher, and getting lower during other updates. For this purpose, the prior covariance matrix of the last model is retained (rather than combining them). This may be an unpleasant surprise for combining parallel updates, so be aware of it.

In these cases, which usually involve complicated models with high autocorrelation in the chains, the user may opt to use parallel processing with the LaplacesDemon.hpc function, or may use the LaplacesDemon function as follows. The user should save (meaning, not overwrite) each object of class demonoid, place multiple objects into a list, and use the Combine function to combine these objects.

For example, suppose a user names the object Fit, as in the LaplacesDemon example. Now, rather than overwriting object Fit, object Fit is renamed, after updating a million iterations, to Fit1. As suggested by Consort, another million iterations are used, but now to create object Fit2. Further suppose this user specified Thinning=1000 in LaplacesDemon, meaning that the million iterations are thinned by 1,000, so only 1,000 iterations are retained in each object, Fit1 and Fit2. In this case, Combine combines the information in Fit1 and Fit2, and returns an object the user names Fit3. Fit3 has only 1,000 iterations, which is the result of appending the iterations in Fit1 and Fit2, and thinning by 2. If 2,000,000 iterations were updated from the beginning, and were thinned by 2,000, then the same information exists now in Fit3. The Consort function can now be applied to Fit3, to see if stationarity is found. If not, then more objects of class demonoid can be collected and combined.

Value

This function returns an object of class demonoid. For more information on an object of class demonoid, see the LaplacesDemon function.

Author(s)

Statisticat, LLC. [email protected]

See Also

caterpillar.plot, Gelman.Diagnostic, LaplacesDemon, LaplacesDemon.hpc, and Thin.

Description

This function provides several styles of conditional plots with base graphics.

Usage

cond.plot(x, y, z, Style="smoothscatter")

Arguments

Details

The cond.plot function provides simple conditional plots with base graphics. All plot styles are conditional upon z. Up to nine conditional plots are produced in a panel.

Plots include:

boxplot: y ~ x | z densover: f(x | z) hist: x | z scatter: x, y | z smoothscatter: x, y | z

The cond.plot function is not intended to try to compete with some of the better graphics packages, but merely to provide simple functionality.

Value

Conditional plots are returned.

Author(s)

Statisticat, LLC. [email protected]

See Also

joint.density.plot and joint.pr.plot.

Examples

library(LaplacesDemon)
x <- rnorm(1000)
y <- runif(1000)
z <- rcat(1000, rep(1/4,4))
cond.plot(x, y, z, Style="smoothscatter")

Description

This may be used to consort with Laplace's Demon regarding an object of class demonoid. Laplace's Demon will offer suggestions.

Usage

Consort(object)

Arguments

Details

First, Consort calls print.demonoid, which prints most of the components to the screen from the supplied object of class demonoid.

Second, Laplace's Demon considers a combination of five conditions when making the largest part of its suggestion. These conditions are: the algorithm, acceptance rate, MCSE, ESS, and stationarity. Other things are considered as well, such as the recommended thinning value is used to suggest a new number of iterations, how fast the algorithm is expected to be, and if the condition of diminishing adaptation (also called the vanishing adaptation condition) was met (for an adaptive algorithm). Diminishing adaptation occurs only when the absolute value of the proposed variances trends downward (toward zero) over the course of all adaptations. When an algorithm is adaptive and it does not have diminishing adaptations, the Consort function will suggest a different adaptive algorithm. The Periodicity argument is suggested to be set equal to the value of Rec.Thinning.

Appeasement applies only when all parameters are continuous.The Hangartner.Diagnostic should be considered for discrete parameters.

Appeasement Conditions

• Algorithm: The final algorithm must be non-adaptive, so that the Markov property holds. This is conservative. A user may have an adaptive (non-final) algorithm in which adaptations in the latest update are stationary, or no longer diminishing. Laplace's Demon is unaware of previous updates, and conservatively interprets this as failing to meet the condition of diminishing adaptation, when the output may be satisfactory. On the other hand, if the adaptive algorithm has essentially stopped adapting, and if there is a non-adaptive version, then the user should consider switching to the non-adaptive algorithm. User discretion is advised.

• Acceptance Rate: The acceptance rate is considered satisfactory if it is within the interval [15%,50%] for most algorithms. Some algorithms have different recommended intervals.

• MCSE: The Monte Carlo Standard Error (MCSE) is considered satisfactory for each target distribution if it is less than 6.27% of the standard deviation of the target distribution. This allows the true mean to be within 5% of the area under a Gaussian distribution around the estimated mean. The MCSE function is used. Toft et al. (2007) propose a stricter criterion of 5%. The criterion of 6.27% for this stopping rule is arbitrary, and may be too lenient or strict, depending on the needs of the user. Nonetheless, it has performed well, and this type of stopping rule has been observed to perform better than MCMC convergence diagnostics (Flegal et al., 2008).

• ESS: The effective sample size (ESS) is considered satisfactory for each target distribution if it is at least 100, which is usually enough to describe 95% probability intervals (see p.interval and LPL.interval for more information). The ESS function is used. When this criterion is unmet, the name of the worst mixing chain in Summary1 appears.

• Stationarity: Each target distribution is considered satisfactory if it is estimated to be stationary with the BMK.Diagnostic function.

Bear in mind that the MCSE, ESS, and stationarity criteria are all univariate measures applied to each marginal posterior distribution. Multivariate forms are not included. By chance alone due to multiple independent tests, 5% of these diagnostics should indicate non-convergence when 'convergence' exists. In contrast, even one non-convergent nuisance parameter is associated with non-convergence in all other parameters. Assessing convergence is difficult.

If all five conditions are satisfactory, then Laplace's Demon is appeased. Otherwise, Laplace's Demon will suggest and supply R code that is ready to be copy/pasted and executed.

To visualize the MCSE-based stopping rule, run the following code:

x <- seq(from=-3, to=3, by=0.1); plot(x, dnorm(x,0,1), type="l"); abline(v=-0.0627); abline(v=0.0627); abline(v=2*-0.0627, col="red"); abline(v=2*0.0627, col="red")

The black vertical lines show the standard error, and the red vertical lines show the 95% interval.

If the user has an object of class demonoid.hpc, then the Consort function may be still be applied, but a particular chain in the object must be specified as a component in a list. For example, with an object called Fit and a goal of consorting over the second chain, the code would be: Consort(Fit[[2]]).

The Demonic Suggestion is usually very helpful, but should not be followed blindly. Do not let it replace critical thinking. For example, Consort may find that diminishing adaptation is unmet, and recommend a different algorithm. However, the user may be convinced that the current algorithm is best, and believe instead that MCMC found a local solution, and is leaving it to find the global solution, in which case adaptations may increase again. Diminishing adaptation may have occurred in a previous run, and is not found in the current run because adaptation is essentially finished. If either of these is true, then it may be best to ignore the newly suggested algorithm, and continue with the current algorithm. The suggested code may be helpful, but it is merely a suggestion.

If achieving the appeasement of Laplace's Demon is difficult, consider ignoring the MCSE criterion and terminate when all other criteria have been met, placing special emphasis on ESS.

Author(s)

Statisticat, LLC. [email protected]

References

Flegal, J.M., Haran, M., and Jones, G.L. (2008). "Markov chain Monte Carlo: Can We Trust the Third Significant Figure?". Statistical Science, 23, p. 250–260.

Toft, N., Innocent, G., Gettinby, G., and Reid, S. (2007). "Assessing the Convergence of Markov Chain Monte Carlo Methods: An Example from Evaluation of Diagnostic Tests in Absence of a Gold Standard". Preventive Veterinary Medicine, 79, p. 244–256.

See Also

BMK.Diagnostic, ESS, Hangartner.Diagnostic, LaplacesDemon, LaplacesDemon.hpc, LPL.interval, MCSE, and p.interval.

Description

The Cumulative Sample Function (CSF) is a visual MCMC diagnostic in which the user may select a measure (such as a variable, summary statistic, or other diagnostic), and observe a plot of how the measure changes over cumulative posterior samples from MCMC, such as the output of LaplacesDemon. This may be considered to be a generalized extension of the cumuplot in the coda package, which is a more restrictive form of the cusum diagnostic introduced by Yu and Myckland (1998).

Yu and Myckland (1998) suggest that CSF plots should be examined after traditional trace plots seem convergent, and assert that faster mixing chains (which are more desirable) result in CSF plots that are more ‘hairy’ (as opposed to smooth), though this is subjective and has been debated. The LaplacesDemon package neither supports nor contradicts the suggestion of mixing and ‘hairiness’, but suggests that CSF plots may be used to provide additional information about a chain. For example, a user may decide on a practical burnin given when a conditional mean obtains a certain standard error.

Usage

CSF(x, name, method="Quantiles", quantiles=c(0.025,0.500,0.975), output=FALSE)

Arguments

Details

When method="ESS", the effective sample size (ESS) is observed as a function of the cumulative samples of x. For more information, see the ESS function.

When method="Geweke.Diagnostic", the Z-score output of the Geweke diagnostic is observed as a function of the cumulative samples of x. For more information, see the Geweke.Diagnostic function.

When method="HPD", the Highest Posterior Density (HPD) interval is observed as a function of the cumulative samples of x. For more information, see the p.interval function.

When method="is.stationary", stationarity is logically tested and the result is observed as a function of the cumulative samples of x. For more information, see the is.stationary function.

When method="Kurtosis", kurtosis is observed as a function of the cumulative samples of x.

When method="MCSE", the Monte Carlo Standard Error (MCSE) estimated with the IMPS method is observed as a function of the cumulative samples of x. For more information, see the MCSE function.

When method="MCSE.bm", the Monte Carlo Standard Error (MCSE) estimated with the batch.means method is observed as a function of the cumulative samples of x. For more information, see the MCSE function.

When method="MCSE.sv", the Monte Carlo Standard Error (MCSE) estimated with the sample.variance method is observed as a function of the cumulative samples of x. For more information, see the MCSE function.

When method="Mean", the mean is observed as a function of the cumulative samples of x.

When method="Mode", the estimated mode is observed as a function of the cumulative samples of x. For more information, see the Mode function.

When method="N.Modes", the estimated number of modes is observed as a function of the cumulative samples of x. For more information, see the Modes function.

When method="Precision", the precision (inverse variance) is observed as a function of the cumulative samples of x.

When method="Quantiles", the quantiles selected with the quantiles argument are observed as a function of the cumulative samples of x.

When method="Skewness", skewness is observed as a function of the cumulative samples of x.

Author(s)

Statisticat, LLC. [email protected]

References

Yu, B. and Myckland, P. (1997). "Looking at Markov Samplers through Cusum Path Plots: A Simple Diagnostic Idea". Statistics and Computing, 8(3), p. 275–286.

See Also

burnin, ESS, Geweke.Diagnostic, is.stationary, LaplacesDemon, MCSE, Mode, Modes, and p.interval.

Examples

#Commented-out because of run-time for package builds
#library(LaplacesDemon)
#x <- rnorm(1000)
#CSF(x, method="ESS")
#CSF(x, method="Geweke.Diagnostic")
#CSF(x, method="HPD")
#CSF(x, method="is.stationary")
#CSF(x, method="Kurtosis")
#CSF(x, method="MCSE")
#CSF(x, method="MCSE.bm")
#CSF(x, method="MCSE.sv")
#CSF(x, method="Mean")
#CSF(x, method="Mode")
#CSF(x, method="N.Modes")
#CSF(x, method="Precision")
#CSF(x, method="Quantiles")
#CSF(x, method="Skewness")

Description

This data set is for discrete choice models and consists of the choice of commuting route to school: arterial, two-lane, or freeway. There were 151 Pennsylvania commuters who started from a residential complex in State College, PA, and commute to downtown State College.

Usage

data(demonchoice)

Format

This data frame contains 151 rows of individual choices and 9 columns. The following data dictionary describes each variable or column.

Choice

This is the route choice: four-lane arterial (35 MPH speed limit), two-lane highway (35 MPH speed limit, with one lane in each direction), or a limited-access four-lane freeway (55 MPH speed liimit.)

HH.Income

This is an ordinal variable of annual household income of the commuter in USD. There are four categories: 1 is less than 20,000 USD, 2 is 20,000-29,999 USD, 3 is 30,000-39,999 USD, and 4 is 40,000 USD or greater.

Vehicle.Age

This is the age in years of the vehicle of the commuter.

Stop.Signs.Arterial

This is the number of stop signs along the arterial route.

Stop.Signs.Two.Lane

This is the number of stop signs along the two-lane route.

Stop.Signs.Freeway

This is the number of stop signs along the freeway route.

Distance.Arterial

This is distance in miles of the arterial route.

Distance.Two.Lane

This is the distance in miles of the two-lane route.

Distance.Freeway

This is the distance in miles of the freeway route.

Source

Washington, S., Congdon, P., Karlaftis, M., and Mannering, F. (2009). "Bayesian Multinomial Logit: Theory and Route Choice Example". Transportation Research Record, 2136, p. 28–36.

Description

This data set consists of daily currency pair prices from 2010 through 2014. Each currency pair has a close, high, and low price.

Usage

data(demonfx)

Format

This data frame contains 1,301 rows as time-periods (with row names) and 39 columns of currency pair prices. The following data dictionary describes each time-series or column.

EURUSD.Close

This is the currency pair closing price.

EURUSD.High

This is the currency pair high price.

EURUSD.Low

This is the currency pair low price.

USDJPY.Close

This is the currency pair closing price.

USDJPY.High

This is the currency pair high price.

USDJPY.Low

This is the currency pair low price.

USDCHF.Close

This is the currency pair closing price.

USDCHF.High

This is the currency pair high price.

USDCHF.Low

This is the currency pair low price.

GBPUSD.Close

This is the currency pair closing price.

GBPUSD.High

This is the currency pair high price.

GBPUSD.Low

This is the currency pair low price.

USDCAD.Close

This is the currency pair closing price.

USDCAD.High

This is the currency pair high price.

USDCAD.Low

This is the currency pair low price.

EURGBP.Close

This is the currency pair closing price.

EURGBP.High

This is the currency pair high price.

EURGBP.Low

This is the currency pair low price.

EURJPY.Close

This is the currency pair closing price.

EURJPY.High

This is the currency pair high price.

EURJPY.Low

This is the currency pair low price.

EURCHF.Close

This is the currency pair closing price.

EURCHF.High

This is the currency pair high price.

EURCHF.Low

This is the currency pair low price.

AUDUSD.Close

This is the currency pair closing price.

AUDUSD.High

This is the currency pair high price.

AUDUSD.Low

This is the currency pair low price.

GBPJPY.Close

This is the currency pair closing price.

GBPJPY.High

This is the currency pair high price.

GBPJPY.Low

This is the currency pair low price.

CHFJPY.Close

This is the currency pair closing price.

CHFJPY.High

This is the currency pair high price.

CHFJPY.Low

This is the currency pair low price.

GBPCHF.Close

This is the currency pair closing price.

GBPCHF.High

This is the currency pair high price.

GBPCHF.Low

This is the currency pair low price.

NZDUSD.Close

This is the currency pair closing price.

NZDUSD.High

This is the currency pair high price.

NZDUSD.Low

This is the currency pair low price.

Source

https://www.global-view.com/forex-trading-tools/forex-history/index.html

Description

These are the monthly number of user sessions at https://web.archive.org/web/20141224051720/http://www.bayesian-inference.com/index by continent. Additional data may be added in the future.

Usage

data(demonsessions)

Format

This data frame contains 26 rows (with row names) and 6 columns. The following data dictionary describes each variable or column.

Africa

This is the African continent.

Americas

This is North and South America.

Asia

This is the Asian continent.

Europe

This is Europe as a continent.

Oceania

This is Oceania, such as Australia.

Not.Set

This includes sessions in which the continent was not set, or is unknown.

Source

https://web.archive.org/web/20141224051720/http://www.bayesian-inference.com/index

Description

Late one night, after witnessing Laplace's Demon in action, I followed him back to what seemed to be his lair. Minutes later, he left again. I snuck inside and saw something labeled 'Demon Snacks'. Hurriedly, I recorded the 39 items, each with a name and 10 nutritional attributes.

Usage

data(demonsnacks)

Format

This data frame contains 39 rows (with row names) and 10 columns. The following data dictionary describes each variable or column.

Serving.Size

This is serving size in grams.

Calories

This is the number of calories.

Total.Fat

This is total fat in grams.

Saturated.Fat

This is saturated fat in grams.

Cholesterol

This is cholesterol in milligrams.

Sodium

This is sodium in milligrams.

Total.Carbohydrate

This is the total carbohydrates in grams.

Dietary.Fiber

This is dietary fiber in grams.

Sugars

This is sugar in grams.

Protein

This is protein in grams.

Source

This data was obtained from the lair of Laplace's Demon!

Description

This data set is for space-time models that require latitude and longitude, or coordinates. This data set consists of the minimum, mean, and maximum temperatures in Texas for 13 months.

Usage

data(demontexas)

Format

This data frame contains 369 rows of sites in Texas and 43 columns. The following data dictionary describes each variable or column.

Elevation

This is the elevation of the site.

Latitude

This is the latitude of the site.

Longitude

This is the longitude of the site.

Gulf

This is a gulf indicator of the site.

Max1

This is the maximum temperature in month 1.

Max2

This is the maximum temperature in month 2.

Max3

This is the maximum temperature in month 3.

Max4

This is the maximum temperature in month 4.

Max5

This is the maximum temperature in month 5.

Max6

This is the maximum temperature in month 6.

Max7

This is the maximum temperature in month 7.

Max8

This is the maximum temperature in month 8.

Max9

This is the maximum temperature in month 9.

Max10

This is the maximum temperature in month 10.

Max11

This is the maximum temperature in month 11.

Max12

This is the maximum temperature in month 12.

Max13

This is the maximum temperature in month 13.

Mean1

This is the mean temperature in month 1.

Mean2

This is the mean temperature in month 2.

Mean3

This is the mean temperature in month 3.

Mean4

This is the mean temperature in month 4.

Mean5

This is the mean temperature in month 5.

Mean6

This is the mean temperature in month 6.

Mean7

This is the mean temperature in month 7.

Mean8

This is the mean temperature in month 8.

Mean9

This is the mean temperature in month 9.

Mean10

This is the mean temperature in month 10.

Mean11

This is the mean temperature in month 11.

Mean12

This is the mean temperature in month 12.

Mean13

This is the mean temperature in month 13.

Min1

This is the minimum temperature in month 1.

Min2

This is the minimum temperature in month 2.

Min3

This is the minimum temperature in month 3.

Min4

This is the minimum temperature in month 4.

Min5

This is the minimum temperature in month 5.

Min6

This is the minimum temperature in month 6.

Min7

This is the minimum temperature in month 7.

Min8

This is the minimum temperature in month 8.

Min9

This is the minimum temperature in month 9.

Min10

This is the minimum temperature in month 10.

Min11

This is the minimum temperature in month 11.

Min12

This is the minimum temperature in month 12.

Min13

This is the minimum temperature in month 13.

Source

http://www.stat.ufl.edu/~winner/datasets.html

Description

The de.Finetti.Game function estimates the interval of a subjective probability regarding a possible event in the near future.

Usage

de.Finetti.Game(width)

Arguments

Details

This function is a variation on the game introduced by de Finetti, who is one of the main developers of subjective probability, along with Ramsey and Savage. In the original context, de Finetti proposed a gamble regarding life on Mars one billion years ago.

The frequentist interpretation of probability defines the probability of an event as the limit of its relative frequency in a large number of trials. Frequentist inference is undefined, for example, when there are no trials from which to calculate a probability. By defining probability relative to frequencies of physical events, frequentists attempt to objectify probability. However, de Finetti asserts that the frequentist (or objective) interpretation always reduces to a subjective interpretation of probability, because probability is a human construct and does not exist independently of humans in nature. Therefore, probability is a degree of belief, and is called subjective or personal probability.

Value

The de.Finetti.Game function returns a vector of length two. The respective elements are the lower and upper bounds of the subjective probability of the participant regarding the possible event in the near future.

Author(s)

Statisticat, LLC. [email protected]

See Also

elicit

Description

The deburn function discards or removes a user-specified number of burn-in iterations from an object of class demonoid.

Usage

deburn(x, BurnIn=0)

Arguments

Details

Documentation for the burnin function provides an introduction to the concept of burn-in as it relates to Markov chains.

The deburn function discards a number of the first posterior samples, as specified by the BurnIn argument. Stationarity is not checked, because it is assumed the user has a reason for using the deburn function, rather than using the results from the object of class demonoid. Therefore, the posterior samples in Posterior1 and Posterior2 are identical, as are Summary1 and Summary2.

Value

The deburn function returns an object of class demonoid.

Author(s)

Statisticat, LLC. [email protected]

See Also

burnin and LaplacesDemon.

Examples

### Assuming the user has Fit which is an object of class demonoid:
#library(LaplacesDemon)
#Fit2 <- deburn(Fit, BurnIn=100)

Description

These functions provide the density, distribution function, quantile function, and random generation for the univariate, asymmetric Laplace distribution with location parameter location, scale parameter scale, and asymmetry or skewness parameter kappa.

Usage

dalaplace(x, location=0, scale=1, kappa=1, log=FALSE)
palaplace(q, location=0, scale=1, kappa=1)
qalaplace(p, location=0, scale=1, kappa=1)
ralaplace(n, location=0, scale=1, kappa=1)

Arguments

Details

• Application: Continuous Univariate

• Density: $p(\theta) = \frac{\kappa \sqrt{2}}{\lambda (1+\kappa^2)} \exp(-|\theta-\mu| \frac{\sqrt{2}}{\lambda} \kappa^{|\theta-\mu|} |\theta-\mu|)$

• Inventor: Kotz, Kozubowski, and Podgorski (2001)

• Notation 1: $\theta \sim \mathcal{AL}(\mu, \lambda, \kappa)$

• Notation 2: $p(\theta) = \mathcal{AL}(\theta | \mu, \lambda, \kappa)$

• Parameter 1: location parameter $\mu$

• Parameter 2: scale parameter $\lambda > 0$

• Parameter 3: skewness parameter $\kappa > 0$

• Mean: $E(\theta) = \mu + \lambda \frac{1/\kappa - \kappa}{\sqrt{2}}$

• Variance: $var(\theta) = \lambda^2 \frac{1 + \kappa^4}{2 \kappa^2}$

• Mode: $mode(\theta) = \mu$

The asymmetric Laplace of Kotz, Kozubowski, and Podgorski (2001), also referred to as AL, is an extension of the univariate, symmetric Laplace distribution to allow for skewness. It is parameterized according to three parameters: location parameter $\mu$, scale parameter $\lambda$, and asymmetry or skewness parameter $\kappa$. The special case of $\kappa=1$ is the symmetric Laplace distribution. Values of $\kappa$ in the intervals $(0, 1)$ and $(1, \infty)$, correspond to positive (right) and negative (left) skewness, respectively. The AL distribution is leptokurtic, and its kurtosis ranges from 3 to 6 as $\kappa$ ranges from 1 to infinity. The skewness of the AL has been useful in engineering and finance. As an example, the AL distribution has been used as a replacement for Gaussian-distributed GARCH residuals. There is also an extension to the asymmetric multivariate Laplace distribution.

The asymmetric Laplace distribution is demonstrated in Kozubowski and Podgorski (2001) to be well-suited for financial modeling, specifically with currency exchange rates.

These functions are similar to those in the VGAM package.

Value

dalaplace gives the density, palaplace gives the distribution function, qalaplace gives the quantile function, and ralaplace generates random deviates.

References

Kotz, S., Kozubowski, T.J., and Podgorski, K. (2001). "The Laplace Distribution and Generalizations: a Revisit with Applications to Communications, Economics, Engineering, and Finance". Boston: Birkhauser.

Kozubowski, T.J. and Podgorski, K. (2001). "Asymmetric Laplace Laws and Modeling Financial Data". Mathematical and Computer Modelling, 34, p. 1003-1021.

See Also

dlaplace and dallaplace

Examples

library(LaplacesDemon)
x <- dalaplace(1,0,1,1)
x <- palaplace(1,0,1,1)
x <- qalaplace(0.5,0,1,1)
x <- ralaplace(100,0,1,1)

#Plot Probability Functions
x <- seq(from=-5, to=5, by=0.1)
plot(x, dalaplace(x,0,1,0.5), ylim=c(0,1), type="l", main="Probability Function",
     ylab="density", col="red")
lines(x, dalaplace(x,0,1,1), type="l", col="green")
lines(x, dalaplace(x,0,1,5), type="l", col="blue")
legend(1, 0.9, expression(paste(mu==0, ", ", lambda==1, ", ", kappa==0.5),
     paste(mu==0, ", ", lambda==1, ", ", kappa==1),
     paste(mu==0, ", ", lambda==1, ", ", kappa==5)),
     lty=c(1,1,1), col=c("red","green","blue"))

Description

These functions provide the density, distribution function, quantile function, and random generation for the univariate, asymmetric, log-Laplace distribution with location parameter $\mu$, scale parameter $\lambda$, and asymmetry or skewness parameter $\kappa$.

Usage

dallaplace(x, location=0, scale=1, kappa=1, log=FALSE)
pallaplace(q, location=0, scale=1, kappa=1)
qallaplace(p, location=0, scale=1, kappa=1)
rallaplace(n, location=0, scale=1, kappa=1)

Arguments

Details

• Application: Continuous Univariate

• Density 1: $p(\theta) = \exp(-\mu)\frac{(\sqrt(2)\kappa / \lambda)(\sqrt(2) / \lambda\kappa)}{(\sqrt(2)\kappa / \lambda)+(\sqrt(2) / (\lambda\kappa))} \exp(-(\frac{\sqrt(2)\kappa}{\lambda})+1), \quad \theta \ge \exp(\mu)$

• Density 2: $p(\theta) = \exp(-\mu) \frac{(\sqrt(2)\kappa / \lambda) (\sqrt(2) / (\lambda\kappa))}{(\sqrt(2)\kappa / \lambda) + (\sqrt(2) / (\lambda\kappa))} \exp(\frac{\sqrt(2)(\log(\theta)-\mu)}{\lambda\kappa} - (\log(\theta)-\mu)), \quad \theta < \exp(\mu)$

• Inventor: Pierre-Simon Laplace

• Notation 1: $\theta \sim \mathcal{ALL}(\mu, \lambda, \kappa)$

• Notation 2: $p(\theta) = \mathcal{ALL}(\theta | \mu, \lambda, \kappa)$

• Parameter 1: location parameter $\mu$

• Parameter 2: scale parameter $\lambda > 0$

• Mean: $E(\theta) =$

• Variance: $var(\theta) =$

• Mode: $mode(\theta) =$

The univariate, asymmetric log-Laplace distribution is derived from the Laplace distribution. Multivariate and symmetric versions also exist.

These functions are similar to those in the VGAM package.

Value

dallaplace gives the density, pallaplace gives the distribution function, qallaplace gives the quantile function, and rallaplace generates random deviates.

References

Kozubowski, T. J. and Podgorski, K. (2003). "Log-Laplace Distributions". International Mathematical Journal, 3, p. 467–495.

See Also

dalaplace, dexp, dlaplace, dlaplacep, dllaplace, dmvl, dnorm, dnormp, dnormv.

Examples

library(LaplacesDemon)
x <- dallaplace(1,0,1,1)
x <- pallaplace(1,0,1,1)
x <- qallaplace(0.5,0,1,1)
x <- rallaplace(100,0,1,1)

#Plot Probability Functions
x <- seq(from=0.1, to=10, by=0.1)
plot(x, dallaplace(x,0,1,0.5), ylim=c(0,1), type="l", main="Probability Function",
     ylab="density", col="red")
lines(x, dallaplace(x,0,1,1), type="l", col="green")
lines(x, dallaplace(x,0,1,5), type="l", col="blue")
legend(5, 0.9, expression(paste(mu==0, ", ", lambda==1, ", ", kappa==0.5),
     paste(mu==0, ", ", lambda==1, ", ", kappa==1),
     paste(mu==0, ", ", lambda==1, ", ", kappa==5)),
     lty=c(1,1,1), col=c("red","green","blue"))

Description

These functions provide the density and random generation for the asymmetric multivariate Laplace distribution with location and skew parameter $\mu$ and covariance $\Sigma$.

Usage

daml(x, mu, Sigma, log=FALSE)
raml(n, mu, Sigma)

Arguments

Details

• Application: Continuous Multivariate

• Density: $p(\theta) = \frac{2\exp(\theta\Omega\theta)}{(2\pi)^{k/2}|\Sigma|^0.5} \frac{\theta\Omega\theta}{2 + \mu\Omega\mu}^{(2-k)/4} K_{(2-k)/2}(\sqrt{(2 + \mu\Omega\mu)(\theta\Omega\theta)})$

• Inventor: Kotz, Kozubowski, and Podgorski (2003)

• Notation 1: $\theta \sim \mathcal{AL}_K(\mu, \Sigma)$

• Notation 2: $p(\theta) = \mathcal{AL}_K(\theta | \mu, \Sigma)$

• Parameter 1: location-skew parameter $\mu$

• Parameter 2: positive-definite covariance matrix $\Sigma$

• Mean: Unknown

• Variance: Unknown

• Mode: $mode(\theta) = \mu$

The asymmetric multivariate Laplace distribution of Kotz, Kozubowski, and Podgorski (2003) is a multivariate extension of the univariate, asymmetric Laplace distribution. It is parameterized according to two parameters: location-skew parameter $\mu$ and positive-definite covariance matrix $\Sigma$. Location and skew occur in the same parameter. When $\mu=0$, the density is the (symmetric) multivariate Laplace of Anderson (1992). As each location deviates from zero, the marginal distribution becomes more skewed. Since location and skew are combined, it is appropriate for zero-centered variables, such as a matrix of centered and scaled dependent variables in cluster analysis, factor analysis, multivariate regression, or multivariate time-series.

The asymmetric multivariate Laplace distribution is also discussed earlier in Kozubowski and Podgorski (2001), and is well-suited for financial modeling via multivariate regression, specifically with currency exchange rates. Cajigas and Urga (2005) fit residuals in a multivariate GARCH model with the asymmetric multivariate Laplace distribution, regarding stocks and bonds. They find that it "overwhelmingly outperforms" normality.

Value

daml gives the density, and raml generates random deviates.

References

Anderson, D.N. (1992). "A Multivariate Linnik Distribution". Statistical Probability Letters, 14, p. 333–336.

Cajigas, J.P. and Urga, G. (2005) "Dynamic Conditional Correlation Models with Asymmetric Laplace Innovations". Centre for Economic Analysis: Cass Business School.

Kotz, S., Kozubowski, T.J., and Podgorski, K. (2003). "An Asymmetric Multivariate Laplace Distribution". Working Paper.

Kozubowski, T.J. and Podgorski, K. (2001). "Asymmetric Laplace Laws and Modeling Financial Data". Mathematical and Computer Modelling, 34, p. 1003–1021.

See Also

dalaplace and dmvl

Examples

library(LaplacesDemon)
x <- daml(c(1,2,3), c(0,1,2), diag(3))
X <- raml(1000, c(0,1,2), diag(3))
joint.density.plot(X[,1], X[,2], color=FALSE)

Description

These functions provide the density, distribution function, quantile function, and random generation for the Bernoulli distribution.

Usage

dbern(x, prob, log=FALSE)
pbern(q, prob, lower.tail=TRUE, log.p=FALSE)
qbern(p, prob, lower.tail=TRUE, log.p=FALSE)
rbern(n, prob)

Arguments

Details

• Application: Continuous Univariate

• Density: $p(\theta) = {p}^{\theta} {(1-p)}^{1-\theta}$, $\theta = 0,1$

• Inventor: Jacob Bernoulli

• Notation 1: $\theta \sim \mathcal{BERN}(p)$

• Notation 2: $p(\theta) = \mathcal{BERN}(\theta | p)$

• Parameter 1: probability parameter $0 \le p \le 1$

• Mean: $E(\theta) = p$

• Variance: $var(\theta) = \frac{p}{1-p}$

• Mode: $mode(\theta) =$

The Bernoulli distribution is a binomial distribution with $n=1$, and one instance of a Bernoulli distribution is called a Bernoulli trial. One coin flip is a Bernoulli trial, for example. The categorical distribution is the generalization of the Bernoulli distribution for variables with more than two discrete values. The beta distribution is the conjugate prior distribution of the Bernoulli distribution. The geometric distribution is the number of Bernoulli trials needed to get one success.

Value

dbern gives the density, pbern gives the distribution function, qbern gives the quantile function, and rbern generates random deviates.

See Also

dbinom

Examples

library(LaplacesDemon)
dbern(1, 0.7)
rbern(10, 0.5)

Description

This is the density and random deviates function for the categorical distribution with probabilities parameter $p$.

Usage

dcat(x, p, log=FALSE)
qcat(pr, p, lower.tail=TRUE, log.pr=FALSE)
rcat(n, p)

Arguments

Details

• Application: Discrete Univariate

• Density: $p(\theta) = \sum \theta p$

• Inventor: Unknown (to me, anyway)

• Notation 1: $\theta \sim \mathcal{CAT}(p)$

• Notation 2: $p(\theta) = \mathcal{CAT}(\theta | p)$

• Parameter 1: probabilities $p$

• Mean: $E(\theta)$ = Unknown

• Variance: $var(\theta)$ = Unknown

• Mode: $mode(\theta)$ = Unknown

Also called the discrete distribution, the categorical distribution describes the result of a random event that can take on one of $k$ possible outcomes, with the probability $p$ of each outcome separately specified. The vector $p$ of probabilities for each event must sum to 1. The categorical distribution is often used, for example, in the multinomial logit model. The conjugate prior is the Dirichlet distribution.

Value

dcat gives the density and rcat generates random deviates.

Author(s)

Statisticat, LLC. [email protected]

See Also

as.indicator.matrix, ddirichlet, and dmultinom.

Examples

library(LaplacesDemon)
dcat(x=1, p=c(0.3,0.3,0.4))
rcat(n=10, p=c(0.1,0.3,0.6))