hpaDist: Probabilities and Moments Hermite Polynomial Approximation

Description

Approximation of truncated, marginal and conditional densities, moments and cumulative probabilities of multivariate distributions via Hermite polynomial based approach proposed by Gallant and Nychka in 1987.

Density approximating function is scale adjusted product of two terms. The first one is squared multivariate polynomial of pol_degrees degrees with pol_coefficients coefficients vector. The second is product of independent normal random variables' densities with expected values and standard deviations given by mean and sd vectors correspondingly. Approximating function satisfies properties of density function thus generating a broad family of distributions. Characteristics of these distributions (moments, quantiles, probabilities and so on) may provide accurate approximations to characteristic of other distributions. Moreover it is usually possible to provide arbitrary close approximation by the means of polynomial degrees increase.

Usage

dhpa(
  x,
  pol_coefficients = numeric(0),
  pol_degrees = numeric(0),
  given_ind = logical(0),
  omit_ind = logical(0),
  mean = numeric(0),
  sd = numeric(0),
  is_parallel = FALSE,
  log = FALSE,
  is_validation = TRUE
)
phpa(
  x,
  pol_coefficients = numeric(0),
  pol_degrees = numeric(0),
  given_ind = logical(0),
  omit_ind = logical(0),
  mean = numeric(0),
  sd = numeric(0),
  is_parallel = FALSE,
  log = FALSE,
  is_validation = TRUE
)
ihpa(
  x_lower = matrix(1, 1),
  x_upper = matrix(1, 1),
  pol_coefficients = numeric(0),
  pol_degrees = numeric(0),
  given_ind = logical(0),
  omit_ind = logical(0),
  mean = numeric(0),
  sd = numeric(0),
  is_parallel = FALSE,
  log = FALSE,
  is_validation = TRUE
)
ehpa(
  x = matrix(1, 1),
  pol_coefficients = numeric(0),
  pol_degrees = numeric(0),
  given_ind = logical(0),
  omit_ind = logical(0),
  mean = numeric(0),
  sd = numeric(0),
  expectation_powers = numeric(0),
  is_parallel = FALSE,
  is_validation = TRUE
)
etrhpa(
  tr_left = matrix(1, 1),
  tr_right = matrix(1, 1),
  pol_coefficients = numeric(0),
  pol_degrees = numeric(0),
  mean = numeric(0),
  sd = numeric(0),
  expectation_powers = numeric(0),
  is_parallel = FALSE,
  is_validation = TRUE
)
dtrhpa(
  x,
  tr_left = matrix(),
  tr_right = matrix(),
  pol_coefficients = numeric(0),
  pol_degrees = numeric(0),
  given_ind = logical(0),
  omit_ind = logical(0),
  mean = numeric(0),
  sd = numeric(0),
  is_parallel = FALSE,
  log = FALSE,
  is_validation = TRUE
)
itrhpa(
  x_lower = matrix(1, 1),
  x_upper = matrix(1, 1),
  tr_left = matrix(1, 1),
  tr_right = matrix(1, 1),
  pol_coefficients = numeric(0),
  pol_degrees = numeric(0),
  given_ind = logical(0),
  omit_ind = logical(0),
  mean = numeric(0),
  sd = numeric(0),
  is_parallel = FALSE,
  log = FALSE,
  is_validation = TRUE
)
dhpaDiff(
  x,
  pol_coefficients = numeric(0),
  pol_degrees = numeric(0),
  given_ind = logical(0),
  omit_ind = logical(0),
  mean = numeric(0),
  sd = numeric(0),
  type = "pol_coefficients",
  is_parallel = FALSE,
  log = FALSE,
  is_validation = TRUE
)
ihpaDiff(
  x_lower = matrix(1, 1),
  x_upper = matrix(1, 1),
  pol_coefficients = numeric(0),
  pol_degrees = numeric(0),
  given_ind = logical(0),
  omit_ind = logical(0),
  mean = numeric(0),
  sd = numeric(0),
  type = "pol_coefficients",
  is_parallel = FALSE,
  log = FALSE,
  is_validation = TRUE
)

Arguments

numeric matrix of density function arguments and conditional values. Note that x rows are points (observations) while random vectors components (variables) are columns.

pol_coefficients

numeric vector of polynomial coefficients.

pol_degrees

non-negative integer vector of polynomial degrees (orders).

given_ind

logical vector indicating whether corresponding random vector component is conditioned. By default it is a logical vector of FALSE values. If give_ind[i] is TRUE then i-th column of x matrix will contain conditional values.

omit_ind

logical vector indicating whether corresponding component is omitted. By default it is a logical vector of FALSE values. If omit_ind[i] is TRUE then values in i-th column of x matrix will be ignored.

mean

numeric vector of expected values.

positive numeric vector of standard deviations.

is_parallel

if TRUE then multiple cores will be used for some calculations. It usually provides speed advantage for large enough samples (about more than 1000 observations).

log

logical; if TRUE then probabilities p are given as log(p) or derivatives will be given respect to log(p)

is_validation

logical value indicating whether function input arguments should be validated. Set it to FALSE for slight performance boost (default value is TRUE).

x_lower

numeric matrix of lower integration limits. Note that x_lower rows are observations while variables are columns.

x_upper

numeric matrix of upper integration limits. Note that x_upper rows are observations while variables are columns.

expectation_powers

integer vector of random vector components powers.

tr_left

numeric matrix of left (lower) truncation limits. Note that tr_left rows are observations while variables are columns. If tr_left and tr_right are single row matrices then the same truncation limits will be applied to all observations that are determined by the first rows of these matrices.

tr_right

numeric matrix of right (upper) truncation limits. Note that tr_right rows are observations while variables are columns. If tr_left and tr_right are single row matrices then the same truncation limits will be applied to all observations that are determined by the first rows of these matrices.

