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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
|
# Density function for the truncated multivariate t-distribution
#
# Author: stefan
###############################################################################
# Density function for the truncated multivariate t-distribution
# @param x
# @param mean
# @param sigma
# @param df degrees of freedom parameter
# @param log
dtmvt <- function(x, mean = rep(0, nrow(sigma)), sigma = diag(length(mean)), df = 1, lower= rep( -Inf, length = length(mean)),
upper = rep( Inf, length = length(mean)), log = FALSE){
# Check of additional inputs like x
if (is.vector(x)) {
x <- matrix(x, ncol = length(x))
}
# Anzahl der Beobachtungen
T = nrow(x)
# check for each row if in support region
insidesupportregion <- logical(T)
for (i in 1:T)
{
insidesupportregion[i] = all(x[i,] >= lower & x[i,] <= upper & !any(is.infinite(x)))
}
# density value for points outside the support region
dv = if (log) { -Inf } else { 0 }
# conditional density
f <- ifelse(insidesupportregion,
dmvt(x, delta=mean, sigma=sigma, df=df, log=log) / pmvt(lower=lower, upper=upper, delta=mean, sigma=sigma, df=df, type="shifted"), dv)
return(f)
}
if (FALSE) {
# Example
x1<-seq(-2, 3, by=0.1)
x2<-seq(-2, 3, by=0.1)
mean=c(0,0)
sigma=matrix(c(1, -0.5, -0.5, 1), 2, 2)
lower=c(-1,-1)
density<-function(x)
{
z=dtmvt(x, mean=mean, sigma=sigma, lower=lower)
z
}
fgrid <- function(x, y, f)
{
z <- matrix(nrow=length(x), ncol=length(y))
for(m in 1:length(x)){
for(n in 1:length(y)){
z[m,n] <- f(c(x[m], y[n]))
}
}
z
}
# compute multivariate-t density d for grid
d=fgrid(x1, x2, function(x) dtmvt(x, mean=mean, sigma=sigma, lower=lower))
# compute multivariate normal density d for grid
d2=fgrid(x1, x2, function(x) dtmvnorm(x, mean=mean, sigma=sigma, lower=lower))
# plot density as contourplot
contour(x1, x2, d, nlevels=5, main="Truncated Multivariate t Density",
xlab=expression(x[1]), ylab=expression(x[2]))
contour(x1, x2, d2, nlevels=5, add=TRUE, col="red")
abline(v=-1, lty=3, lwd=2)
abline(h=-1, lty=3, lwd=2)
}
|
