Chapter 6 The amt-library
This section gives an overview of the different probability density/mass functions, in addition to some utility functions available in the amt-library. The library provides an automatic
implementation of the methodology of Kleppe (2022b) for automatically computing a “metric tensors” suitable Riemannian sampling algorithms for arbitrary models.
The design of the library allows the user to use the same model specification for both Riemannian and regular fixed metric sampling (see Chapter 7).
The methodology of Kleppe (2022b) generally requires that all probability densities involved in the model where the argument depends on \(\boldsymbol \theta\) (generally the conditional densities of priors and latent variables) to have continuous first order derivatives. To achieve such continuity, many of the below distributions are known probability distributions that have been transformed to have support on the whole real line in such a way that the resulting probability distribution is sufficiently smooth. These transformations are also desirable from a numerical point of view.
As a basic example of such transformations involves the \(ExpGamma\)-distribution (see below for proper definition):
// inside operator()()
PARAMETER_SCALAR(logSigma,0.0);
model__+=expGamma_ld(logSigma,1.0,1.0);
varType sigma=exp(logSigma);
model__.generated(asDouble(sigma),"sigma");In the above example, sigma will be exponentially distributed with mean 1 (until further conditioning is done) and will be recorded as a generated quantity.
6.1 Utilities for symmetric positive definite matrices
The amt-package provides a set of utilities for representing symmetric positive definite (SPD) matrices. Internally, this class template SPDmatrix represents a \(d \times d\) SPD matrix as an “internal” vector of \(d(d+1)/2\) varTypes.
SPDmatrix objects may be used as e.g. covariance- or precision matrices in multivariate normal distributions, and it also possible to assign e.g. Wishart distributions on such objects (though in practice, such assignments are really distributions on the internal vector, consistent with the matrix-valued distribution).
The easiest way to construct a SPDmatrix is via the PARAMETER_SPD_MATRIX(name,dim,...) macro. Typical usage involves something like:
// inside operator()() :
PARAMETER_SPD_MATRIX(P,3);
// Note: PARAMETER_SPD_MATRIX(P,3) really makes a parameter vector
// named P_internal, and subsequently constructs SPDmatrix<varType> P
Eigen::VectorXd scaleDiag; scaleDiag.setConstant(3,4.0);
model__ += wishartDiagScale_ld(P,scaleDiag,10.0);
// P has a Wishart(4*I,10.0) distribution a priori
model__.generated(P,"P"); // full matrix stored
model__.generated(P.coeff(0,2),"P13"); // single element storedA complete class where the above code block is the body of the operator()()-method is stored in the file wishart_example.cpp, which is built an run via:
mdl <- pdmphmc::build("wishart_example.cpp")## model name : model## process type : HMCProcessConstr## Runge Kutta step type : RKDP54## Transport map type : diagLinearTM_VARI## compilation exited successfullyfit <- pdmphmc::run(mdl)
pdmphmc::clean.model(mdl,remove.all = TRUE)## [1] TRUEfit## run output for model: model## # of chains : 4## Summary based on discrete samples:## mean se_mean sd n_eff Rhat
## P_internal[1] 3.578 0.010 0.485 2492 1.005
## P_internal[2] 3.466 0.010 0.497 2431 1.003
## P_internal[3] 3.330 0.009 0.539 3209 1.002
## P_internal[4] 0.000 0.007 0.370 3153 1.002
## P_internal[5] -0.009 0.007 0.365 3057 1.001
## P_internal[6] 0.006 0.008 0.378 2498 1.001
## P[1,1] 39.964 0.364 18.438 2565 1.005
## P[2,1] 0.124 0.236 13.323 3176 1.002
## P[3,1] -0.229 0.236 13.037 3037 1.000
## P[2,2] 40.288 0.353 17.848 2559 1.005
## P[3,2] 0.341 0.245 12.695 2680 1.001
## P[3,3] 40.141 0.312 18.027 3330 1.001
## P13 -0.229 0.236 13.037 3037 1.000## summary based on integrated samples:## estimate se_estimate n_eff Rhat
## P[1,1] 39.673 0.190 1690 1.006
## P[2,1] 0.025 0.054 7448 1.040
## P[3,1] 0.052 0.054 6768 1.033
## P[2,2] 40.231 0.182 1733 1.004
## P[3,2] 0.156 0.110 1972 1.003
## P[3,3] 40.250 0.166 2207 1.002
## P13 0.052 0.054 6768 1.033## NOTE: integrated samples do NOT reflect the complete target distibution, only indicated moments with respect to the target distributionNotice that by default, the “internal” parameter vector P_internal is stored3 in addition to the generated quantities. Note, in this case, we would expect \(E(\mathbf P)=40 \mathbf I\).
6.2 Univariate continuous distributions
6.2.1 expGamma_ld(arg,shape,scale)
The distribution obtained by applying the (natural) logarithm to a gamma-distributed random variable with shape parameter \(\alpha=\)shape and scale parameter \(\beta=\)scale. Argument of resulting PDF \(x=\)arg. Log-density:
\[
\log p(x|\alpha,\beta) = -a\log(b) - \log(\Gamma(a)) + ax - \exp(x)/b.
\]
6.2.2 invLogitBeta_ld(x,a,b)
The distribution obtained by applying the logit-function to a Beta distributed random variable with shape parameters a and b. Argument of resulting PDF: \(x=\)x. Log-density:
\[
\log p(x|a,b) = \log(\Gamma(a+b)) - \log(\Gamma(a)\Gamma(b)) + a\log\left(\frac{\exp(x)}{1+\exp(x)}\right) - b\log(1+\exp(x)).
\]
6.3 Univariate discrete distributions
6.3.1 bernoulli_logit_lm(y,alpha)
Bernoulli-distribution with \(P(y=1)=\exp(\alpha)/(1+\exp(\alpha))\). Argument \(y\) is of integer type.
6.4 Distributions on unbounded vectors
6.4.1 multi_normal_prec_ld(arg,mu,Prec)
Multivariate normal distribution where the precision matrix Prec is a SPDmatrix (see Section 6.1). Both the argument arg and mean are vectors with dimension equal to Prec.dim().
6.4.2 iid_multi_normal_prec_ld(arg,mu,Prec)
Same as above, but arg is now a matrix, where the columns are iid realizations from \(N(\mu,P^{-1})\), where \(P=\)Prec.
6.4.3 multi_normal_ld(arg,mu,Sigma)
Multivariate normal distribution where the covariance matrix Sigma is a SPDmatrix (see Section 6.1). Both the argument arg and mean are vectors with dimension equal to Sigma.dim().
To do: iid-variant of this.
6.4.4 normalAR1_ld(arg,mu,phi,sigma)
Stationary Gaussian first order autoregressive process
\[
x_1 \sim N\left (\mu,\frac{\sigma^2}{1-\phi^2} \right), x_{t} = \mu + \phi ( x_{t-1}-\mu) + N(0,\sigma^2),t=2,\dots,T\;
\]
for scalar parameters \(\mu\), \(-1<\phi<1\) and \(\sigma>0\). \(\mathbf x=\)arg is a \(T\)-dimensional vector.
6.4.5 normalRW1_ld(x,sigma)
First order (intrinsic) Gaussian random walk model:
\[
x_t \sim N(x_{t-1},\sigma^2),\; t=2,\dots,T,
\]
for \(\sigma>0\). \(\mathbf x=\)x is a \(T\)-dimensional vector. Note that this model specifies a degenerate probability distribution, and \(\mathbf x\) cannot be sampled without further conditioning.
6.5 Distributions on symmetric positive definite matrices
The notation for the Wishart distribution is \(\mathcal W(\mathbf M,\nu)\) where \(\mathbf M\) is the scale matrix and \(\nu\) is the degrees of freedom parameter. Consequently, if \(\mathbf P \sim \mathcal W(\mathbf M,\nu)\), then \(E(\mathbf P)=\nu \mathbf M\).
References
This may be avoided by using the
store.pars-option inpdmphmc::run().↩︎