type

determines the partial derivatives to be included into the gradient. If type="pol_coefficients" then gradient will contain partial derivatives respect to polynomial coefficients listed in the same order as pol_coefficients. Other available types are type = "mean" and type = "sd". For function dhpaDiff it is possible to take gradient respect to the x points setting type="x". For function ihpaDiff it is possible to take gradient respect to the x lower and upper points setting type = "x_lower" or type = "upper" correspondingly. In order to get full gradient please set type="all".

Value

Functions dhpa, phpa and dtrhpa return vector of probabilities of length nrow(x).

Functions ihpa and itrhpa return vector of probabilities of length nrow(x_upper).

If x argument has not been provided or is a single row matrix then function ehpa returns moment value. Otherwise it returns vector of length nrow(x) containing moments values.

If tr_left and tr_right arguments are single row matrices then function etrhpa returns moment value. Otherwise it returns vector of length max(nrow(tr_left), nrow(tr_right)) containing moments values.

Functions dhpaDiff and ihpaDiff return Jacobin matrix. The number of columns depends on type argument. The number of rows is nrow(x) for dhpaDiff and nrow(x_upper) for ihpaDiff

If mean or sd are not specified they assume the default values of $m$ -dimensional vectors of 0 and 1, respectively. If x_lower is not specified then it is the matrix of the same size as x_upper containing negative infinity values only. If expectation_powers is not specified then it is $m$ -dimensional vector of 0 values.

Please see 'Details' section for additional information.

Details

It is possible to approximate densities dhpa, cumulative probabilities phpa, ihpa, moments ehpa as well as their truncated dtrhpa, itrhpa, etrhpa forms and gradients dhpaDiff, ihpaDiff. Note that phpa is special of ihpa where x corresponds to x_upper while x_lower is matrix of negative infinity values. So phpa intended to approximate cumulative distribution functions while ihpa approximates probabilities that random vector components will be between values determined by rows of x_lower and x_upper matrices. Further details are given below.

Since density approximating function is non-negative and integrates to 1 it is density function for some $m$ -variate random vector $ξ$ . Approximating function $f_{ξ} (x)$ has the following form: $f_{ξ} (x) = f_{ξ} (x; μ, σ, α) = \frac{1}{ψ} \prod_{t = 1}^{m} ϕ (x_{t}; μ_{t}, σ_{t}) {(\sum_{i_{1} = 0}^{K_{1}} . . . \sum_{i_{m} = 0}^{K_{m}} α_{(i_{1}, . . ., i_{m})} \prod_{r = 1}^{m} x_{r}^{i_{r}})}^{2}$ $ψ = \sum_{i_{1} = 0}^{K_{1}} . . . \sum_{i_{m} = 0}^{K_{m}} \sum_{j_{1} = 0}^{K_{1}} . . . \sum_{j_{m} = 0}^{K_{m}} α_{(i_{1}, \dots, i_{m})} α_{(j_{1}, \dots, j_{m})} \prod_{r = 1}^{m} M (i_{r} + j_{r}; μ_{r}, σ_{r}),$ where:

$x = (x_{1}, . . . x_{m})$ - is vector of arguments i.e. rows of x matrix in dhpa.

$α_{(i_{1}, \dots, i_{m})}$ - is polynomial coefficient corresponding to pol_coefficients[k] element. In order to investigate correspondence between k and $(i_{1}, \dots, i_{m})$ values please see 'Examples' section below or polynomialIndex function 'Details', 'Value' and 'Examples' sections. Note that if $m = 1$ then pol_coefficients[k] simply corresponds to $α_{k - 1}$ .

$(K_{1}, . . ., K_{m})$ - are polynomial degrees (orders) provided via pol_degrees argument so pol_degrees[i] determines $K_{i}$ .

$ϕ (.; μ_{t}, σ_{t})$ - is normal random variable density function where $μ_{t}$ and $σ_{t}$ are mean and standard deviation determined by mean[t] and sd[t] arguments values.

$M (q; μ_{r}, σ_{r})$ - is $q$ -th order moment of normal random variable with mean $μ_{r}$ and standard deviation $σ_{r}$ . Note that function normalMoment allows to calculate and differentiate normal random variable's moments.

$ψ$ - constant term insuring that $f_{ξ} (x)$ is density function.

Therefore dhpa allows to calculate $f_{ξ} (x)$ values at points determined by rows of x matrix given polynomial degrees pol_degrees ( $K$ ) as well as mean ( $μ$ ), sd ( $σ$ ) and pol_coefficients ( $α$ ) parameters values. Note that mean, sd and pol_degrees are $m$ -variate vectors while pol_coefficients has prod(pol_degrees + 1) elements.

Cumulative probabilities could be approximated as follows: $P ({\underset{―}{x}}_{1} \leq ξ_{1} \leq {\overset{―}{x}}_{1}, . . ., {\underset{―}{x}}_{m} \leq ξ_{m} \leq {\overset{―}{x}}_{m}) =$ $= {\bar{F}}_{ξ} (\underset{―}{x}, \bar{x}) = {\bar{F}}_{ξ} (\underset{―}{x}, \bar{x}; μ, σ, α) = \frac{1}{ψ} \prod_{t = 1}^{m} (Φ ({\bar{x}}_{t}; μ_{t}, σ_{t}) - Φ ({\underset{―}{x}}_{t}; μ_{t}, σ_{t})) *$ $* \sum_{i_{1} = 0}^{K_{1}} . . . \sum_{i_{m} = 0}^{K_{m}} \sum_{j_{1} = 0}^{K_{1}} . . . \sum_{j_{m} = 0}^{K_{m}} α_{(i_{1}, . . ., i_{m})} α_{(j_{1}, . . ., j_{m})} \prod_{r = 1}^{m} M_{T R} (i_{r} + j_{r}; {\underset{―}{x}}_{r}, {\overset{―}{x}}_{r}, μ_{r}, σ_{r})$

