lsConstrain.fit: Minimize Inequality Constrained Mahalanobis Distance

Description Usage Arguments Value References Examples

View source: R/lsConstrain.fit.q

Description

Find the vector z that solves:

min{ (x - z)'inv(S)(x - z); Az <= b },

where x is an input vector, S its covariance matrix, A is a matrix of known contrasts, and b is a vector of known constraint constants.

Usage

1
lsConstrain.fit(x, b, s, a, iflag, itmax=4000, eps=1e-06, eps2=1e-06)

Arguments

x

vector of length n

b

vector of length k, containing constraint constants

s

matrix of dim n x n, the covariance matrix for vector x

a

matrix of dim k x n, for the contraints

iflag

vector of length k; an item = 0 if inequality constraint, 1 if equality constraint

itmax

scalar for number of max interations

eps

scalar of accuracy for convergence

eps2

scalar to determine close to zero

Value

List with the following components:

itmax: (defined above)

eps: (defined above)

eps2: (defined above)

iflag: (defined above)

xkt: vector of length k, for the Kuhn-Tucker coefficients.

iter: number of completed iterations.

supdif: greatest difference between estimates across a full cycle

ifault: error indicator: 0 = no error 1 = itmax exceeded 3 = invalid constraint function for some row ASA'=0.

a: (defined above)

call: function call

x.init: input vector x.

x.final: the vector "z" that solves the equation (see z in description).

s: (defind above)

min.dist: the minimum value of the function (see description).

References

Wollan PC, Dykstra RL. Minimizing inequality constrained mahalanobis distances. Applied Statistics Algorithm AS 225 (1987).

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
# An simulation example with linear regression with 3 beta's, 
# where we have the contraints:
#
# b[1] > 0
# b[2] - b[1] < 0
# b[3] < 0


set.seed(111)

n <- 100
x <- rep(1:3,rep(n,3))
x <- 1*outer(x,1:3,"==")

beta <- c(2,1,1)

y <- x%*%beta + rnorm(nrow(x))

fit <- lm(y ~-1 + x)

s <- solve( t(x) %*% x )

bhat <- fit$coef


a <-  rbind(c(-1, 0, 0),
            c(-1, 1, 0),
            c( 0, 0, 1))

# View expected constraints (3rd not met):

a %*% bhat
#            [,1] 
# [1,] -1.8506811
# [2,] -0.9543320
# [3,]  0.8590827

b <- rep(0, nrow(a))
iflag <- rep(0,length(b))

save <- lsConstrain.fit(x=bhat,b=b, s=s, a=a, iflag=iflag, itmax=500, 
                        eps=1e-6, eps2=1e-6)

save

Example output

          [,1]
[1,] -1.987181
[2,] -1.004048
[3,]  1.136383
$itmax
[1] 500

$eps
[1] 1e-06

$eps2
[1] 1e-06

$b
[1] 0 0 0

$iflag
[1] 0 0 0

$xkt
[1]   0.0000   0.0000 227.2766

$iter
[1] 2

$supdif
[1] 0

$ifault
[1] 0

$a
     [,1] [,2] [,3]
[1,]   -1    0    0
[2,]   -1    1    0
[3,]    0    0    1

$call
lsConstrain.fit(x = bhat, b = b, s = s, a = a, iflag = iflag, 
    itmax = 500, eps = 1e-06, eps2 = 1e-06)

$x.init
[1] 1.9871809 0.9831329 1.1363830

$x.final
[1] 1.9871809 0.9831329 0.0000000

$s
     [,1] [,2] [,3]
[1,] 0.01 0.00 0.00
[2,] 0.00 0.01 0.00
[3,] 0.00 0.00 0.01

$min.dist
         [,1]
[1,] 129.1366

ibdreg documentation built on Sept. 5, 2021, 5:53 p.m.