Parameters of conditional multivariate normal distribution

Description

This function calculates mean (expectation) and covariance matrix of conditional multivariate normal distribution.

Usage

cmnorm(
  mean,
  sigma,
  given_ind,
  given_x,
  dependent_ind = numeric(),
  is_validation = TRUE,
  is_names = TRUE,
  control = NULL,
  n_cores = 1L
)

Arguments

Details

Consider m-dimensional multivariate normal vector X=(X_{1},...,X_{m})^{T}~\sim N(\mu,\Sigma), where E(X)=\mu and Cov(X)=\Sigma are expectation (mean) and covariance matrix respectively.

Let's denote vectors of indexes of conditioned and unconditioned elements of X by I_{g} and I_{d} respectively. By x^{(g)} denote deterministic (column) vector of given values of X_{I_{g}}. The function calculates expected value and covariance matrix of conditioned multivariate normal vector X_{I_{d}} | X_{I_{g}} = x^{(g)}. For example if I_{g}=(1, 3) and x^{(g)}=(-1, 1) then I_{d}=(2, 4, 5) so the function calculates:

\mu_{c}=E\left(\left(X_{2}, X_{4}, X_{5}\right) | X_{1} = -1, X_{3} = 1\right)

\Sigma_{c}=Cov\left(\left(X_{2}, X_{4}, X_{5}\right) | X_{1} = -1, X_{3} = 1\right)

In general case:

\mu_{c} = E\left(X_{I_{d}} | X_{I_{g}} = x^{(g)}\right) = \mu_{I_{d}} + \left(x^{(g)} - \mu_{I_{g}}\right) \left(\Sigma_{(I_{d}, I_{g})} \Sigma_{(I_{g}, I_{g})}^{-1}\right)^{T}

\Sigma_{c} = Cov\left(X_{I_{d}} | X_{I_{g}} = x^{(g)}\right) = \Sigma_{(I_{d}, I_{d})} - \Sigma_{(I_{d}, I_{g})} \Sigma_{(I_{g}, I_{g})}^{-1} \Sigma_{(I_{g}, I_{d})}

Note that \Sigma_{(A, B)}, where A,B\in\{d, g\}, is a submatrix of \Sigma generated by intersection of I_{A} rows and I_{B} columns of \Sigma.

Below there is a correspondence between aforementioned theoretical (mathematical) notations and function arguments:

• mean - \mu.

• sigma - \Sigma.

• given_ind - I_{g}.

• given_x - x^{(g)}.

• dependent_ind - I_{d}

Moreover \Sigma_{(I_{g}, I_{d})} is a theoretical (mathematical) notation for sigma[given_ind, dependent_ind]. Similarly \mu_{g} represents mean[given_ind].

By default dependent_ind contains all indexes that are not in given_ind. It is possible to omit and duplicate indexes of dependent_ind. But at least single index should be provided for given_ind without any duplicates. Also dependent_ind and given_ind should not have the same elements. Moreover given_ind should not be of the same length as mean so at least one component should be unconditioned.

If given_x is a vector then (if possible) it will be treated as a matrix with the number of columns equal to the length of mean.

Currently control has no input arguments intended for the users. This argument is used for some internal purposes of the package.

Value

This function returns an object of class "mnorm_cmnorm".

An object of class "mnorm_cmnorm" is a list containing the following components:

• mean - conditional mean.

• sigma - conditional covariance matrix.

• sigma_d - covariance matrix of unconditioned elements.

• sigma_g - covariance matrix of conditioned elements.

• sigma_dg - matrix of covariances between unconditioned and conditioned elements.

• s12s22 - equals to the matrix product of sigma_dg and solve(sigma_g).

Note that mean corresponds to \mu_{c} while sigma represents \Sigma_{c}. Moreover sigma_d is \Sigma_{I_{d}, I_{d}}, sigma_g is \Sigma_{I_{g}, I_{g}} and sigma_dg is \Sigma_{I_{d}, I_{g}}.

Since \Sigma_{c} do not depend on X^{(g)} the output sigma does not depend on given_x. In particular output sigma remains the same independent of whether given_x is a matrix or vector. Oppositely if given_x is a matrix then output mean is a matrix which rows correspond to conditional means associated with given values provided by corresponding rows of given_x.

The order of elements of output mean and output sigma depends on the order of dependet_ind elements that is ascending by default. The order of given_ind elements does not matter. But, please, check that the order of given_ind match the order of given values i.e. the order of given_x columns.

Examples

# Consider multivariate normal vector:
# X = (X1, X2, X3, X4, X5) ~ N(mean, sigma)

# Prepare multivariate normal vector parameters
  # expected value
mean <- c(-2, -1, 0, 1, 2)
n_dim <- length(mean)
  # correlation matrix
cor <- c(   1,  0.1,  0.2,   0.3,  0.4,
          0.1,    1, -0.1,  -0.2, -0.3,
          0.2, -0.1,    1,   0.3,  0.2,
          0.3, -0.2,  0.3,     1, -0.05,
          0.4, -0.3,  0.2, -0.05,     1)
cor <- matrix(cor, ncol = n_dim, nrow = n_dim, byrow = TRUE)
  # covariance matrix
sd_mat <- diag(c(1, 1.5, 2, 2.5, 3))
sigma <- sd_mat %*% cor %*% t(sd_mat)

# Estimate parameters of conditional distribution i.e.
# when the first and the third components of X are conditioned:
# (X2, X4, X5 | X1 = -1, X3 = 1)
given_ind <- c(1, 3)
given_x <- c(-1, 1)
par <- cmnorm(mean = mean, sigma = sigma,
              given_ind = given_ind,
              given_x = given_x)
  # E(X2, X4, X5 | X1 = -1, X3 = 1)
par$mean
  # Cov(X2, X4, X5 | X1 = -1, X3 = 1)
par$sigma

# Additionally calculate E(X2, X4, X5 | X1 = 2, X3 = 3)
given_x_mat <- rbind(given_x, c(2, 3))
par1 <- cmnorm(mean = mean, sigma = sigma,
               given_ind = given_ind,
               given_x = given_x_mat)
par1$mean