where:

$Φ (.; μ_{t}, σ_{t})$ - is normal random variable's cumulative distribution function where $μ_{t}$ and $σ_{t}$ are mean and standard deviation determined by mean[t] and sd[t] arguments values.

$M_{T R} (q; {\underset{―}{x}}_{r}, {\overset{―}{x}}_{r}, μ_{r}, σ_{r})$ - is $q$ -th order moment of truncated (from above by ${\overset{―}{x}}_{r}$ and from below by ${\underset{―}{x}}_{r}$ ) normal random variable with mean $μ_{r}$ and standard deviation $σ_{r}$ . Note that function truncatedNormalMoment allows to calculate and differentiate truncated normal random variable's moments.

$\overset{―}{x} = ({\overset{―}{x}}_{1}, . . ., {\overset{―}{x}}_{m})$ - vector of upper integration limits i.e. rows of x_upper matrix in ihpa.

$\underset{―}{x} = ({\underset{―}{x}}_{1}, . . ., {\underset{―}{x}}_{m})$ - vector of lower integration limits i.e. rows of x_lower matrix in ihpa.

Therefore ihpa allows to calculate interval distribution function ${\bar{F}}_{ξ} (\underset{―}{x}, \bar{x})$ values at points determined by rows of x_lower ( $\underset{―}{x}$ ) and x_upper ( $\overset{―}{x}$ ) matrices. The rest of the arguments are similar to dhpa.

Expected value powered product approximation is as follows:

$E (\prod_{t = 1}^{m} ξ_{t}^{k_{t}}) = \frac{1}{ψ} \sum_{i_{1} = 0}^{K_{1}} . . . \sum_{i_{m} = 0}^{K_{m}} \sum_{j_{1} = 0}^{K_{1}} . . . \sum_{j_{m} = 0}^{K_{m}} α_{(i_{1}, . . ., i_{m})} α_{(j_{1}, . . ., j_{m})} \prod_{r = 1}^{m} M (i_{r} + j_{r} + k_{t}; μ_{r}, σ_{r})$

where $(k_{1}, . . ., k_{m})$ are integer powers determined by expectation_powers argument of ehpa so expectation_powers[t] assigns $k_{t}$ . Note that argument x in ehpa allows to determined conditional values.

Expanding polynomial degrees $(K_{1}, . . ., K_{m})$ it is possible to provide arbitrary close approximation to density of some $m$ -variate random vector $ξ^{⋆}$ . So actually $f_{ξ} (x)$ approximates $f_{ξ^{⋆}} (x)$ . Accurate approximation requires appropriate mean, sd and pol_coefficients values selection. In order to get sample estimates of these parameters please apply hpaML function.

In order to perform calculation for marginal distribution of some $ξ$ components please provide omitted components via omit_ind argument. For examples if ones assume $m = 5$ -variate distribution and wants to deal with $1$ -st, $3$ -rd, and $5$ -th components only i.e. $(ξ_{1}, ξ_{3}, ξ_{5})$ then set omit_ind = c(FALSE, TRUE, FALSE, TRUE, FALSE) indicating that $ξ_{2}$ and $ξ_{4}$ should be 'omitted' from $ξ$ since $2$ -nd and $4$ -th values of omit_ind are TRUE. Then x still should be $5$ column matrix but values in $2$ -nd and $4$ -th columns will not affect calculation results. Meanwhile note that marginal distribution of t components of $ξ$ usually do not coincide with any marginal distribution generated by t-variate density approximating function.

In order to perform calculation for conditional distribution i.e. given fixed values for some $ξ$ components please provide these components via given_ind argument. For example if ones assume $m = 5$ -variate distribution and wants to deal with $1$ -st, $3$ -rd, and $5$ -th components given fixed values (suppose 8 and 10) for the other two components i.e. $(ξ | ξ_{2} = 8, ξ_{4} = 10)$ then set given_ind = c(FALSE, TRUE, FALSE, TRUE, FALSE) and x[2] = 8, x[4] = 10 where for simplicity it is assumed that x is single row $5$ column matrix; it is possible to provide different conditional values for the same components simply setting different values to different x rows.

Note that it is possible to combine given_ind and omit_ind arguments. However it is wrong to set both given_ind[i] and omit_ind[i] to TRUE. Also at least one value should be FALSE both for given_ind and omit_ind.

In order to consider truncated distribution of $ξ$ i.e. $(ξ | {\overset{―}{a}}_{1} \leq ξ_{1} \leq {\overset{―}{b}}_{1}, \dots, {\overset{―}{a}}_{m} \leq ξ_{m} \leq {\overset{―}{b}}_{m})$ please set lower (left) truncation points $\overset{―}{a}$ and upper (right) truncation points $\overset{―}{b}$ via tr_left and tr_right arguments correspondingly. Note that if lower truncation points are negative infinite and upper truncation points are positive infinite then dtrhpa, itrhpa and etrhpa are similar to dhpa, ihpa and ehpa correspondingly.

In order to calculate Jacobian of $f_{ξ} (x; μ, σ, α)$ and ${\bar{F}}_{ξ} (\underset{―}{x}, \bar{x}; μ, σ, α)$ w.r.t all ore some particular parameters please apply dhpaDiff and ihpaDiff functions correspondingly specifying parameters of interest via type argument. If x or x_lower and x_upper are single row matrices then gradients will be calculated.

For further information please see 'Examples' section. Note that examples are given separately for each function.

Densities Hermite polynomial approximation approach has been proposed by A. Gallant and D. W. Nychka in 1987. The main idea is to approximate unknown distribution density with scaled Hermite polynomial. For more information please refer to the literature listed below.

