← Back to branch summary

~ubuntu-branches/ubuntu/quantal/mnormt/quantal

« back to all changes in this revision

Viewing changes to R/mnormt.R

  • Committer: Bazaar Package Importer
  • Author(s): Dirk Eddelbuettel
  • Date: 2007-10-03 21:25:21 UTC
  • Revision ID: james.westby@ubuntu.com-20071003212521-e6zr9fhv1oowj396
Tags: upstream-1.2-1
Import upstream version 1.2-1

Show diffs side-by-side

added added

removed removed

Lines of Context:
R/mnormt.R
 
1
 
 
2
# if(.Platform$OS.type=="unix") dyn.load("mnormt/src/sadmvnt.so")
 
3
# if(.Platform$OS.type=="windows") dyn.load("mnormt/src/sadmvnt.dll")
 
4
 
 
5
 
 
6
dmnorm <- function(x, mean=rep(0,d), varcov, log=FALSE)
 
7
{
 
8
  d  <- if(is.matrix(varcov)) ncol(varcov) else 1
 
9
  if(d>1 & is.vector(x)) x <- matrix(x, 1, d)
 
10
  n  <- if(d==1)  length(x) else nrow(x) 
 
11
  X  <- t(matrix(x, nrow=n, ncol=d)) - mean
 
12
  Q  <- apply((solve(varcov)%*% X)* X, 2, sum) 
 
13
  logDet <- sum(logb(abs(diag(qr(varcov)[[1]]))))
 
14
  logPDF <- as.vector(Q + d*logb(2*pi)+logDet)/(-2)
 
15
  if(log) logPDF else exp(logPDF)
 
16
}
 
17
 
 
18
 
 
19
rmnorm <- function(n=1, mean=rep(0,d), varcov)
 
20
 {
 
21
  d <- if(is.matrix(varcov)) ncol(varcov) else 1
 
22
  z <- matrix(rnorm(n*d),n,d) %*% chol(varcov)
 
23
  y <- t(mean+t(z))
 
24
  return(y)
 
25
 }
 
26
 
 
27
 
 
28
pmnorm <- function(x, mean=rep(0,length(x)), varcov, ...)
 
29
   sadmvn(lower=rep(-Inf, length(x)), upper=x, mean, varcov, ...)
 
30
 
 
31
 
 
32
sadmvn <- function(lower, upper, mean, varcov,  
 
33
                   maxpts=2000*d, abseps=1e-6, releps=0)
 
34
{
 
35
  if(any(lower > upper)) stop("lower>upper integration limits")
 
36
  if(any(lower == upper)) return(0)
 
37
  d <- as.integer(if(is.matrix(varcov)) ncol(varcov) else 1)
 
38
  varcov <- matrix(varcov, d, d)
 
39
  sd  <- sqrt(diag(varcov))
 
40
  rho <- diag(1/sd,d,d) %*% varcov %*% diag(1/sd,d,d)
 
41
  lower <- as.double((lower-mean)/sd)
 
42
  upper <- as.double((upper-mean)/sd)
 
43
  infin <- rep(2,d)
 
44
  infin <- replace(infin, (upper == Inf) & (lower > -Inf), 1)
 
45
  infin <- replace(infin, (upper < Inf) & (lower == -Inf), 0)
 
46
  infin <- replace(infin, (upper == Inf) & (lower == -Inf), -1)
 
47
  infin <- as.integer(infin)
 
48
  lower <- replace(lower, lower == -Inf, 0)
 
49
  upper <- replace(upper, upper == Inf, 0)
 
50
  correl <- as.double(rho[upper.tri(rho, diag=FALSE)])
 
51
  maxpts <- as.integer(maxpts)
 
52
  abseps <- as.double(abseps)
 
53
  releps <- as.double(releps)
 
54
  error  <- as.double(0)
 
55
  value  <- as.double(0)
 
56
  inform <- as.integer(0)
 
57
  result <- .Fortran("sadmvn", d, lower, upper, infin, correl, maxpts,
 
58
               abseps, releps, error, value, inform, PACKAGE="mnormt")
 
59
  prob <- result[[10]]
 
60
  attr(prob,"error")  <- result[[9]]
 
61
  attr(prob,"status") <- switch(1+result[[11]], 
 
62
                "normal completion", "accuracy non achieved", "oversize")
 
63
  return(prob)
 
64
}
 
65
 
 
66
#----
 