# Duplicates and omitted indexes are allowed for dependent_ind
# For given_ind duplicates are not allowed
# Let's calculate conditional parameters for (X5, X2, X5 | X1 = -1, X3 = 1):
dependent_ind <- c(5, 2, 5)
par2 <- cmnorm(mean = mean, sigma = sigma,
               given_ind = given_ind,
               given_x = given_x,
               dependent_ind = dependent_ind)
  # E(X5, X2, X5 | X1 = -1, X3 = 1)
par2$mean
  # Cov(X5, X2, X5 | X1 = -1, X3 = 1)
par2$sigma

Density of (conditional) multivariate normal distribution

Description

This function calculates and differentiates density of (conditional) multivariate normal distribution.

Usage

dmnorm(
  x,
  mean,
  sigma,
  given_ind = numeric(),
  log = FALSE,
  grad_x = FALSE,
  grad_sigma = FALSE,
  is_validation = TRUE,
  control = NULL,
  n_cores = 1L
)

Arguments

Details

Consider notations from the Details section of cmnorm. The function calculates density f(x^{(d)}|x^{(g)}) of conditioned multivariate normal vector X_{I_{d}} | X_{I_{g}} = x^{(g)}. Where x^{(d)} is a subvector of x associated with X_{I_{d}} i.e. unconditioned components. Therefore x[given_ind] represents x^{(g)} while x[-given_ind] is x^{(d)}.

If grad_x is TRUE then function additionally estimates the gradient respect to both unconditioned and conditioned components:

\nabla f(x^{(d)}|x^{(g)})= \left(\frac{\partial f(x^{(d)}|x^{(g)})}{\partial x_{1}} ,..., \frac{\partial f(x^{(d)}|x^{(g)})}{\partial x_{m}}\right),

where each x_{i} belongs either to x^{(d)} or x^{(g)} depending on whether i\in I_{d} or i\in I_{g} correspondingly. In particular subgradients of density function respect to x^{(d)} and x^{(g)} are of the form:

\nabla_{x^{(d)}}\ln f(x^{(d)}|x^{(g)}) = -\left(x^{(d)}-\mu_{c}\right)\Sigma_{c}^{-1}

\nabla_{x^{(g)}}\ln f(x^{(d)}|x^{(g)}) = -\nabla_{x^{(d)}}f(x^{(d)}|x^{(g)})\Sigma_{d,g}\Sigma_{g,g}^{-1}

If grad_sigma is TRUE then function additionally estimates the gradient respect to the elements of covariance matrix \Sigma. For i\in I_{d} and j\in I_{d} the function calculates:

\frac{\partial \ln f(x^{(d)}|x^{(g)})}{\partial \Sigma_{i, j}} = \left(\frac{\partial \ln f(x^{(d)}|x^{(g)})}{\partial x_{i}} \times \frac{\partial \ln f(x^{(d)}|x^{(g)})}{\partial x_{j}} - \Sigma_{c,(i, j)}^{-1}\right) / \left(1 + I(i=j)\right),

where I(i=j) is an indicator function which equals 1 when the condition i=j is satisfied and 0 otherwise.

For i\in I_{d} and j\in I_{g} the following formula is used:

\frac{\partial \ln f(x^{(d)}|x^{(g)})}{\partial \Sigma_{i, j}} = -\frac{\partial \ln f(x^{(d)}|x^{(g)})}{\partial x_{i}} \times \left(\left(x^{(g)}-\mu_{g}\right)\Sigma_{g,g}^{-1}\right)_{q_{g}(j)}-

-\sum\limits_{k=1}^{n_{d}}(1+I(q_{d}(i)=k))\times (\Sigma_{d,g}\Sigma_{g,g}^{-1})_{k,q_{g}(j)}\times \frac{\partial \ln f(x^{(d)}|x^{(g)})}{\partial \Sigma_{i, q^{-1}_{d}(k)}},

where q_{g}(j)=\sum\limits_{k=1}^{m} I\left(I_{g,k} \leq j\right) and q_{d}(i)=\sum\limits_{k=1}^{m} I\left(I_{d,k} \leq i\right) represent the order of the i-th and j-th elements in I_{g} and I_{d} correspondingly i.e. x_{i}=x^{(d)}_{q_{d}(i)}=x_{I_{d, q_{d}(i)}} and x_{j}=x^{(g)}_{q_{g}(j)}=x_{I_{g, q_{g}(j)}}. Note that q_{g}(j)^{-1} and q_{d}(i)^{-1} are inverse functions. Number of conditioned and unconditioned components are denoted by n_{g}=\sum\limits_{k=1}^{m}I(k\in I_{g}) and n_{d}=\sum\limits_{k=1}^{m}I(k\in I_{d}) respectively. For the case i\in I_{g} and j\in I_{d} the formula is similar.

For i\in I_{g} and j\in I_{g} the following formula is used:

\frac{\partial \ln f(x^{(d)}|x^{(g)})}{\partial \Sigma_{i, j}} = -\nabla_{x^{(d)}}\ln f(x^{(d)}|x^{(g)})\times \left(x^{(g)}\times(\Sigma_{d,g} \times \Sigma_{g,g}^{-1} \times I_{g}^{*} \times \Sigma_{g,g}^{-1})^{T}\right)^T -

-\sum\limits_{k_{1}=1}^{n_{d}}\sum\limits_{k_{2}=k_{1}}^{n_{d}} \frac{\partial \ln f(x^{(d)}|x^{(g)})} {\partial \Sigma_{q_{d}(k_{1})^{-1}, q_{d}(k_{2})^{-1}}} \left(\Sigma_{d,g} \times \Sigma_{g,g}^{-1} \times I_{g}^{*} \times \Sigma_{g,g}^{-1}\times\Sigma_{d,g}^T\right)_{q_{d}(k_{1})^{-1}, q_{d}(k_{2})^{-1}},

where I_{g}^{*} is a square n_{g}-dimensional matrix of zeros except I_{g,(i,j)}^{*}=I_{g,(j,i)}^{*}=1.

Argument given_ind represents I_{g} and it should not contain any duplicates. The order of given_ind elements does not matter so it has no impact on the output.

More details on abovementioned differentiation formulas could be found in the appendix of E. Kossova and B. Potanin (2018).

Currently control has no input arguments intended for the users. This argument is used for some internal purposes of the package.