References

A. Gallant and D. W. Nychka (1987) <doi:10.2307/1913241>

Examples

Run this code

# NOT RUN {
## Example demonstrating dhpa function application.
## Let's approximate some three random variables (i.e. X1, X2 and X3) 
## joint density function at points  x = (0,1, 0.2, 0.3) and 
## y = (0.5, 0.8, 0.6) with Hermite polynomial of (1, 2, 3) degrees which 
## polynomial coefficients equal 1 except coefficient related  to x1*(x^3) 
## polynomial element which equals 2. Also suppose that normal density 
## related mean vector equals (1.1, 1.2, 1.3) while standard deviations 
## vector is (2.1, 2.2, 2.3).

# Prepare initial values
x <- matrix(c(0.1, 0.2, 0.3), nrow = 1)   # x point as a single row matrix
y <- matrix(c(0.5, 0.8, 0.6), nrow = 1)   # y point as a single row matrix
x_y <- rbind(x, y)                        # matrix which rows are x and y
mean <- c(1.1, 1.2, 1.3)
sd <- c(2.1, 2.2, 2.3)
pol_degrees <- c(1, 2, 3)

# Create polynomial powers and indexes correspondence matrix
pol_ind <- polynomialIndex(pol_degrees)
# Set all polynomial coefficients to 1
pol_coefficients <- rep(1, ncol(pol_ind))
pol_degrees_n <- length(pol_degrees)

# Assign coefficient 2 to the polynomial element (x1 ^ 1)*(x2 ^ 0)*(x3 ^ 2)
pol_coefficients[apply(pol_ind, 2, function(x) all(x == c(1, 0, 2)))] <- 2

# Visualize correspondence between polynomial 
# elements and their coefficients
as.data.frame(rbind(pol_ind, pol_coefficients),
   row.names = c("x1 power", "x2 power", "x3 power", "coefficients"),
   optional = TRUE)
printPolynomial(pol_degrees, pol_coefficients)

# Calculate density approximation 
   # at point x (note that x should be a matrix)
dhpa(x = x,
     pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd)
   # at points x and y
dhpa(x = x_y,
     pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd)

# Condition second component to be 0.5 i.e. X2 = 0.5.
# Substitute x and y second components with conditional value 0.5
x <- matrix(c(0.1, 0.5, 0.3), nrow = 1) # or simply x[2] <- 0.5
y <- matrix(c(0.4, 0.5, 0.6), nrow = 1) # or simply y[2] <- 0.5
x_y <- rbind(x, y) 
# Set TRUE to the second component indicating that it is conditioned
given_ind <- c(FALSE, TRUE, FALSE)

# Calculate conditional (on X2 = 0.5) density approximation 
   # at point x
dhpa(x = x,
     pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd,
     given_ind = given_ind)
   # at points x and y
dhpa(x = x_y,
     pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd,
     given_ind = given_ind)	
	
# Consider third component marginal distribution conditioned on the
# second component 0.5 value i.e. (X3 | X2 = 0.5).
# Set TRUE to the first component indicating that it is omitted
omit_ind <- c(TRUE, FALSE, FALSE)

# Calculate conditional (on x2 = 0.5) marginal (for x3) density approximation
   # at point x
dhpa(x = x,
     pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd,
     given_ind = given_ind, 
     omit_ind = omit_ind)
   # at points x and y
dhpa(x = x_y,
     pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd,
     given_ind = given_ind, 
     omit_ind = omit_ind)  
        
## Example demonstrating phpa function application.
## Let's approximate some three random variables (X1, X2, X3) 
## joint cumulative distribution function (cdf) at point (0,1, 0.2, 0.3)
## with Hermite polynomial of (1, 2, 3) degrees which polynomial 
## coefficients equal 1 except coefficient related to x1*(x^3) polynomial 
## element which equals 2. Also suppose that normal density related
## mean vector equals (1.1, 1.2, 1.3) while standard deviations
## vector is (2.1, 2.2, 2.3).

## Prepare initial values
x <- matrix(c(0.1, 0.2, 0.3), nrow = 1)
mean <- c(1.1, 1.2, 1.3)
sd <- c(2.1, 2.2, 2.3)
pol_degrees <- c(1, 2, 3)

# Create polynomial powers and indexes correspondence matrix
pol_ind <- polynomialIndex(pol_degrees)

# Set all polynomial coefficients to 1
pol_coefficients <- rep(1, ncol(pol_ind))
pol_degrees_n <- length(pol_degrees)

# Assign coefficient 2 to the polynomial element (x1 ^ 1)*(x2 ^ 0)*(x3 ^ 2)
pol_coefficients[apply(pol_ind, 2, function(x) all(x == c(1, 0, 2)))] <- 2

# Visualize correspondence between polynomial
# elements and their coefficients
as.data.frame(rbind(pol_ind, pol_coefficients),
	           row.names = c("x1 power", "x2 power", 
	                         "x3 power", "coefficients"),
          	 optional = TRUE)
printPolynomial(pol_degrees, pol_coefficients)

# Calculate cdf approximation at point x
phpa(x = x,
     pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd)
	
# Condition second component to be 0.5
# Substitute x second component with conditional value 0.5
x <- matrix(c(0.1, 0.5, 0.3), nrow = 1) # or simply x[2] <- 0.5

# Set TRUE to the second component indicating that it is conditioned
given_ind <- c(FALSE, TRUE, FALSE)

# Calculate conditional (on X2 = 0.5) cdf approximation at point x
phpa(x = x,
     pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd,
     given_ind = given_ind)
	
# Consider third component marginal distribution
# conditioned on the second component 0.5 value

