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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
|
# The following functions were originally in the matrixcalc package
# written by Frederick Novomestky
# and distributed under the GPL license version 2 or higher
# last modified 2021-06-02
duplication.matrix <- function (n = 1) {
return(D.matrix(n))
}
D.matrix <- function (n) {
if (missing(n))
stop("argument n is missing")
if (!is.numeric(n))
stop("argument n is not numeric")
if (n != trunc(n))
stop("argument n is not an integer")
if (n < 2)
stop("argument n is less than 2")
p <- n * (n + 1)/2
nsq <- n * n
Dt <- matrix(0, nrow = p, ncol = nsq)
T <- T.matrices(n)
u <- u.vectors(n)
k <- u$k
I <- u$I
for (j in 1:n) {
for (i in j:n) {
Dt <- Dt + I[, k[i, j]] %*% t(vec(T[[i]][[j]]))
}
}
return(t(Dt))
}
vech <- function (x) {
if (!is.square.matrix(x))
stop("argument x is not a square numeric matrix")
return(t(t(x[!upper.tri(x)])))
}
is.square.matrix <- function (x) {
if (!is.matrix(x))
stop("argument x is not a matrix")
return(nrow(x) == ncol(x))
}
T.matrices <- function (n) {
if (missing(n))
stop("argument n is missing")
if (!is.numeric(n))
stop("argument n is not numeric")
if (n != trunc(n))
stop("argument n is not an integer")
if (n < 2)
stop("argument n is less than 2")
E <- E.matrices(n)
T <- list()
for (i in 1:n) {
T[[i]] <- list()
for (j in 1:n) {
if (i == j) {
T[[i]][[j]] <- E[[i]][[j]]
}
else {
T[[i]][[j]] <- E[[i]][[j]] + E[[j]][[i]]
}
}
}
return(T)
}
u.vectors <- function (n) {
if (n != trunc(n))
stop("argument n is not an integer")
if (n < 2)
stop("argument n is less than 2")
p <- n * (n + 1)/2
I <- diag(rep(1, p))
k <- matrix(0, nrow = n, ncol = n)
for (j in 1:n) {
for (i in j:n) {
k[i, j] <- (j - 1) * n + i - 0.5 * j * (j - 1)
}
}
return(list(k = k, I = I))
}
vec <- function (x) {
if (!is.matrix(x)) {
stop("argument x is not a matrix")
}
if (!is.numeric(x)) {
stop("argument x is not a numeric matrix")
}
return(t(t(as.vector(x))))
}
E.matrices <- function (n) {
if (missing(n))
stop("argument n is missing")
if (!is.numeric(n))
stop("argument n is not numeric")
if (n != trunc(n))
stop("argument n is not an integer")
if (n < 2)
stop("argument n is less than 2")
I <- diag(rep(1, n))
E <- list()
for (i in 1:n) {
E[[i]] <- list()
for (j in 1:n) {
E[[i]][[j]] <- I[i, ] %o% I[j, ]
}
}
return(E)
}
|