Value

This function returns an object of class "mnorm_dmnorm".

An object of class "mnorm_dmnorm" is a list containing the following components:

• den - density function value at x.

• grad_x - gradient of density respect to x if grad_x or grad_sigma input argument is set to TRUE.

• grad_sigma - gradient respect to the elements of sigma if grad_sigma input argument is set to TRUE.

If log is TRUE then den is a log-density so output grad_x and grad_sigma are calculated respect to the log-density.

Output grad_x is a Jacobian matrix which rows are gradients of the density function calculated for each row of x. Therefore grad_x[i, j] is a derivative of the density function respect to the j-th argument at point x[i, ].

Output grad_sigma is a 3D array such that grad_sigma[i, j, k] is a partial derivative of the density function respect to the sigma[i, j] estimated for the observation x[k, ].

References

E. Kossova., B. Potanin (2018). Heckman method and switching regression model multivariate generalization. Applied Econometrics, vol. 50, pages 114-143.

Examples

# Consider multivariate normal vector:
# X = (X1, X2, X3, X4, X5) ~ N(mean, sigma)

# Prepare multivariate normal vector parameters
  # expected value
mean <- c(-2, -1, 0, 1, 2)
n_dim <- length(mean)
  # correlation matrix
cor <- c(   1,  0.1,  0.2,   0.3,  0.4,
          0.1,    1, -0.1,  -0.2, -0.3,
          0.2, -0.1,    1,   0.3,  0.2,
          0.3, -0.2,  0.3,     1, -0.05,
          0.4, -0.3,  0.2, -0.05,     1)
cor <- matrix(cor, ncol = n_dim, nrow = n_dim, byrow = TRUE)
  # covariance matrix
sd_mat <- diag(c(1, 1.5, 2, 2.5, 3))
sigma <- sd_mat %*% cor %*% t(sd_mat)

# Estimate the density of X at point (-1, 0, 1, 2, 3)
x <- c(-1, 0, 1, 2, 3)
d.list <- dmnorm(x = x, mean = mean, sigma = sigma)
d <- d.list$den
print(d)

# Estimate the density of X at points
# x=(-1, 0, 1, 2, 3) and y=(-1.2, -1.5, 0, 1.2, 1.5)
y <- c(-1.5, -1.2, 0, 1.2, 1.5)
xy <- rbind(x, y)
d.list.1 <- dmnorm(x = xy, mean = mean, sigma = sigma)
d.1 <- d.list.1$den
print(d.1)

# Estimate the density of Xc=(X2, X4, X5 | X1 = -1, X3 = 1) at 
# point xd=(0, 2, 3) given conditioning values xg=(-1, 1)
given_ind <- c(1, 3)
d.list.2 <- dmnorm(x = x, mean = mean, sigma = sigma, 
                   given_ind = given_ind)
d.2 <- d.list.2$den
print(d.2)

# Estimate the gradient of density respect to the argument and 
# covariance matrix at points 'x' and 'y'
d.list.3 <- dmnorm(x = xy, mean = mean, sigma = sigma,
                   grad_x = TRUE, grad_sigma = TRUE)
# Gradient respect to the argument
grad_x.3 <- d.list.3$grad_x
   # at point 'x'
print(grad_x.3[1, ])
   # at point 'y'
print(grad_x.3[2, ])
# Partial derivative at point 'y' respect 
# to the 3-rd argument
print(grad_x.3[2, 3])
# Gradient respect to the covariance matrix
grad_sigma.3 <- d.list.3$grad_sigma
# Partial derivative respect to sigma(3, 5) at 
# point 'y'
print(grad_sigma.3[3, 5, 2])

# Estimate the gradient of the log-density function of 
# Xc=(X2, X4, X5 | X1 = -1, X3 = 1) and Yc=(X2, X4, X5 | X1 = -1.5, X3 = 0)
# respect to the argument and covariance matrix at 
# points xd=(0, 2, 3) and yd=(-1.2, 0, 1.5) respectively given
# conditioning values xg=(-1, 1) and yg=(-1.5, 0) correspondingly
d.list.4 <- dmnorm(x = xy, mean = mean, sigma = sigma,
                   grad_x = TRUE, grad_sigma = TRUE,
                   given_ind = given_ind, log = TRUE)
# Gradient respect to the argument
grad_x.4 <- d.list.4$grad_x
   # at point 'xd'
print(grad_x.4[1, ])
   # at point 'yd'
print(grad_x.4[2, ])
# Partial derivative at point 'xd' respect to 'xg[2]'
print(grad_x.4[1, 3])
# Partial derivative at point 'yd' respect to 'yd[5]'
print(grad_x.4[2, 5])
# Gradient respect to the covariance matrix
grad_sigma.4 <- d.list.4$grad_sigma
# Partial derivative respect to sigma(3, 5) at 
# point 'yd'
print(grad_sigma.4[3, 5, 2])

# Compare analytical gradients from the previous example with
# their numeric (forward difference) analogues at point 'xd'
# given conditioning 'xg'
delta <- 1e-6
grad_x.num <- rep(NA, 5)
grad_sigma.num <- matrix(NA, nrow = 5, ncol = 5)
for (i in 1:5)
{
  x.delta <- x
  x.delta[i] <- x[i] + delta
  d.list.delta <- dmnorm(x = x.delta, mean = mean, sigma = sigma,
                         grad_x = TRUE, grad_sigma = TRUE,
                         given_ind = given_ind, log = TRUE)
  grad_x.num[i] <- (d.list.delta$den - d.list.4$den[1]) / delta
   for(j in 1:5)
   {
     sigma.delta <- sigma
     sigma.delta[i, j] <- sigma[i, j] + delta 
     sigma.delta[j, i] <- sigma[j, i] + delta 
     d.list.delta <- dmnorm(x = x, mean = mean, sigma = sigma.delta,
                            grad_x = TRUE, grad_sigma = TRUE,
                            given_ind = given_ind, log = TRUE)
     grad_sigma.num[i, j] <- (d.list.delta$den - d.list.4$den[1]) / delta
   }
}
# Comparison of gradients respect to the argument
h.x <- cbind(analytical = grad_x.4[1, ], numeric = grad_x.num)
rownames(h.x) <- c("xg[1]", "xd[1]", "xg[2]", "xd[3]", "xd[4]")
print(h.x)
# Comparison of gradients respect to the covariance matrix
h.sigma <- list(analytical = grad_sigma.4[, , 1], numeric = grad_sigma.num)
print(h.sigma)