# Set TRUE to the first component indicating that it is omitted
omit_ind <- c(TRUE, FALSE, FALSE)

# Calculate conditional (on X2 = 0.5) marginal (for X3) cdf 
# approximation at point x
phpa(x = x,
      pol_coefficients = pol_coefficients, 
      pol_degrees = pol_degrees,
      mean = mean, sd = sd,
      given_ind = given_ind, 
      omit_ind = omit_ind)

## Example demonstrating ihpa function application.
## Let's approximate some three random variables (X1, X2, X3) joint interval 
## distribution function (intdf) at lower and upper points (0,1, 0.2, 0.3) 
## and (0,4, 0.5, 0.6) correspondingly with Hermite polynomial of (1, 2, 3) 
## degrees which polynomial coefficients equal 1 except coefficient related 
## to x1*(x^3) polynomial element which equals 2. Also suppose that normal
## density related mean vector equals (1.1, 1.2, 1.3) while standard
## deviations vector is (2.1, 2.2, 2.3).

## Prepare initial values
x_lower <- matrix(c(0.1, 0.2, 0.3), nrow=1)
x_upper <- matrix(c(0.4, 0.5, 0.6), nrow=1)
mean <- c(1.1, 1.2, 1.3)
sd <- c(2.1, 2.2, 2.3)
pol_degrees <- c(1, 2, 3)

# Create polynomial powers and indexes correspondence matrix
pol_ind <- polynomialIndex(pol_degrees)

# Set all polynomial coefficients to 1
pol_coefficients <- rep(1, ncol(pol_ind))
pol_degrees_n <- length(pol_degrees)

# Assign coefficient 2 to the polynomial element (x1 ^ 1)*(x2 ^ 0)*(x3 ^ 2)
pol_coefficients[apply(pol_ind, 2, function(x) all(x == c(1, 0, 2)))] <- 2

# Visualize correspondence between polynomial 
# elements and their coefficients
as.data.frame(rbind(pol_ind, pol_coefficients),
	           row.names = c("x1 power", "x2 power", 
	                         "x3 power", "coefficients"),
              optional = TRUE)
printPolynomial(pol_degrees, pol_coefficients)

# Calculate intdf approximation at points x_lower and x_upper
ihpa(x_lower = x_lower, x_upper = x_upper, 
     pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd)
	
# Condition second component to be 0.7
# Substitute x second component with conditional value 0.7
x_upper <- matrix(c(0.4, 0.7, 0.6), nrow = 1)

# Set TRUE to the second component indicating that it is conditioned
given_ind <- c(FALSE, TRUE, FALSE)

# Calculate conditional (on X2 = 0.5) intdf approximation 
# at points x_lower and x_upper
ihpa(x_lower = x_lower, x_upper = x_upper,
     pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd,
     given_ind = given_ind)
	
# Consider third component marginal distribution
# conditioned on the second component 0.7 value
# Set TRUE to the first component indicating that it is omitted
omit_ind <- c(TRUE, FALSE, FALSE)

# Calculate conditional (on X2 = 0.5) marginal (for X3) 
# intdf approximation at points x_lower and x_upper
ihpa(x_lower = x_lower, x_upper = x_upper,
     pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd,
     given_ind = given_ind, omit_ind = omit_ind)
	
## Example demonstrating ehpa function application.
## Let's approximate some three random variables (X1, X2, X3) powered product 
## expectation for powers (3, 2, 1) with Hermite polynomial of (1, 2, 3) 
## degrees which polynomial coefficients equal 1 except coefficient 
## related to x1*(x^3) polynomial element which equals 2.
## Also suppose that normal density related mean vector equals 
## (1.1, 1.2, 1.3) while standard deviations vector is (2.1, 2.2, 2.3).

# Prepare initial values
expectation_powers = c(3,2,1)
mean <- c(1.1, 1.2, 1.3)
sd <- c(2.1, 2.2, 2.3)
pol_degrees <- c(1, 2, 3)

# Create polynomial powers and indexes correspondence matrix
pol_ind <- polynomialIndex(pol_degrees)

# Set all polynomial coefficients to 1
pol_coefficients <- rep(1, ncol(pol_ind))
pol_degrees_n <- length(pol_degrees)

#Assign coefficient 2 to the polynomial element (x1 ^ 1)*(x2 ^ 0)*(x3 ^ 2)
pol_coefficients[apply(pol_ind, 2, function(x) all(x == c(1, 0, 2)))] <- 2

# Visualize correspondence between polynomial elements and their coefficients
as.data.frame(rbind(pol_ind, pol_coefficients),
              row.names = c("x1 power", "x2 power", 
                            "x3 power", "coefficients"),
              optional = TRUE)
printPolynomial(pol_degrees, pol_coefficients)

# Calculate expected powered product approximation
ehpa(pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd, 
     expectation_powers = expectation_powers)
	
# Condition second component to be 0.5
# Substitute x second component with conditional value 0.5
x <- matrix(c(NA, 0.5, NA), nrow = 1)
#Set TRUE to the second component indicating that it is conditioned
given_ind <- c(FALSE, TRUE, FALSE)

# Calculate conditional (on X2 = 0.5) expected powered product approximation
ehpa(x = x,
     pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd, 
     expectation_powers = expectation_powers,
     given_ind = given_ind)
	
# Consider third component marginal distribution
# conditioned on the second component 0.5 value
# Set TRUE to the first component indicating that it is omitted
omit_ind <- c(TRUE, FALSE, FALSE)

# Calculate conditional (on X2 = 0.5) marginal (for X3) expected powered 
# product approximation at points x_lower and x_upper
ehpa(x = x,
     pol_coefficients = pol_coefficients, 
     pol_degrees = pol_degrees,
     mean = mean, sd = sd, 
     expectation_powers = expectation_powers,
     given_ind = given_ind, 
     omit_ind = omit_ind)
	
