MDI_solve: Solve MDI with linear equality constraints using a Monte...
In miguelbiron/mcMDISolver: Monte Carlo Minimum Discrimination Information Solver

Description Usage Arguments Details Value See Also Examples

View source: R/solver.R

This function solves a Minimum Discrimination Information problem with linear equality constraints, as described in the package vignette.

1 2	MDI_solve(X, f, ..., m, l_start = rep(0, k + 1), method = "Newton", control = list(trace = 1))

`X`	an S x n matrix containing S samples from p
`f`	a function returning a vector of size k that evaluates the LHS of restrictions
`m`	a vector of size k with the RHS of the restrictions
`l_start`	a guess of the optimal k + 1 lambdas. Default is `rep(0, k + 1)`
`control`	a list of parameters passed to `nleqslv`

The non - linear set of equations is solved using package nleqslv. By default, the method is set to Newton, as this usually results in faster convergences. The rest of the parameters of the solver that are not passed through the control list are left in their default values.

This method assumes that you are able to sample from the p distribution in order to produce the matrix X. If this is not the case, then a robust alternative is to use an MCMC method (such as the Metropolis algorithm), assuming that you can at least evaluate the density of p.

The function f should take as input a vector x of size n (the dimension of the space where x lives) and return a vector of size k (the number of non-trivial conditions), such that the first component of f(x) corresponds to the value of the first restriction function evaluated at x, the second component to the second restriction, and so on.

A list as returned by nleqslv

nleqslv for details on this non - linear equation solver, and mcmc for an implementation of the Metropolis algorithm.

#################################################################
# Toy example:
#
# p: multivariate normal (0, \Sigma)
# Restrictions: force 0.5 mass of the marginals to lie above some
#               randomly chosen thresholds (chi).
# The above implies:
#                   n = k
#                   f(x) = (x > chi)
#                   m = rep(0.5, n)
#################################################################

library(mcMDISolver)

set.seed(1313) # for reproducibility

n = 9 # number of variables
k = n # number of restrictions
S = 50000 # number of samples

# restrictions
chi = rnorm(n, 1.96, 0.1) # critical values
f = function(x, c = chi){ # vectorized over the size of x (n)
  return(x > c)
}
m = rep(.5, k) # rhs of restrictions

# p-sampler
# p is mvnorm with average cor = rho
rho = 0.6 # average correlation
R = diag(n)
R[upper.tri(R)] = pmin(pmax(rnorm(n*(n-1)/2, rho, 0.01), -1), 1)
R[lower.tri(R)] = t(R)[lower.tri(R)]

## SOLVE: serial implementation
start.time = proc.time()
fit_MDI = mcMDISolver::MDI_solve(mvtnorm::rmvnorm(n = S, sigma = R), f = f, m = m)
print(proc.time() - start.time)
print(fit_MDI$x)

## Examine q distribution
# Draw samples using Metropolis algorithm

# log density of q
log_dq = function(x, l = fit_MDI$x){
  return(mvtnorm::dmvnorm(x, sigma = R, log = T) - as.numeric(crossprod(l, c(1, f(x)))))
}

# draw from q
mcmc_S = 2 * S
burn = trunc(0.1 * mcmc_S) # number of samples to burn
sample_q = mcmc::metrop(log_dq, initial = chi, nbatch = mcmc_S + burn)
X_q = sample_q$batch[-(1:burn), ]

# histograms
# note the sudden increase in density around chi
par(mfrow = c(3, 3))
for(i in 1:n) hist(X_q[,i], main = paste("chi =", round(chi[i], 2)), xlab = "")
par(mfrow = c(1, 1))

# check restrictions close to zero
rowMeans(apply(X_q, 1, f)) - m