Convert base representation of a number into integer

Description

Converts base representation of a number into integer.

Usage

fromBase(x, base = 2L)

Arguments

Value

The function returns a positive integer that is a conversion from base under given coefficients x.

Examples

fromBase(c(1, 2, 0, 2, 3), 5)

Halton sequence

Description

Calculate elements of the Halton sequence and of some other pseudo-random sequences.

Usage

halton(
  n = 1L,
  base = as.integer(c(2)),
  start = 1L,
  random = "NO",
  type = "halton",
  scrambler = "NO",
  is_validation = TRUE,
  n_cores = 1L
)

Arguments

Details

Function seqPrimes could be used to provide the prime numbers for the base input argument.

Value

The function returns a matrix which i-th column is a sequence with base base[i] and elements with indexes from start to start + n.

References

J. Halton (1964) <doi:10.2307/2347972>

S. Kolenikov (2012) <doi:10.1177/1536867X1201200103>

Examples

halton(n = 100, base = c(2, 3, 5), start = 10)

Differentiate Regularized Incomplete Beta Function.

Description

Calculate derivatives of the regularized incomplete beta function that is a cumulative distribution function of the beta distribution.

Usage

pbetaDiff(x, p = 10, q = 0.5, n = 10L, is_validation = TRUE, control = NULL)

Arguments

Details

The function implements differentiation algorithm of R. Boik and J. Robinson-Cox (1998). Currently only first-order derivatives are considered.

Value

The function returns a list which has the following elements:

• dx - numeric vector of derivatives respect to each element of x.

• dp - numeric vector of derivatives respect to p for each element of x.

• dq - numeric vector of derivatives respect to q for each element of x.

References

Boik, R. J. and Robinson-Cox, J. F. (1998). Derivatives of the Incomplete Beta Function. Journal of Statistical Software, 3 (1), pages 1-20.

Examples

# Some values from Table 1 of R. Boik and J. Robinson-Cox (1998)
pbetaDiff(x = 0.001,  p = 1.5, q = 11)
pbetaDiff(x = 0.5,  p = 1.5, q = 11)

# Compare analytical and numeric derivatives
delta <- 1e-6
x <- c(0.01, 0.25, 0.5, 0.75, 0.99)
p <- 5
q <- 10
out <- pbetaDiff(x = x, p = p, q = q)
p0 <- pbeta(q = x, shape1 = p, shape2 = q)

# Derivatives respect to x
p1 <- pbeta(q = x + delta, shape1 = p, shape2 = q)
data.frame(numeric = (p1 - p0) / delta, analytical = out$dx)
  
# Derivatives respect to p
p1 <- pbeta(q = x, shape1 = p + delta, shape2 = q)
data.frame(numeric = (p1 - p0) / delta, analytical = out$dp)

# Derivatives respect to q
p1 <- pbeta(q = x, shape1 = p, shape2 = q + delta)
data.frame(numeric = (p1 - p0) / delta, analytical = out$dq)

Probabilities of (conditional) multivariate normal distribution

Description

This function calculates and differentiates probabilities of (conditional) multivariate normal distribution.

Usage

pmnorm(
  lower,
  upper,
  given_x = numeric(),
  mean = numeric(),
  sigma = matrix(),
  given_ind = numeric(),
  n_sim = 1000L,
  method = "default",
  ordering = "mean",
  log = FALSE,
  grad_lower = FALSE,
  grad_upper = FALSE,
  grad_sigma = FALSE,
  grad_given = FALSE,
  is_validation = TRUE,
  control = NULL,
  n_cores = 1L,
  marginal = NULL,
  grad_marginal = FALSE,
  grad_marginal_prob = FALSE
)

Arguments

Details

Consider notations from the Details sections of cmnorm and dmnorm. The function calculates probabilities of the form:

P\left(x^{(l)}\leq X_{I_{d}}\leq x^{(u)}|X_{I_{g}}=x^{(g)}\right)

where x^{(l)} and x^{(u)} are lower and upper integration limits respectively i.e. lower and upper correspondingly. Also x^{(g)} represents given_x. Note that lower and upper should be matrices of the same size. Also given_x should have the same number of rows as lower and upper.

To calculate bivariate probabilities the function applies the method of Drezner and Wesolowsky described in A. Genz (2004). In contrast to the classical implementation of this method the function applies Gauss-Legendre quadrature with 30 sample points to approximate integral (1) of A. Genz (2004). Classical implementations of this method use up to 20 points but requires some additional transformations of (1). During preliminary testing it has been found that approach with 30 points provides similar accuracy being slightly faster because of better vectorization capabilities.

To calculate trivariate probabilities the function uses Drezner method following formula (14) of A. Genz (2004). Similarly to bivariate case 30 points are used in Gauss-Legendre quadrature.

The function may apply the method of Gassmann (2003) for estimation of m>3 dimensional normal probabilities. It uses matrix 5 representation of Gassmann (2003) and 30 points of Gauss-Legendre quadrature.

For m-variate probabilities, where m > 1, the function may apply GHK algorithm described in section 4.2 of A. Genz and F. Bretz (2009). The implementation of GHK is based on deterministic Halton sequence with n_sim draws and uses variable reordering suggested in section 4.1.3 of A. Genz and F. Bretz (2009). The ordering algorithm may be determined via ordering argument. Available options are "NO", "mean" (default), and "variance".

Univariate probabilities are always calculated via standard approach so in this case method will not affect the output. If method = "Gassmann" then the function uses fast (aforementioned) algorithms for bivariate and trivariate probabilities or the method of Gassmann for m>3 dimensional probabilities. If method = "GHK" then GHK algorithm is used. If method = "default" then "Gassmann" is used for bivariate and trivariate probabilities while "GHK" is used for m>3 dimensional probabilities. During future updates "Gassmann" may become a default method for calculation of the 4-5 dimensional probabilities.

We are going to provide alternative estimation algorithms during future updates. These methods will be available via method argument.