## Example demonstrating etrhpa function application.
## Let's approximate some three truncated random variables (X1, X2, X3) 
## powered product expectation for powers (3, 2, 1) with Hermite polynomial 
## of (1,2,3) degrees which polynomial coefficients equal 1 except 
## coefficient related to x1*(x^3) polynomial element which equals 2. Also
## suppose that normal density related mean vector equals (1.1, 1.2, 1.3) 
## while standard deviations vector is (2.1, 2.2, 2.3). Suppose that lower  
## and upper truncation points are (-1.1,-1.2,-1.3) and (1.1,1.2,1.3) 
## correspondingly.

# Prepare initial values
expectation_powers = c(3,2,1)
tr_left = matrix(c(-1.1,-1.2,-1.3), nrow = 1)
tr_right = matrix(c(1.1,1.2,1.3), nrow = 1)
mean <- c(1.1, 1.2, 1.3)
sd <- c(2.1, 2.2, 2.3)
pol_degrees <- c(1, 2, 3)

# Create polynomial powers and indexes correspondence matrix
pol_ind <- polynomialIndex(pol_degrees)
# Set all polynomial coefficients to 1
pol_coefficients <- rep(1, ncol(pol_ind))
pol_degrees_n <- length(pol_degrees)

# Assign coefficient 2 to the polynomial element (x1 ^ 1)*(x2 ^ 0)*(x3 ^ 2)
pol_coefficients[apply(pol_ind, 2, function(x) all(x == c(1, 0, 2)))] <- 2

# Visualize correspondence between polynomial elements and their coefficients
as.data.frame(rbind(pol_ind, pol_coefficients),
              row.names = c("x1 power", "x2 power", 
                            "x3 power", "coefficients"),
              optional = TRUE)
printPolynomial(pol_degrees, pol_coefficients)

# Calculate expected powered product approximation for truncated distribution
etrhpa(pol_coefficients = pol_coefficients, 
       pol_degrees = pol_degrees,
       mean = mean, sd = sd, 
       expectation_powers = expectation_powers,
       tr_left = tr_left, tr_right = tr_right)
       
## Example demonstrating dtrhpa function application.
## Let's approximate some three random variables (X1, X2, X3) joint density 
## function at point (0,1, 0.2, 0.3) with Hermite polynomial of (1,2,3)  
## degrees which polynomial coefficients equal 1 except coefficient related 
## to x1*(x^3) polynomial element which equals 2. Also suppose that normal 
## density related mean vector equals (1.1, 1.2, 1.3) while standard 
## deviations vector is (2.1, 2.2, 2.3). Suppose that lower and upper 
## truncation points are (-1.1,-1.2,-1.3) and (1.1,1.2,1.3) correspondingly.

# Prepare initial values
x <- matrix(c(0.1, 0.2, 0.3), nrow=1)
tr_left = matrix(c(-1.1,-1.2,-1.3), nrow = 1)
tr_right = matrix(c(1.1,1.2,1.3), nrow = 1)
mean <- c(1.1, 1.2, 1.3)
sd <- c(2.1, 2.2, 2.3)
pol_degrees <- c(1, 2, 3)

# Create polynomial powers and indexes correspondence matrix
pol_ind <- polynomialIndex(pol_degrees)

# Set all polynomial coefficients to 1
pol_coefficients <- rep(1, ncol(pol_ind))
pol_degrees_n <- length(pol_degrees)

# Assign coefficient 2 to the polynomial element (x1 ^ 1)*(x2 ^ 0)*(x3 ^ 2)
pol_coefficients[apply(pol_ind, 2, function(x) all(x == c(1, 0, 2)))] <- 2

# Visualize correspondence between polynomial elements and their coefficients
as.data.frame(rbind(pol_ind, pol_coefficients),
              row.names = c("x1 power", "x2 power", 
                            "x3 power", "coefficients"),
              optional = TRUE)
printPolynomial(pol_degrees, pol_coefficients)

# Calculate density approximation at point x
dtrhpa(x = x,
       pol_coefficients = pol_coefficients, 
       pol_degrees = pol_degrees,
       mean = mean, sd = sd,
       tr_left = tr_left, 
       tr_right = tr_right)
	
# Condition second component to be 0.5
# Substitute x second component with conditional value 0.5
x <- matrix(c(0.1, 0.5, 0.3), nrow = 1)
# Set TRUE to the second component indicating that it is conditioned
given_ind <- c(FALSE, TRUE, FALSE)
# Calculate conditional (on x2 = 0.5) density approximation at point x
dtrhpa(x = x,
       pol_coefficients = pol_coefficients, 
       pol_degrees = pol_degrees,
       mean = mean, sd = sd,
       given_ind = given_ind,
       tr_left = tr_left, tr_right = tr_right)
	
# Consider third component marginal distribution
# conditioned on the second component 0.5 value
# Set TRUE to the first component indicating that it is omitted
omit_ind <- c(TRUE, FALSE, FALSE)

# Calculate conditional (on X2 = 0.5) marginal (for X3) 
# density approximation at point x
dtrhpa(x = x,
       pol_coefficients = pol_coefficients, 
       pol_degrees = pol_degrees,
       mean = mean, sd = sd,
       given_ind = given_ind, omit_ind = omit_ind,
       tr_left = tr_left, tr_right = tr_right)
       
## Example demonstrating itrhpa function application.
## Let's approximate some three truncated random variables (X1, X2, X3) joint 
## interval distribution function at lower and upper points (0,1, 0.2, 0.3) 
## and (0,4, 0.5, 0.6) correspondingly with Hermite polynomial of (1 ,2, 3) 
## degrees which polynomial coefficients equal 1 except coefficient
## related to x1*(x^3) polynomial element which equals 2. Also suppose 
## that normal density related mean vector equals (1.1, 1.2, 1.3) while 
## standard deviations vector is (2.1, 2.2, 2.3). Suppose that lower and 
## upper truncation are (-1.1,-1.2,-1.3) and (1.1,1.2,1.3) correspondingly.

