¡@
Some results on the truncated multivariate t distribution
by Hsiu J. Ho, Tsung-I Lin, Hsuan-Yu Chen and Wan-Lun Wang
Journal of Statistical Planning and Inference (2012), v142, pp.25-40
# R commands: slice sampling for the truncated multivariate t (TMVT) distribution
TT.GS = function(n,mu,S,nu,lower,upper)
{
require(mvtnorm)
p=length(mu)
TT.GS.sp = function(n,R,nu,a,b)
{
#initial value by using rejection sampling
al = pmvt(lower=a,upper=b,sigma=R,df=nu)[1]
repeat
{
x = rmvt(round(1/al),sigma=R,df=nu)
index = colSums(t(x)>a & t(x)<b) == p
if(sum(index)>0) break
}
x = x[which(index)[1],]
if(n<1) return(t(x))
#######
X = matrix(NA, n, p)
R.inv = solve(R)
for(i in 1:n)
{
delta = sum(colSums(x*R.inv)*x)
y = runif(1,0,exp(-.5*(nu+p)*log(1+delta/nu)))
kap = nu*(y^(-2/(nu+p))-1)
for(j in 1:p)
{
ss = x[-j]%*%R.inv[-j,-j]%*%x[-j]
mj = - sum(R.inv[-j,j]*x[-j]) / R.inv[j,j]
tj = sqrt(mj^2 + (kap - ss) / R.inv[j,j])
xij = runif(1,max(a[j],mj-tj),min(b[j],mj+tj))
X[i,j] = xij
x[j] = xij
}
}
return(X)
}
s = sqrt(diag(S))
R = S/outer(s,s,"*")
Z = TT.GS.sp(n,R,nu,a=(lower-mu)/s,b=(upper-mu)/s)
X = t(mu + t(Z) * s)
return(X)
}
# R commands: calculation of the first two moments of the TMVT distribution
TT.moment = function(a,b,R,nu)
{
require(mvtnorm)
GB = GenzBretz(maxpts = 5e4, abseps = 1e-9, releps = 0)
p = length(a)
a = ifelse(a==-Inf,rep(-1e12,p),a)
b = ifelse(b== Inf,rep( 1e12,p),b)
al0 = pmvt(lower = a, upper = b, sigma = R, df = nu, algorithm = GB)[1]
### pdf & cdf
la1 = (nu-2)/nu; la2 = (nu-4)/nu
da = (nu-1)/(nu+a^2); db = (nu-1)/(nu+b^2)
f1a = sqrt(la1)*dt(sqrt(la1)*a,df=nu-2)
f1b = sqrt(la1)*dt(sqrt(la1)*b,df=nu-2)
f2 = matrix(NA, p, p)
G1a = G1b = rep(NA, p)
G2 = matrix(NA, p, p)
for(r in 1:p)
{
temp = R[-r,r]
S1 = R[-r,-r] - temp %*% t(R[r,-r])
mua = temp * a[r]; low = a[-r]-mua; upp = b[-r]-mua
G1a[r] = ifelse(p==2,pt(upp/sqrt(S1/da[r]),df=nu-1)-pt(low/sqrt(S1/da[r]),df=nu-1)
,pmvt(lower = low, upper = upp, sigma = S1/da[r], df = nu-1, algorithm = GB)[1])
mub = temp * b[r]; low = a[-r]-mub; upp = b[-r]-mub
G1b[r] = ifelse(p==2,pt(upp/sqrt(S1/db[r]),df=nu-1)-pt(low/sqrt(S1/db[r]),df=nu-1)
,pmvt(lower = low, upper = upp, sigma = S1/db[r], df = nu-1, algorithm = GB)[1])
}
qa = f1a*G1a; qb = f1b*G1b
EX = c(R %*% (qa-qb)) / al0 / la1
H = matrix(0,p,p)
for(r in 1:(p-1))
{
for(s in (r+1):p)
{
rs = c(r,s)
pdf.aa = dmvt(c(a[r],a[s]),sigma=R[rs,rs]/la2,df=nu-4, log =F)
pdf.ab = dmvt(c(a[r],b[s]),sigma=R[rs,rs]/la2,df=nu-4, log =F)
pdf.ba = dmvt(c(b[r],a[s]),sigma=R[rs,rs]/la2,df=nu-4, log =F)
pdf.bb = dmvt(c(b[r],b[s]),sigma=R[rs,rs]/la2,df=nu-4, log =F)
if(p==2){cdf.aa=cdf.ab=cdf.ba=cdf.bb=1}
if(p>2)
{
tmp = R[-rs,rs]%*%solve(R[rs,rs])
mu.aa = c(tmp%*%c(a[r],a[s]))
mu.ab = c(tmp%*%c(a[r],b[s]))
mu.ba = c(tmp%*%c(b[r],a[s]))
mu.bb = c(tmp%*%c(b[r],b[s]))
daa = (nu-2)/(nu+(a[r]^2-2*R[r,s]*a[r]*a[s]+a[s]^2)/(1-R[r,s]^2))
dab = (nu-2)/(nu+(a[r]^2-2*R[r,s]*a[r]*b[s]+b[s]^2)/(1-R[r,s]^2))
dba = (nu-2)/(nu+(b[r]^2-2*R[r,s]*b[r]*a[s]+a[s]^2)/(1-R[r,s]^2))
dbb = (nu-2)/(nu+(b[r]^2-2*R[r,s]*b[r]*b[s]+b[s]^2)/(1-R[r,s]^2))
R21 = R[-rs,-rs] - R[-rs,rs]%*%solve(R[rs,rs]) %*% R[rs,-rs]
cdf.aa = ifelse(p==3,pt((b[-rs]-mu.aa)/sqrt(R21/daa),df=nu-2)-pt((a[-rs]-mu.aa)/sqrt(R21/daa),df=nu-2)
,pmvt(lower = a[-rs]-mu.aa, upper = b[-rs]-mu.aa, sigma = R21/daa, df=nu-2, algorithm = GB)[1])
cdf.ab = ifelse(p==3,pt((b[-rs]-mu.ab)/sqrt(R21/dab),df=nu-2)-pt((a[-rs]-mu.ab)/sqrt(R21/dab),df=nu-2)
,pmvt(lower = a[-rs]-mu.ab, upper = b[-rs]-mu.ab, sigma = R21/dab, df=nu-2, algorithm = GB)[1])
cdf.ba = ifelse(p==3,pt((b[-rs]-mu.ba)/sqrt(R21/dba),df=nu-2)-pt((a[-rs]-mu.ba)/sqrt(R21/dba),df=nu-2)
,pmvt(lower = a[-rs]-mu.ba, upper = b[-rs]-mu.ba, sigma = R21/dba, df=nu-2, algorithm = GB)[1])
cdf.bb = ifelse(p==3,pt((b[-rs]-mu.bb)/sqrt(R21/dbb),df=nu-2)-pt((a[-rs]-mu.bb)/sqrt(R21/dbb),df=nu-2)
,pmvt(lower = a[-rs]-mu.bb, upper = b[-rs]-mu.bb, sigma = R21/dbb, df=nu-2, algorithm = GB)[1])
}
H[r,s] = H[s,r] = pdf.aa*cdf.aa - pdf.ab*cdf.ab - pdf.ba*cdf.ba + pdf.bb*cdf.bb
}
}
H = H / la2
D = matrix(0,p,p)
diag(D) = a * qa - b * qb - diag(R%*%H)
al1 = pmvt(lower = a, upper = b, sigma = R/la1, df=nu-2, algorithm = GB)[1]
EXX = (al1 * R + R %*% (H + D) %*% R) / al0 / la1
return(list(EX=EX,EXX=EXX))
}
# A test example
rho=0.9
S=matrix(c(1, rho ,rho, 1),2,2)
nu=5
p=2
a=c(-2,2)
b=c(1,3)
mu = rep(0, p)
Y= TT.GS(n=10000, mu, S, nu, lower=a, upper=b)
M.Y=TT.moment(a, b, R=S, nu)
# First Moment
M.Y$EX
y.bar=colMeans(Y)
y.bar
# Second Moment
M.Y$EXX-M.Y$EX%*%t(M.Y$EX)
S.y=cov(Y)
S.y
¡@
¡@
¡@