The function is optimized to perform much faster when all upper integration limits upper are finite while all lower integration limits lower are negative infinite. The derivatives could be also calculated much faster when some integration limits are infinite.

For simplicity of notations further let's consider unconditioned probabilities. Derivatives respect to conditioned components are similar to those mentioned in Details section of dmnorm. We also provide formulas for m \geq 3. But the function may calculate derivatives for m \leq 2 using some simplifications of the formulas mentioned below.

If grad_upper is TRUE then function additionally estimates the gradient respect to upper:

\frac{\partial P\left(x^{(l)}\leq X\leq x^{(u)}\right)}{\partial x^{(u)}_{i}}= P\left(x^{(l)}_{(-i)}\leq X_{(-i)}\leq x^{(u)}_{(-i)}| X_{i}=x^{(u)}_{i}\right) f_{X_{i}}\left(x^{(u)}_{i};\mu_{i},\Sigma_{i,i}\right)

If grad_lower is TRUE then function additionally estimates the gradient respect to lower:

\frac{\partial P\left(x^{(l)}\leq X\leq x^{(u)}\right)}{\partial x^{(l)}_{i}}= -P\left(x^{(l)}_{(-i)}\leq X_{(-i)}\leq x^{(u)}_{(-i)}| X_{i}=x^{(l)}_{i}\right) f_{X_{i}}\left(x^{(l)}_{i};\mu_{i},\Sigma_{i,i}\right)

If grad_sigma is TRUE then function additionally estimates the gradient respect to sigma. For i\ne j the function calculates derivatives respect to the covariances:

\frac{\partial P\left(x^{(l)}\leq X\leq x^{(u)}\right)}{\partial \Sigma_{i, j}}=

=P\left(x^{(l)}_{(-(i, j))}\leq X_{-(i, j)}\leq x^{(u)}_{(-(i, j))}| X_{i}=x^{(u)}_{i}, X_{j}=x^{(u)}_{j}\right) f_{X_{i}, X_{j}}\left(x^{(u)}_{i}, x^{(u)}_{j}; \mu_{(i,j)},\Sigma_{(i, j),(i, j)}\right) -

-P\left(x^{(l)}_{(-(i, j))}\leq X_{-(i, j)}\leq x^{(u)}_{(-(i, j))}| X_{i}=x^{(l)}_{i}, X_{j}=x^{(u)}_{j}\right) f_{X_{i}, X_{j}}\left(x^{(l)}_{i}, x^{(u)}_{j}; \mu_{(i,j)},\Sigma_{(i, j),(i, j)}\right) -

-P\left(x^{(l)}_{(-(i, j))}\leq X_{-(i, j)}\leq x^{(u)}_{(-(i, j))}| X_{i}=x^{(u)}_{i}, X_{j}=x^{(l)}_{j}\right) f_{X_{i}, X_{j}}\left(x^{(u)}_{i}, x^{(l)}_{j}; \mu_{(i,j)},\Sigma_{(i, j),(i, j)}\right) +

+P\left(x^{(l)}_{(-(i, j))}\leq X_{-(i, j)}\leq x^{(u)}_{(-(i, j))}| X_{i}=x^{(l)}_{i}, X_{j}=x^{(l)}_{j}\right) f_{X_{i}, X_{j}}\left(x^{(l)}_{i}, x^{(l)}_{j}; \mu_{(i,j)},\Sigma_{(i, j),(i, j)}\right)

Note that if some of integration limits are infinite then some elements of this equation converge to zero which highly simplifies the calculations.

Derivatives respect to variances are calculated using derivatives respect to covariances and integration limits:

\frac{\partial P\left(x^{(l)}\leq X\leq x^{(u)}\right)}{\partial \Sigma_{i, i}}=

-\frac{\partial P\left(x^{(l)}\leq X\leq x^{(u)}\right)}{\partial x^{(l)}_{i}} \frac{x^{(l)}_{i}}{2\Sigma_{i, i}} -\frac{\partial P\left(x^{(l)}\leq X\leq x^{(u)}\right)}{\partial x^{(u)}_{i}} \frac{x^{(u)}_{i}}{2\Sigma_{i, i}}-

-\sum\limits_{j\ne i}\frac{\partial P\left(x^{(l)}\leq X\leq x^{(u)}\right)}{\partial \Sigma_{i, j}} \frac{\Sigma_{i, j}}{2\Sigma_{i, i}}

If grad_given is TRUE then function additionally estimates the gradient respect to given_x using formulas similar to those described in Details section of dmnorm.

More details on abovementioned differentiation formulas could be found in the appendix of E. Kossova and B. Potanin (2018).

If marginal is not empty then Gaussian copula will be used instead of a classical multivariate normal distribution. Without loss of generality let's assume that \mu is a vector of zeros and introduce some additional notations:

q_{i}^{(u)} = \Phi^{-1}\left(P_{i}\left(\frac{x_{i}^{(u)}}{\sigma_{i}}\right)\right)

q_{i}^{(l)} = \Phi^{-1}\left(P_{i}\left(\frac{x_{i}^{(l)}}{\sigma_{i}}\right)\right)

where \Phi(.)^{-1} is a quantile function of a standard normal distribution and P_{i} is a cumulative distribution function of the standartized (to zero mean and unit variance) marginal distribution which name and parameters are defined by names(marginal)[i] and marginal[[i]] correspondingly. For example if marginal[i] = "logistic" then:

P_{i}(t) = \frac{1}{1+e^{-\pi t / \sqrt{3}}}

Let's denote by X^{*} random vector which is distributed with Gaussian (its covariance matrix is \Sigma) copula and marginals defined by marginal. Then probabilities for these random vector are calculated as follows:

P\left(x^{(l)}\leq X^{*}\leq x^{(u)}\right) = P\left(\sigma q^{(l)}\leq X\leq \sigma q^{(u)}\right) = P_{0}\left(\sigma q^{(l)}, \sigma q^{(u)}\right)