# Prepare initial values
x_lower <- matrix(c(0.1, 0.2, 0.3), nrow=1)
x_upper <- matrix(c(0.4, 0.5, 0.6), nrow=1)
tr_left = matrix(c(-1.1,-1.2,-1.3), nrow = 1)
tr_right = matrix(c(1.1,1.2,1.3), nrow = 1)
mean <- c(1.1, 1.2, 1.3)
sd <- c(2.1, 2.2, 2.3)
pol_degrees <- c(1, 2, 3)

# Create polynomial powers and indexes correspondence matrix
pol_ind <- polynomialIndex(pol_degrees)
# Set all polynomial coefficients to 1
pol_coefficients <- rep(1, ncol(pol_ind))
pol_degrees_n <- length(pol_degrees)

# Assign coefficient 2 to the polynomial element (x1 ^ 1)*(x2 ^ 0)*(x3 ^ 2)
pol_coefficients[apply(pol_ind, 2, function(x) all(x == c(1, 0, 2)))] <- 2
# Visualize correspondence between polynomial 
# elements and their coefficients
as.data.frame(rbind(pol_ind, pol_coefficients),
              row.names = c("x1 power", "x2 power", 
                            "x3 power", "coefficients"),
              optional = TRUE)
printPolynomial(pol_degrees, pol_coefficients)

# Calculate intdf approximation at points x_lower and x_upper
itrhpa(x_lower = x_lower, x_upper = x_upper, 
       pol_coefficients = pol_coefficients, 
       pol_degrees = pol_degrees,
       mean = mean, sd = sd,
       tr_left = tr_left, tr_right = tr_right)
    
# Condition second component to be 0.7
# Substitute x second component with conditional value 0.7
x_upper <- matrix(c(0.4, 0.7, 0.6), nrow = 1)
# Set TRUE to the second component indicating that it is conditioned
given_ind <- c(FALSE, TRUE, FALSE)

# Calculate conditional (on X2 = 0.5) intdf 
# approximation at points x_lower and x_upper
itrhpa(x_lower = x_lower, x_upper = x_upper,
       pol_coefficients = pol_coefficients, 
       pol_degrees = pol_degrees,
       mean = mean, sd = sd,
       given_ind = given_ind,
       tr_left = tr_left, tr_right = tr_right)
    
# Consider third component marginal distribution
# conditioned on the second component 0.7 value
# Set TRUE to the first component indicating that it is omitted
omit_ind <- c(TRUE, FALSE, FALSE)

# Calculate conditional (on X2 = 0.5) marginal (for X3) intdf 
# approximation at points x_lower and x_upper
itrhpa(x_lower = x_lower, x_upper = x_upper,
       pol_coefficients = pol_coefficients, 
       pol_degrees = pol_degrees,
       mean = mean, sd = sd,
       given_ind = given_ind, omit_ind = omit_ind,
       tr_left = tr_left, tr_right = tr_right)
       
## Example demonstrating dhpaDiff function application.
## Let's approximate some three random variables (X1, X2, X3) joint density
## function at point (0,1, 0.2, 0.3) with Hermite polynomial of (1,2,3)
## degrees which polynomial coefficients equal 1 except coefficient related
## to x1*(x^3) polynomial element which equals 2. Also suppose that normal
## density related mean vector equals (1.1, 1.2, 1.3) while standard
## deviations vector is (2.1, 2.2, 2.3). In this example let's calculate
## density approximating function's gradient respect to various parameters

# Prepare initial values
x <- matrix(c(0.1, 0.2, 0.3), nrow = 1)
mean <- c(1.1, 1.2, 1.3)
sd <- c(2.1, 2.2, 2.3)
pol_degrees <- c(1, 2, 3)

# Create polynomial powers and indexes correspondence matrix
pol_ind <- polynomialIndex(pol_degrees)
# Set all polynomial coefficients to 1
pol_coefficients <- rep(1, ncol(pol_ind))
pol_degrees_n <- length(pol_degrees)

# Assign coefficient 2 to the polynomial element (x1 ^ 1)*(x2 ^ 0)*(x3 ^ 2)
pol_coefficients[apply(pol_ind, 2, function(x) all(x == c(1, 0, 2)))] <- 2

# Visualize correspondence between polynomial
# elements and their coefficients
as.data.frame(rbind(pol_ind, pol_coefficients),
              row.names = c("x1 power", "x2 power", 
                            "x3 power", "coefficients"),
              optional = TRUE)
printPolynomial(pol_degrees, pol_coefficients)

# Calculate density approximation gradient 
# respect to polynomial coefficients at point x
dhpaDiff(x = x,
         pol_coefficients = pol_coefficients, 
         pol_degrees = pol_degrees,
         mean = mean, sd = sd)
	
# Condition second component to be 0.5
# Substitute x second component with conditional value 0.5
x <- matrix(c(0.1, 0.5, 0.3), nrow = 1)
# Set TRUE to the second component indicating that it is conditioned
given_ind <- c(FALSE, TRUE, FALSE)

# Calculate conditional (on x2 = 0.5) density approximation's 
# gradient respect to polynomial coefficients at point x
dhpaDiff(x = x,
         pol_coefficients = pol_coefficients, 
         pol_degrees = pol_degrees,
         mean = mean, sd = sd,
         given_ind = given_ind)
	
# Consider third component marginal distribution
# conditioned on the second component 0.5 value
# Set TRUE to the first component indicating that it is omitted
omit_ind <- c(TRUE, FALSE, FALSE)

