# R/mvorpb.R In mvprpb: Orthant Probability of the Multivariate Normal Distribution

#### Documented in mvorpb

```mvorpb <-
function( dim.p , m.tgt , v.tgt , n.itr , it.rg ){

#	Function test.orthant
#		Evaluate an orthant probability
#	Arguments:
#		dim.p: Dimension to evaluate (Integer, Scalar)
#		m.tgt: Mean (Real vector, length: dim.p)
#		v.tgt: Covariance Matrix (Real square matrix of size dim.p)
#		n.itr: Number of intervals for numerical integration
#		it.rg: Maximum point of the numerical integration range
#
#	In this function, two functions,  f.num.cumint.wdiff and f.comb.gen.01, are called, which are
#	defined in this file.

f.num.cumint.wdiff <-  function( f.val , diff.val , int.wid ){
#	Numericl Integration using Derivitives at Grid Points
#	This function is called from test.orthant

val.low <-   f.val[ - length(f.val) ]
val.up  <-   f.val[ - 1 ]
val.accm <-  c( 0 , 0.5*cumsum( ( val.low + val.up )*int.wid )
- cumsum(diff(diff.val)*int.wid^2)/12
)
return(val.accm)
}	# enf of function f.num.cumint.wdiff

f.comb.gen.01 <-  function( q.in , candi , pre.val ){
#	Generate Combinations
#	The technique of recursive call is used.
#	This function is called from test.orthant

if(q.in<=1){
rep.count <-  length(candi)
ret.val <-  rbind( matrix(rep(pre.val,rep.count),ncol=rep.count) , matrix(candi,nrow=1) )
}else{
if(length(candi)<=q.in){
ret.val <-  matrix( c(pre.val ,candi) , ncol=1)
}else{
ret.val <- cbind(	f.comb.gen.01(q.in-1 , candi[-1] , c(pre.val,candi[1])) ,
f.comb.gen.01(q.in   , candi[-1] , pre.val )   )
}}   # end if

return( ret.val )
} # enf of function f.comb.gen.01

#
#	Step 1
#		Initial Setting

scale.go <-  1 / sqrt( diag( v.tgt ) )
m.wrk <-  m.tgt * scale.go
v.go <-  diag( scale.go ) %*% v.tgt %*% diag( scale.go )

#	Cholesky Decomposition
#			t(b.vec) %*% b.vec === v.go

v.go <-  0.5 * ( v.go + t(v.go) )
b.vec <-  chol( v.go )

m.a  <-  solve(t(b.vec))%*%m.wrk
c.vec.pre <-  t(solve( b.vec ))

int.plane.vec <-  apply( b.vec , 1 , sum )
int.plane.vec <-   int.plane.vec / sqrt(sum(int.plane.vec^2))
c.vec <-  sweep( c.vec.pre , 2 , c.vec.adj , "/")

intg.lat.p <-  ( 0:n.itr )* it.rg / n.itr
intg.lat.d <-  diff( intg.lat.p )

wk.mat <-  matrix( NA , nrow=length(intg.lat.p) , ncol=(2^dim.p-1) )
wk.mat.diff <- wk.mat

wk.orthpnt.lensq <- rep( NA ,  2^dim.p-1)
wk.cross.mean <- rep( NA ,  2^dim.p-1)

#
#	Step 2
#		Evaluate vector w^J for each subspace

w.base <- matrix(NA , ncol=dim.p , nrow=dim.p )
w.plane.vec <- int.plane.vec
wk.orthpnt <- matrix(NA,nrow=dim.p,ncol=(2^dim.p))
wk.orthpnt[,2^dim.p-1] <- int.plane.vec
w.pos.comb <- 0
w.dim <- 0
i.w <- 0
w.stack <- rep(NA , dim.p)

while(w.dim >= 0){
i.w <- i.w + 1
w.u <- b.vec[ , i.w ]

if(w.dim>0){
w.seq <- 1:w.dim
w.coef1 <- t( w.base[,w.seq,drop=FALSE]) %*% matrix( w.u , ncol=1)
w.u <- w.u - w.base[,w.seq,drop=FALSE] %*% matrix(w.coef1 , ncol=1 )
}
w.u <- w.u / sqrt(sum(w.u^2))
w.coef2 <- sum( w.plane.vec * w.u )
w.pos.comb.new <- w.pos.comb + 2^(i.w - 1)
wk.orthpnt[ ,  2^dim.p -1- w.pos.comb.new ] <- w.plane.vec - w.coef2 * w.u

if( i.w < dim.p ){
w.plane.vec <-  wk.orthpnt[ ,  2^dim.p -1- w.pos.comb.new ]
w.dim <- w.dim + 1
w.stack[ w.dim ] <- i.w
w.base[ , w.dim ] <- w.u
w.pos.comb <- w.pos.comb.new

} else {
w.dim <- w.dim - 1
if(w.dim > -1){
i.w <- w.stack[ w.dim +1]
w.pos.comb <- w.pos.comb - 2^(i.w-1)
w.plane.vec <- wk.orthpnt[,2^dim.p-1-w.pos.comb]
}
} # end if
}	# end while

wk.orthpnt <- sweep(wk.orthpnt , 2,c(t(wk.orthpnt)%*% int.plane.vec ),"/"  )

#	Step 3
#		Evaluation of Two Dimensional Nodes

for( w.i in seq( 2 , dim.p )){
for( w.j in seq( 1 , w.i - 1 )){
w.vi <-  c.vec[ , w.i ]
w.vj <-  c.vec[ , w.j ]
w.v.dif <-  w.vi - w.vj
w.v.dif.len.s <-  sum( w.v.dif^2 )
w.v.dif.len  <-  sqrt( w.v.dif.len.s )

wk.1 <-  (-1) * sum( w.vj *w.v.dif )/w.v.dif.len.s
w.mean.posval <-  sum(w.v.dif * m.a )/w.v.dif.len
w.vec.base <-   wk.1*w.vi + (1-wk.1)*w.vj

w.speed.i <-  sum((w.vi - w.vec.base)*w.v.dif)/w.v.dif.len
w.speed.j <-  sum((w.vj - w.vec.base)*w.v.dif)/w.v.dif.len
w.pos.i <-  intg.lat.p * w.speed.i - w.mean.posval
w.pos.j <-  intg.lat.p * w.speed.j-  w.mean.posval

w.prb.sec <-  (pnorm(  w.pos.i ) - pnorm(  w.pos.j  ) )

w.prb.sec.diff <- (
w.speed.i * dnorm( w.pos.i )
-
w.speed.j * dnorm( w.pos.j )
)

w.respos <-  2^(w.i-1) + 2^(w.j - 1)
wk.mat[    , w.respos ] <-  w.prb.sec
wk.mat.diff[ , w.respos ] <-  w.prb.sec.diff

wk.orthpnt.lensq[ w.respos] <-  sum(w.vec.base^2)
wk.cross.mean[ w.respos] <- sum(  w.vec.base * m.a )

}}	#	next w.i , w.j

#	Step 4
#		Evaluation of Intermediate Dimension

if( dim.p > 2){
for(w.dim in seq(3,dim.p)){

w.pat.accm  <-  f.comb.gen.01(w.dim , seq(dim.p) , NULL)

for(w.i in seq(1 , ncol(w.pat.accm))){

w.pat.go <-  w.pat.accm[ , w.i ]
w.mat.x <-  c.vec[ , w.pat.go ]

w.prb.accm <-  0
w.prb.accm.diff <-  0

w.highdim.pos <-  sum( 2^(w.pat.go - 1 ) )
w.highdim.orthpnt <- wk.orthpnt[ , w.highdim.pos ]

wk.orthpnt.lensq[ w.highdim.pos ] <- sum( w.highdim.orthpnt^2 )
wk.cross.mean[ w.highdim.pos ] <-sum(  w.highdim.orthpnt * m.a )

for(w.j in seq(w.dim)){
#		Evaluate each integral in Step 4

w.valno.ineq <-  w.pat.go[ - w.j ]
w.lowdim.pos <-  sum(2^(w.valno.ineq-1))
w.lowdim.base <-  wk.orthpnt[,w.lowdim.pos]

w.shiftspeed <- sqrt( wk.orthpnt.lensq[w.lowdim.pos] - wk.orthpnt.lensq[ w.highdim.pos ]  )

if(abs(w.shiftspeed)>1e-6){
w.mean.posval <-  (wk.cross.mean[w.lowdim.pos]-wk.cross.mean[ w.highdim.pos ])/w.shiftspeed

w.val.lowdim <-  wk.mat[,w.lowdim.pos]
w.phi.arg <- intg.lat.p*w.shiftspeed-w.mean.posval
w.phi <- dnorm(w.phi.arg )
w.prb.intg.pre   <-  w.shiftspeed * w.val.lowdim * w.phi
w.prb.intg.diff.pre <-   w.shiftspeed * w.phi * (wk.mat.diff[ ,w.lowdim.pos ]
- w.phi.arg * w.shiftspeed * w.val.lowdim
)

w.prb.sec <-   f.num.cumint.wdiff( w.prb.intg.pre , w.prb.intg.diff.pre ,  intg.lat.d )

if(sum(b.vec[,w.pat.go[w.j]] * w.highdim.orthpnt)>0){
w.prb.accm <-  w.prb.accm + w.prb.sec
w.prb.accm.diff  <-  w.prb.accm.diff  + w.prb.intg.pre
} else {
w.prb.accm <-  w.prb.accm -w.prb.sec
w.prb.accm.diff  <-  w.prb.accm.diff  - w.prb.intg.pre

}
}	# end if
}	# next w.j

wk.mat[ , w.highdim.pos ] <-  w.prb.accm
wk.mat.diff[ , w.highdim.pos ] <-  w.prb.accm.diff

}	# next w.i
}}	# next w.dim  end if

#
#	Step 5
#		Highest Dimension

wk.respos <-  2^dim.p - 1
wk.3 <-  wk.mat[ , wk.respos ]
wk.mean.posval <-  sum( int.plane.vec * m.a )

ww.phi.arg <- intg.lat.p - wk.mean.posval
w.phi <- dnorm( ww.phi.arg )

w.prb.intg.pre   <-  wk.3* w.phi
w.prb.intg.diff.pre <-   w.phi * (wk.mat.diff[ , wk.respos ] - ww.phi.arg *  wk.3 )

res.prb <-   f.num.cumint.wdiff( w.prb.intg.pre ,w.prb.intg.diff.pre  , intg.lat.d)
ret.val <-  res.prb[ length(res.prb) ]
attr(ret.val , "error-itg-rg" ) <- pnorm( - (it.rg - w.mean.posval) )
return( ret.val  )
}	# enf of function test.orthant
```

## Try the mvprpb package in your browser

Any scripts or data that you put into this service are public.

mvprpb documentation built on May 2, 2019, 6:36 a.m.