where q^{(l)} = (q_{1}^{(l)},...,q_{m}^{(l)}), q^{(u)} = (q_{1}^{(u)},...,q_{m}^{(u)}) and \sigma = (\sqrt{\Sigma_{1, 1}},...,\sqrt{\Sigma_{m, m}}). Therefore probabilities of X^{*} may be calculated using probabilities of multivariate normal random vector X (with the same covariance matrix) by substituting lower and upper integration limits x^{(l)} and x^{(u)} with \sigma q^{(l)} and \sigma q^{(u)} correspondingly. Therefore differentiation formulas are similar to those mentioned above and will be automatically adjusted if marginal is not empty (including conditional probabilities).

Argument control is a list with the following input parameters:

• random_sequence – numeric matrix of uniform pseudo random numbers (like Halton sequence). The number of columns should be equal to the number of dimensions of multivariate random vector. If omitted than this matrix will be generated automatically using n_sim number of simulations.

Value

This function returns an object of class "mnorm_pmnorm".

An object of class "mnorm_pmnorm" is a list containing the following components:

• prob - probability that multivariate normal random variable will be between lower and upper bounds.

• grad_lower - gradient of probability respect to lower if grad_lower or grad_sigma input argument is set to TRUE.

• grad_upper - gradient of probability respect to upper if grad_upper or grad_sigma input argument is set to TRUE.

• grad_sigma - gradient respect to the elements of sigma if grad_sigma input argument is set to TRUE.

• grad_given - gradient respect to the elements of given_x if grad_given input argument is set to TRUE.

• grad_marginal - gradient respect to the elements of marginal_par if grad_marginal input argument is set to TRUE. Currently only derivatives respect to the parameters of "PGN" distribution are available.

If log is TRUE then prob is a log-probability so output grad_lower, grad_upper, grad_sigma and grad_given are calculated respect to the log-probability.

Output grad_lower and grad_upper are Jacobian matrices which rows are gradients of the probabilities calculated for each row of lower and upper correspondingly. Similarly grad_given is a Jacobian matrix respect to given_x.

Output grad_sigma is a 3D array such that grad_sigma[i, j, k] is a partial derivative of the probability function respect to the sigma[i, j] estimated for k-th observation.

Output grad_marginal is a list such that grad_marginal[[i]] is a Jacobian matrice which rows are gradients of the probabilities calculated for each row of lower and upper correspondingly respect to the elements of marginal_par[[i]].

References

Genz, A. (2004), Numerical computation of rectangular bivariate and trivariate normal and t-probabilities, Statistics and Computing, 14, 251-260.

Genz, A. and Bretz, F. (2009), Computation of Multivariate Normal and t Probabilities. Lecture Notes in Statistics, Vol. 195. Springer-Verlag, Heidelberg.

E. Kossova, B. Potanin (2018). Heckman method and switching regression model multivariate generalization. Applied Econometrics, vol. 50, pages 114-143.

H. I. Gassmann (2003). Multivariate Normal Probabilities: Implementing an Old Idea of Plackett's. Journal of Computational and Graphical Statistics, vol. 12 (3), pages 731-752.

Examples

# Consider multivariate normal vector:
# X = (X1, X2, X3, X4, X5) ~ N(mean, sigma)

# Prepare multivariate normal vector parameters
  # expected value
mean <- c(-2, -1, 0, 1, 2)
n_dim <- length(mean)
  # correlation matrix
cor <- c(   1,  0.1,  0.2,   0.3,  0.4,
          0.1,    1, -0.1,  -0.2, -0.3,
          0.2, -0.1,    1,   0.3,  0.2,
          0.3, -0.2,  0.3,     1, -0.05,
          0.4, -0.3,  0.2, -0.05,     1)
cor <- matrix(cor, ncol = n_dim, nrow = n_dim, byrow = TRUE)
  # covariance matrix
sd_mat <- diag(c(1, 1.5, 2, 2.5, 3))
sigma <- sd_mat %*% cor %*% t(sd_mat)

# Estimate probability:
# P(-3 < X1 < 1, -2.5 < X2 < 1.5, -2 < X3 < 2, -1.5 < X4 < 2.5, -1 < X5 < 3)
lower <- c(-3, -2.5, -2, -1.5, -1)
upper <- c(1, 1.5, 2, 2.5, 3)
p.list <- pmnorm(lower = lower, upper = upper,
                 mean = mean, sigma = sigma)
p <- p.list$prob
print(p)

# Additionally estimate a probability
lower.1 <- c(-Inf, 0, -Inf, 1, -Inf)
upper.1 <- c(Inf, Inf, 3, 4, 5)
lower.mat <- rbind(lower, lower.1)
upper.mat <- rbind(upper, upper.1)
p.list.1 <- pmnorm(lower = lower.mat, upper = upper.mat,
                   mean = mean, sigma = sigma)
p.1 <- p.list.1$prob
print(p.1)

# Estimate the probabilities
# P(-1 < X1 < 1, -3 < X3 < 3, -5 < X5 < 5 | X2 = -2, X4 = 4)
lower.2 <- c(-1, -3, -5)
upper.2 <- c(1, 3, 5)
given_ind <- c(2, 4)
given_x <- c(-2, 4)
p.list.2 <- pmnorm(lower = lower.2, upper = upper.2,
                   mean = mean, sigma = sigma,
                   given_ind = given_ind, given_x = given_x)
p.2 <- p.list.2$prob
print(p.2)

# Additionally estimate the probability
# P(-Inf < X1 < 1, -3 < X3 < Inf, -Inf < X5 < Inf | X2 = 4, X4 = -2)
lower.3 <- c(-Inf, -3, -Inf)
upper.3 <- c(1, Inf, Inf)
given_x.1 <- c(-2, 4)
lower.mat.2 <- rbind(lower.2, lower.3)
upper.mat.2 <- rbind(upper.2, upper.3)
given_x.mat <- rbind(given_x, given_x.1)
p.list.3 <- pmnorm(lower = lower.mat.2, upper = upper.mat.2,
                   mean = mean, sigma = sigma,
                   given_ind = given_ind, given_x = given_x.mat)
p.3 <- p.list.3$prob
print(p.3)

# Estimate the gradient of previous probabilities respect various arguments
p.list.4 <- pmnorm(lower = lower.mat.2, upper = upper.mat.2,
                   mean = mean, sigma = sigma,
                   given_ind = given_ind, given_x = given_x.mat,
                   grad_lower = TRUE, grad_upper = TRUE,
                   grad_sigma = TRUE, grad_given = TRUE)