67
 
 
68
dmt <- function (x, mean=rep(0,d), S, df = Inf, log = FALSE) 
 
69
{
 
70
  if (df == Inf)  return(dmnorm(x, mean, S, log = log))
 
71
  d  <- if(is.matrix(S)) ncol(S) else 1
 
72
  if(d>1 & is.vector(x)) x <- matrix(x, 1, d)
 
73
  n  <- if(d==1) length(x) else nrow(x) 
 
74
  X <- t(matrix(x, nrow = n, ncol = d)) - mean
 
75
  Q <- apply((solve(S) %*% X) * X, 2, sum)
 
76
  logDet <- sum(logb(abs(diag(qr(S)$qr))))
 
77
  logPDF <- (lgamma((df + d)/2) - 0.5 * (d * logb(pi * df) + logDet)
 
78
             - lgamma(df/2) - 0.5 * (df + d) * logb(1 + Q/df))
 
79
  if(log) logPDF else exp(logPDF)
 
80
}
 
81
 
 
82
rmt <- function(n=1, mean=rep(0,d), S, df=Inf)
 
83
{ 
 
84
  d <- if(is.matrix(S)) ncol(S) else 1 
 
85
  if(df==Inf) x <-1 else x <- rchisq(n,df)/df
 
86
  z <- rmnorm(n, rep(0,d), S)
 
87
  y <- t(mean + t(z/sqrt(x)))
 
88
  return(y)
 
89
}
 
90
 
 
91
 
 
92
pmt <- function(x, mean=rep(0,length(x)), S, df=Inf, ...)
 
93
   sadmvt(df, lower=rep(-Inf, length(x)), upper=x, mean, S, ...)
 
94
 
 
95
 
 
96
sadmvt <- function(df, lower, upper, mean, S, 
 
97
                   maxpts=2000*d, abseps=1e-6, releps=0)
 
98
{
 
99
  if(df == Inf) return(sadmvn(lower, upper, mean, S, maxpts, abseps, releps))
 
100
  if(any(lower > upper)) stop("lower>upper integration limits")
 
101
  if(any(lower == upper)) return(0)
 
102
  if(round(df) != df) warning("non integer df is rounded to integer") 
 
103
  df <- as.integer(round(df))
 
104
  d  <- as.integer(if(is.matrix(S)) ncol(S) else 1)
 
105
  S  <- matrix(S, d, d)
 
106
  sd  <- sqrt(diag(S))
 
107
  rho <- diag(1/sd,d,d) %*% S %*% diag(1/sd,d,d)
 
108
  lower <- as.double((lower-mean)/sd)
 
109
  upper <- as.double((upper-mean)/sd)
 
110
  infin <- rep(2,d)
 
111
  infin <- replace(infin, (upper == Inf) & (lower > -Inf), 1)
 
112
  infin <- replace(infin, (upper < Inf) & (lower == -Inf), 0)
 
113
  infin <- replace(infin, (upper == Inf) & (lower == -Inf), -1)
 
114
  infin <- as.integer(infin)
 
115
  lower <- replace(lower, lower == -Inf, 0)
 
116
  upper <- replace(upper, upper == Inf, 0)
 
117
  correl <- rho[upper.tri(rho, diag=FALSE)]
 
118
  maxpts <- as.integer(maxpts)
 
119
  abseps <- as.double(abseps)
 
120
  releps <- as.double(releps)
 
121
  error  <- as.double(0)
 
122
  value  <- as.double(0)
 
123
  inform <- as.integer(0)
 
124
  result <- .Fortran("sadmvt", d, df, lower, upper, infin, correl, maxpts,
 
125
                   abseps, releps, error, value, inform, PACKAGE="mnormt")
 
126
  prob <- result[[11]]
 
127
  attr(prob,"error")  <- result[[10]]
 
128
  attr(prob,"status") <- switch(1+result[[12]], 
 
129
                "normal completion", "accuracy non achieved", "oversize")
 
130
  return(prob)
 
131
}
 
132
 
 
133
.First.lib <- function(library, pkg)
 
134
{ 
 
135
   Rv <- R.Version()
 
136
   if(Rv$major < 2 |(Rv$major == 2 && Rv$minor < 2.0))
 
137
        stop("This package requires R 2.2.0 or later")
 
138
   library.dynam("mnormt", pkg, library)
 
139
   invisible()
 
140
}
 
141
 
 
142
 
 
143