File: varcov.R

package info (click to toggle)
r-cran-eco 4.0-1-3
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 772 kB
  • sloc: ansic: 4,214; makefile: 2
file content (71 lines) | stat: -rw-r--r-- 1,816 bytes parent folder | download | duplicates (7)
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
varcov <- function(object, ...)
  UseMethod("varcov")

varcov.eco <- function(object, subset = NULL, ...) {
  if (is.null(subset))
    subset <- 1:nrow(object$Sigma)
  else if (max(subset) > nrow(object$Sigma))
    stop(paste("invalid input for `subset.' only", nrow(object$Sigma),
               "draws are stored.")) 

  p <- ncol(object$mu)
  n <- length(subset)
  Sigma <- array(0, c(p, p, n))
  cov <- object$Sigma[subset,]

  for (i in 1:n) {
    count <- 1
    for (j in 1:p) {
      Sigma[j,j:p,i] <- cov[i,count:(count+p-j)]
      count <- count + p - j + 1
    }
    diag(Sigma[,,i]) <- diag(Sigma[,,i]/2)
    Sigma[,,i] <- Sigma[,,i] + t(Sigma[,,i])
  }
  if (n > 1)
    return(Sigma)
  else
    return(Sigma[,,1])  
}

varcov.ecoNP <- function(object, subset = NULL, obs = NULL, ...) {
  if (is.null(subset))
    subset <- 1:nrow(object$Sigma)
  else if (max(subset) > nrow(object$Sigma))
    stop(paste("invalid input for `subset.' only", nrow(object$Sigma),
               "draws are stored.")) 

  if (is.null(obs))
    obs <- 1:dim(object$Sigma)[3]
  else if (max(subset) > dim(object$Sigma)[3])
    stop(paste("invalid input for `obs.' only", dim(object$Sigma)[3],
               "draws are stored."))
  
  p <- ncol(object$mu)
  n <- length(subset)
  m <- length(obs)
  Sigma <- array(0, c(p, p, n, m))
  cov <- object$Sigma[subset,,obs]

  for (k in 1:m) {
    for (i in 1:n) {
      count <- 1
      for (j in 1:p) {
        Sigma[j,j:p,i,k] <- cov[i,count:(count+p-j),k]
        count <- count + p - j + 1
      }
      diag(Sigma[,,i,k]) <- diag(Sigma[,,i,k]/2)
      Sigma[,,i,k] <- Sigma[,,i,k] + t(Sigma[,,i,k])
    }
  }
  if (n > 1)
    if (m > 1)
      return(Sigma)
    else
      return(Sigma[,,,1])
  else
    if (m > 1)
      return(Sigma[,,1,])
    else
      return(Sigma[,,1,1])
}