p.4 <- p.list.4$prob
print(p.4)
# Gradient respect to 'lower'
grad_lower <- p.list.4$grad_lower
   # for the first probability
print(grad_lower[1, ])
   # for the second probability
print(grad_lower[2, ])
# Gradient respect to 'upper'
grad_upper <- p.list.4$grad_upper
   # for the first probability
print(grad_upper[1, ])
   # for the second probability
print(grad_upper[2, ])
# Gradient respect to 'given_x'
grad_given <- p.list.4$grad_given
   # for the first probability
print(grad_given[1, ])
   # for the second probability
print(grad_given[2, ])
# Gradient respect to 'sigma'
grad_given <- p.list.4$grad_given
   # for the first probability
print(grad_given[1, ])
   # for the second probability
print(grad_given[2, ])


# Compare analytical gradients from the previous example with
# their numeric (forward difference) analogues for the first probability
n_dependent <- length(lower.2)
n_given <- length(given_x)
n_dim <- n_dependent + n_given
delta <- 1e-6
grad_lower.num <- rep(NA, n_dependent)
grad_upper.num <- rep(NA, n_dependent)
grad_given.num <- rep(NA, n_given)
grad_sigma.num <- matrix(NA, nrow = n_dim, ncol = n_dim)
for (i in 1:n_dependent)
{
  # respect to lower
  lower.delta <- lower.2
  lower.delta[i] <- lower.2[i] + delta
  p.list.delta <- pmnorm(lower = lower.delta, upper = upper.2,
                         given_x = given_x,
                         mean = mean, sigma = sigma,
                         given_ind = given_ind)
  grad_lower.num[i] <- (p.list.delta$prob - p.list.4$prob[1]) / delta
  # respect to upper
  upper.delta <- upper.2
  upper.delta[i] <- upper.2[i] + delta
  p.list.delta <- pmnorm(lower = lower.2, upper = upper.delta,
                         given_x = given_x,
                         mean = mean, sigma = sigma,
                         given_ind = given_ind)
  grad_upper.num[i] <- (p.list.delta$prob - p.list.4$prob[1]) / delta
}
for (i in 1:n_given)
{
  # respect to lower
  given_x.delta <- given_x
  given_x.delta[i] <- given_x[i] + delta
  p.list.delta <- pmnorm(lower = lower.2, upper = upper.2,
                         given_x = given_x.delta,
                         mean = mean, sigma = sigma,
                         given_ind = given_ind)
  grad_given.num[i] <- (p.list.delta$prob - p.list.4$prob[1]) / delta
}
for (i in 1:n_dim)
{
  for(j in 1:n_dim)
  {
    # respect to sigma
    sigma.delta <- sigma
    sigma.delta[i, j] <- sigma[i, j] + delta 
    sigma.delta[j, i] <- sigma[j, i] + delta 
    p.list.delta <- pmnorm(lower = lower.2, upper = upper.2,
                           given_x = given_x,
                           mean = mean, sigma = sigma.delta,
                           given_ind = given_ind)
    grad_sigma.num[i, j] <- (p.list.delta$prob - p.list.4$prob[1]) / delta
  }
}
# Comparison of gradients respect to lower integration limits
h.lower <- cbind(analytical = p.list.4$grad_lower[1, ], 
                 numeric = grad_lower.num)
print(h.lower)
# Comparison of gradients respect to upper integration limits
h.upper <- cbind(analytical = p.list.4$grad_upper[1, ], 
                 numeric = grad_upper.num)
print(h.upper)
# Comparison of gradients respect to given values
h.given <- cbind(analytical = p.list.4$grad_given[1, ], 
                 numeric = grad_given.num)
print(h.given)
# Comparison of gradients respect to the covariance matrix
h.sigma <- list(analytical = p.list.4$grad_sigma[, , 1], 
                numeric = grad_sigma.num)
print(h.sigma)

# Let's again estimate probability
# P(-1 < X1 < 1, -3 < X3 < 3, -5 < X5 < 5 | X2 = -2, X4 = 4)
# But assume that standardized (to zero mean and unit variance): 
# 1) X1 and X2 have standardized PGN distribution with coefficients vectors
#    pc1 = (0.5, -0.2, 0.1) and pc2 = (0.2, 0.05) correspondingly.
# 2) X3 has standardized student distribution with 5 degrees of freedom
# 3) X4 has standardized logistic distribution
# 4) X5 has standard normal distribution
marginal <- list(PGN = c(0.5, -0.2, 0.1), hpa = c(0.2, 0.05), 
                 student = 5, logistic = numeric(), normal = NULL)
p.list.5 <- pmnorm(lower = lower.2, upper = upper.2,
                   mean = mean, sigma = sigma,
                   given_ind = given_ind, given_x = given_x,
                   grad_lower = TRUE, grad_upper = TRUE,
                   grad_sigma = TRUE, grad_given = TRUE,
                   marginal = marginal, grad_marginal = TRUE)     
# Lets investigate derivatives respect to parameters
# of marginal distributions
  # respect to pc1 of X1
p.list.5$grad_marginal[[1]]              
  # respect to pc2 of X2
p.list.5$grad_marginal[[2]]  
  # derivative respect to degrees of freedom of X5 is
  # currently unavailable and will be set to 0
p.list.5$grad_marginal[[3]]                    

Quantile function of a normal distribution

Description

Calculate quantile of a normal distribution using one of the available methods.

Usage

qnormFast(
  p,
  mean = 0L,
  sd = 1L,
  method = "Voutier",
  is_validation = TRUE,
  n_cores = 1L
)

Arguments

Details

If method = "Voutier" then the method of P. Voutier (2010) is used which maximum absolute error is about 0.000025. If method = "Shore" then the approach proposed by H. Shore (1982) is applied which maximum absolute error is about 0.026 for quantiles of level between 0.0001 and 0.9999.

Value

The function returns a vector of p-level quantiles of a normal distribution with mean equal to mean and standard deviation equal to sd.

References

H. Shore (1982) <doi:10.2307/2347972>

P. Voutier (2010) <doi:10.48550/arXiv.1002.0567>

Examples