# Calculate conditional (on X2 = 0.5) marginal (for X3) density 
# approximation's gradient respect to: 
  # polynomial coefficients
dhpaDiff(x = x,
         pol_coefficients = pol_coefficients, 
         pol_degrees = pol_degrees,
         mean = mean, sd = sd,
         given_ind = given_ind, 
         omit_ind = omit_ind)
  # mean
dhpaDiff(x = x,
         pol_coefficients = pol_coefficients, 
         pol_degrees = pol_degrees,
         mean = mean, sd = sd,
         given_ind = given_ind, 
         omit_ind = omit_ind,
         type = "mean")
  # sd
dhpaDiff(x = x,
         pol_coefficients = pol_coefficients, 
         pol_degrees = pol_degrees,
         mean = mean, sd = sd,
         given_ind = given_ind, 
         omit_ind = omit_ind,
         type = "sd")
 # x
dhpaDiff(x = x,
         pol_coefficients = pol_coefficients, 
         pol_degrees = pol_degrees,
         mean = mean, sd = sd,
         given_ind = given_ind, 
         omit_ind = omit_ind,
         type = "x")
         
## Example demonstrating ihpaDiff function application.
## Let's approximate some three random variables (X1, X2, X3 ) joint interval 
## distribution function (intdf) at lower and upper points (0,1, 0.2, 0.3) 
## and (0,4, 0.5, 0.6) correspondingly with Hermite polynomial of (1, 2, 3) 
## degrees which polynomial coefficients equal 1 except coefficient 
## related to x1*(x^3) polynomial element which equals 2.
## Also suppose that normal density related mean vector equals 
## (1.1, 1.2, 1.3) while standard deviations vector is (2.1, 2.2, 2.3).
## In this example let's calculate interval distribution approximating 
## function gradient respect to polynomial coefficients.

# Prepare initial values
x_lower <- matrix(c(0.1, 0.2, 0.3), nrow=1)
x_upper <- matrix(c(0.4, 0.5, 0.6), nrow=1)
mean <- c(1.1, 1.2, 1.3)
sd <- c(2.1, 2.2, 2.3)
pol_degrees <- c(1, 2, 3)

# Create polynomial powers and indexes correspondence matrix
pol_ind <- polynomialIndex(pol_degrees)

# Set all polynomial coefficients to 1
pol_coefficients <- rep(1, ncol(pol_ind))
pol_degrees_n <- length(pol_degrees)

# Assign coefficient 2 to the polynomial element (x1 ^ 1)*(x2 ^ 0)*(x3 ^ 2)
pol_coefficients[apply(pol_ind, 2, function(x) all(x == c(1, 0, 2)))] <- 2

# Visualize correspondence between polynomial 
# elements and their coefficients
as.data.frame(rbind(pol_ind, pol_coefficients),
              row.names = c("x1 power", "x2 power", 
                            "x3 power", "coefficients"),
              optional = TRUE)
printPolynomial(pol_degrees, pol_coefficients)

# Calculate intdf approximation gradient respect to 
# polynomial coefficients at points x_lower and x_upper
ihpaDiff(x_lower = x_lower, x_upper = x_upper, 
         pol_coefficients = pol_coefficients, 
         pol_degrees = pol_degrees,
         mean = mean, sd = sd)
	
# Condition second component to be 0.7
# Substitute x second component with conditional value 0.7
x_upper <- matrix(c(0.4, 0.7, 0.6), nrow = 1)

# Set TRUE to the second component indicating that it is conditioned
given_ind <- c(FALSE, TRUE, FALSE)

# Calculate conditional (on X2 = 0.5) intdf approximation
# respect to polynomial coefficients at points x_lower and x_upper
ihpaDiff(x_lower = x_lower, x_upper = x_upper,
         pol_coefficients = pol_coefficients, 
         pol_degrees = pol_degrees,
         mean = mean, sd = sd,
         given_ind = given_ind)
	
# Consider third component marginal distribution
# conditioned on the second component 0.7 value
# Set TRUE to the first component indicating that it is omitted
omit_ind <- c(TRUE, FALSE, FALSE)

# Calculate conditional (on X2 = 0.5) marginal (for X3) intdf approximation
# respect to:
  # polynomial coefficients
ihpaDiff(x_lower = x_lower, x_upper = x_upper,
         pol_coefficients = pol_coefficients, 
         pol_degrees = pol_degrees,
         mean = mean, sd = sd,
         given_ind = given_ind, omit_ind = omit_ind)
  # mean
ihpaDiff(x_lower = x_lower, x_upper = x_upper,
         pol_coefficients = pol_coefficients, 
         pol_degrees = pol_degrees,
         mean = mean, sd = sd,
         given_ind = given_ind, omit_ind = omit_ind,
         type = "mean")
  # sd
ihpaDiff(x_lower = x_lower, x_upper = x_upper,
         pol_coefficients = pol_coefficients, 
         pol_degrees = pol_degrees,
         mean = mean, sd = sd,
         given_ind = given_ind, omit_ind = omit_ind,
         type = "sd")
  # x_lower
ihpaDiff(x_lower = x_lower, x_upper = x_upper,
         pol_coefficients = pol_coefficients, 
         pol_degrees = pol_degrees,
         mean = mean, sd = sd,
         given_ind = given_ind, omit_ind = omit_ind,
         type = "x_lower")
  # x_upper
ihpaDiff(x_lower = x_lower, x_upper = x_upper,
          pol_coefficients = pol_coefficients, 
          pol_degrees = pol_degrees,
          mean = mean, sd = sd,
          given_ind = given_ind, omit_ind = omit_ind,
          type = "x_upper")
          
# }

Run the code above in your browser using DataLab