qnormFast(c(0.1, 0.9), mean = 1, sd = 2)

Random number generator for (conditional) multivariate normal distribution

Description

This function generates random numbers (i.e. variates) from (conditional) multivariate normal distribution.

Usage

rmnorm(
  n,
  mean,
  sigma,
  given_ind = numeric(),
  given_x = numeric(),
  dependent_ind = numeric(),
  is_validation = TRUE,
  n_cores = 1L
)

Arguments

Details

This function uses Cholesky decomposition to generate multivariate normal variates from independent standard normal variates.

Value

This function returns a numeric matrix which rows a random variates from (conditional) multivariate normal distribution with mean equal to mean and covariance equal to sigma. If given_x and given_ind are also provided then random variates will be from conditional multivariate normal distribution. Please, see details section of cmnorm to get additional information on the conditioning procedure.

Examples

# Consider multivariate normal vector:
# X = (X1, X2, X3, X4, X5) ~ N(mean, sigma)

# Prepare multivariate normal vector parameters
  # expected value
mean <- c(-2, -1, 0, 1, 2)
n_dim <- length(mean)
  # correlation matrix
cor <- c(   1,  0.1,  0.2,   0.3,  0.4,
          0.1,    1, -0.1,  -0.2, -0.3,
          0.2, -0.1,    1,   0.3,  0.2,
          0.3, -0.2,  0.3,     1, -0.05,
          0.4, -0.3,  0.2, -0.05,     1)
cor <- matrix(cor, ncol = n_dim, nrow = n_dim, byrow = TRUE)
  # covariance matrix
sd_mat <- diag(c(1, 1.5, 2, 2.5, 3))
sigma <- sd_mat %*% cor %*% t(sd_mat)

# Simulate random variates from this distribution
rmnorm(n = 3, mean = mean, sigma = sigma)

# Simulate random variate from (X1, X3, X5 | X1 = -1, X4 = 1)
given_x <- c(-1, 1)
given_ind = c(1, 4)
rmnorm(n = 1, mean = mean, sigma = sigma, 
       given_x = given_x, given_ind = given_ind)

# Simulate random variate from (X1, X3, X5 | X1 = -1, X4 = 1)
# and (X1, X3, X5 | X1 = 2, X4 = 3)
given_x = rbind(c(-1, 1), c(2, 3))
rmnorm(n = nrow(given_x), mean = mean, sigma = sigma, 
       given_x = given_x, given_ind = given_ind)

Sequence of prime numbers

Description

Calculates the sequence of prime numbers.

Usage

seqPrimes(n)

Arguments

Value

The function returns a numeric vector containing first n prime numbers. The current (naive) implementation of the algorithm is not efficient in terms of speed so it is suited for low n < 10000 but requires just O(n) memory usage.

Examples

seqPrimes(10)

Standardized Student t Distribution

Description

These functions calculate and differentiate a cumulative distribution function and density function of the standardized (to zero mean and unit variance) Student distribution. Quantile function and random numbers generator are also provided.

Usage

dt0(x, df = 10, log = FALSE, grad_x = FALSE, grad_df = FALSE)

pt0(x, df = 10, log = FALSE, grad_x = FALSE, grad_df = FALSE, n = 10L)

rt0(n = 1L, df = 10)

qt0(x = 1L, df = 10)

Arguments

Details

Standardized (to zero mean and unit variance) Student distribution has the following density and cumulative distribution functions:

f(x) = \frac{\Gamma\left(\frac{v + 1}{2}\right)}{ \sqrt{(v - 2)\pi}\Gamma\left(\frac{v}{2}\right)} \left(1 + \frac{x^2}{v - 2}\right)^{-\frac{v+1}{2}},

F(x) = \begin{cases} 1 - \frac{1}{2}I(\frac{v - 2}{x^2 + v - 2}, \frac{v}{2}, \frac{1}{2})\text{, if }x \geq 0\\ \frac{1}{2}I(\frac{v - 2}{x^2 + v - 2}, \frac{v}{2}, \frac{1}{2})\text{, if }x < 0 \end{cases},

where v > 2 is the number of degrees of freedom df and I(.) is a cumulative distribution function of beta distribution which is calculated by pbeta function.

Value

Function rt0 returns a numeric vector of random numbers. Function qt0 returns a numeric vector of quantiles. Functions pt0 and dt0 return a list which may contain the following elements:

• prob - numeric vector of probabilities calculated for each element of x. Exclusively for pt0 function.

• den - numeric vector of densities calculated for each each element of x. Exclusively for dt0 function.

• grad_x - numeric vector of derivatives respect to p for each element of x. This element appears only if input argument grad_x is TRUE.

• grad_df - numeric vector of derivatives respect to q for each element of x. This element appears only if input argument grad_df is TRUE.

Examples

# Simple examples
pt0(x = 0.3, df = 10, log = FALSE, grad_x = TRUE, grad_df = TRUE)
dt0(x = 0.3, df = 10, log = FALSE, grad_x = TRUE, grad_df = TRUE)
qt0(x = 0.3, df = 10)

# Compare analytical and numeric derivatives
delta <- 1e-6
x <- c(-2, -1, 0, 1, 2)
df <- 5

# For probabilities
out <- pt0(x, df = df, grad_x = TRUE, grad_df = TRUE)
p0 <- out$prob
  # grad_x
p1 <- pt0(x + delta, df = df)$prob
data.frame(numeric = (p1 - p0) / delta, analytical = out$grad_x)
  # grad_df
p1 <- pt0(x, df = df + delta)$prob
data.frame(numeric = (p1 - p0) / delta, analytical = out$grad_df)

# For densities
out <- dt0(x, df = df, grad_x = TRUE, grad_df = TRUE)
p0 <- out$den
  # grad_x
p1 <- dt0(x + delta, df = df)$den
data.frame(numeric = (p1 - p0) / delta, analytical = out$grad_x)
  # grad_df
p1 <- dt0(x, df = df + delta)$den
data.frame(numeric = (p1 - p0) / delta, analytical = out$grad_df)

Convert integer value to other base

Description

Converts integer value to other base.

Usage

toBase(x, base = 2L)

Arguments

Value

The function returns a numeric vector containing representation of x in a base given in base.

Examples

toBase(888, 5)