In stla/matrixsampling: Simulations of Matrix Variate Distributions
knitr::opts_chunk$set(collapse = TRUE, fig.width = 4, fig.height = 4, cache = TRUE)
library(matrixsampling)
Wishart central
phiW <- function(nu, Sigma, z = seq(0.01, 4, length.out = 20)){
p <- nrow(Sigma)
sapply(z,
function(z) complexplus::Det(diag(p) - 2*1i*(z*diag(p)+matrix(z,p,p))%*%Sigma)^(-nu/2))
}
tr <- function(x) sum(diag(x))
Phi <- function(sims, z = seq(0.01, 4, length.out = 20)){
p <- nrow(sims[,,1])
sapply(z,
function(z) mean(apply(sims, 3,
function(x) exp(1i*tr((z*diag(p)+matrix(z,p,p))%*%x)))))
}
nsims <- 30000
p <- 3
Sigma <- toeplitz(p:1)/100
# nu > p-1
nu <- p-1 + 3
W <- rwishart(nsims, nu, Sigma)
phisims <- Phi(W)
layout(t(1:2))
plot(Re(phiW(nu, Sigma)), type="o", pch=19)
lines(Re(phisims), type="o", col="red", pch=19)
plot(Im(phiW(nu, Sigma)), type="o", pch=19)
lines(Im(phisims), type="o", col="red", pch=19)
# nu integer <p
nu <- p-2
W <- rwishart(nsims, nu, Sigma)
phisims <- Phi(W)
layout(t(1:2))
plot(Re(phiW(nu, Sigma)), type="o", pch=19)
lines(Re(phisims), type="o", col="red", pch=19)
plot(Im(phiW(nu, Sigma)), type="o", pch=19)
lines(Im(phisims), type="o", col="red", pch=19)
# check singular case
Sigma <- cbind(c(3,2,0),c(2,2,0), c(0,0,0))
r <- 2 # rank Sigma
SigmaPlus <- MASS::ginv(Sigma)
nu <- 2
sims <- rwishart(nsims, nu, Sigma)
sims2 <- numeric(nsims)
for(i in 1:nsims){
y <- rnorm(3)
sims2[i] <- (t(y) %*% SigmaPlus %*% y) /
(t(y) %*% MASS::ginv(sims[,,i]) %*% y)
}
curve(ecdf(sims2)(x), to=10)
curve(pchisq(x, nu-r+1), add=TRUE, col="red")
# check tcrossprod(matrix normal)
Sigma <- toeplitz(2:1)
nsims <- 50000
sims_normal <- rmatrixnormal(nsims, M = matrix(0, 2, 3), U = Sigma, V = diag(3))
sims_wishart <- rwishart(nsims, nu = 3, Sigma = Sigma)
sims <- array(NA_real_, dim = c(2,2,nsims))
for(i in 1:nsims){
sims[,,i] <- tcrossprod(sims_normal[,,i])
}
phisims <- Phi(sims)
nu <- 3
layout(t(1:2))
plot(Re(phiW(nu, Sigma)), type="o", pch=19)
lines(Re(phisims), type="o", col="red", pch=19)
plot(Im(phiW(nu, Sigma)), type="o", pch=19)
lines(Im(phisims), type="o", col="red", pch=19)
Wishart noncentral
phiW <- function(nu, Sigma, Theta, z = seq(0.01, 4, length.out = 20)){
p <- nrow(Sigma)
sapply(z,
function(z){
complexplus::Det(diag(p) - 2*1i*(z*diag(p)+matrix(z,p,p))%*%Sigma)^(-nu/2) *
exp(1i*tr(solve(diag(p) - 2*1i*(z*diag(p)+matrix(z,p,p))%*%Sigma) %*% (z*diag(p)+matrix(z,p,p)) %*% Theta))
})
}
nsims <- 30000
p <- 3
Sigma <- toeplitz(p:1)/100
Theta <- matrix(0.1, p, p)
# nu > 2*p-1
nu <- 2*p-1 + 3
W <- rwishart(nsims, nu, Sigma, Theta)
phisims <- Phi(W)
layout(t(1:2))
plot(Re(phiW(nu, Sigma, Theta)), type="o", pch=19)
lines(Re(phisims), type="o", col="red", pch=19)
plot(Im(phiW(nu, Sigma, Theta)), type="o", pch=19)
lines(Im(phisims), type="o", col="red", pch=19)
# nu <= 2*p-1 integer
nu <- p+1
W <- rwishart(nsims, nu, Sigma, Theta)
phisims <- Phi(W)
layout(t(1:2))
plot(Re(phiW(nu, Sigma, Theta)), type="o", pch=19)
lines(Re(phisims), type="o", col="red", pch=19)
plot(Im(phiW(nu, Sigma, Theta)), type="o", pch=19)
lines(Im(phisims), type="o", col="red", pch=19)
# nu < 2*p-1 noninteger
nu <- p+0.5
W <- rwishart(nsims, nu, Sigma, Theta)
phisims <- Phi(W)
layout(t(1:2))
plot(Re(phiW(nu, Sigma, Theta)), type="o", pch=19)
lines(Re(phisims), type="o", col="red", pch=19)
plot(Im(phiW(nu, Sigma, Theta)), type="o", pch=19)
lines(Im(phisims), type="o", col="red", pch=19)
# nu = p-1
nu <- p-1
W <- rwishart(nsims, nu, Sigma, Theta)
phisims <- Phi(W)
layout(t(1:2))
plot(Re(phiW(nu, Sigma, Theta)), type="o", pch=19)
lines(Re(phisims), type="o", col="red", pch=19)
plot(Im(phiW(nu, Sigma, Theta)), type="o", pch=19)
lines(Im(phisims), type="o", col="red", pch=19)
Beta central
rU <- function(nsims, p, e, h, theta=0){
betasims <- rbeta(nsims, e/2, h/2, ncp=theta)
for(i in seq_len(p-1)){
betasims <- betasims * rbeta(nsims, (e-i)/2, h/2)
}
betasims
}
rUprime <- function(nsims, p, e, h){
betasims <- brr::rbeta2(nsims, e/2, h/2, scale=1)
for(i in seq_len(p-1)){
betasims <- betasims * brr::rbeta2(nsims, (e-i)/2, (h-i)/2, scale=1)
}
betasims
}
nsims <- 40000
# Beta I
p <- 3
a <- 2; b <- 2
detsims <- apply(rmatrixbeta(nsims, p, a, b), 3, det)
curve(ecdf(detsims)(x), to=quantile(detsims, 0.8))
betasims <- rU(nsims, p, 2*a, 2*b)
curve(ecdf(betasims)(x), add=TRUE, col="red")
# Beta II
detsims <- apply(rmatrixbetaII(nsims, p, a, b), 3, det)
curve(ecdf(detsims)(x), to=quantile(detsims, 0.8))
betasims <- rUprime(nsims, p, 2*a, 2*b)
curve(ecdf(betasims)(x), add=TRUE, col="red")
Beta I noncentral rank one
theta <- 2
detsims <- apply(rmatrixbeta(nsims, p, a, b, Theta1 = diag(c(theta, rep(0,p-1)))), 3, det)
curve(ecdf(detsims)(x), to=quantile(detsims, 0.8))
betasims <- rU(nsims, p, 2*a, 2*b, theta)
curve(ecdf(betasims)(x), add=TRUE, col="red")
# check p=1
Theta <- 1
sims <- rmatrixbeta(nsims, 1, a, b, Theta1=Theta)[1,1,]
curve(ecdf(sims)(x))
curve(pbeta(x, a, b, ncp=Theta), add=TRUE, col="red")
#
sims <- rmatrixbeta(nsims, 1, a, b, Theta2=Theta)[1,1,]
curve(ecdf(sims)(x))
curve(1-pbeta(1-x, b, a, ncp=Theta), add=TRUE, col="blue")
Beta - the two definitions
Phi <- function(sims, z = seq(0.01, 4, length.out = 20)){
sapply(z,
function(z) mean(apply(sims, 3,
function(x) exp(1i*tr((z*diag(p)+matrix(z,p,p))%*%x)))))
}
VtoU <- function(Vsims){
array(apply(Vsims, 3,
function(V){
chol2inv(chol(diag(p)+V))%*%V
}), dim=dim(Vsims))
}
p <- 5
a <- 3; b <- 3
Theta1 <- diag(p)
Theta2 <- 3*diag(p)
U1 <- rmatrixbeta(nsims, p, a, b, def=1, Theta1=Theta1, Theta2=Theta2)
V1 <- rmatrixbetaII(nsims, p, a, b, def=1,Theta1=Theta1, Theta2=Theta2)
U2 <- rmatrixbeta(nsims, p, a, b, def=2, Theta1=Theta1, Theta2=Theta2)
V2 <- rmatrixbetaII(nsims, p, a, b, def=2, Theta1=Theta1, Theta2=Theta2)
phiU1 <- Phi(U1)
phiU2 <- Phi(U2)
phiV1 <- Phi(VtoU(V1))
phiV2 <- Phi(VtoU(V2))
layout(t(1:2))
plot(Re(phiU1), type="o", pch=19)
lines(Re(phiU2), col="red", type="o", pch=19)
lines(Re(phiV1), col="green", type="o", pch=19)
lines(Re(phiV2), col="blue", type="o", pch=19)
plot(Im(phiU1), type="o", pch=19)
lines(Im(phiU2), col="red", type="o", pch=19)
lines(Im(phiV1), col="green", type="o", pch=19)
lines(Im(phiV2), col="blue", type="o", pch=19)
# def 2
Theta1 <- diag(p) + matrix(1, p, p)
Theta2 <- toeplitz(p:1)
U2 <- rmatrixbeta(nsims, p, a, b, def=2, Theta1=Theta1, Theta2=Theta2)
V2 <- rmatrixbetaII(nsims, p, a, b, def=2, Theta1=Theta1, Theta2=Theta2)
phiU2 <- Phi(U2)
phiV2 <- Phi(VtoU(V2))
layout(t(c(1,2)))
plot(Re(phiU2), type="o", pch=19)
lines(Re(phiV2), col="blue", type="o", pch=19)
plot(Im(phiU2), type="o", pch=19)
lines(Im(phiV2), col="blue", type="o", pch=19)
Student
nsims <- 30000
nu <- 6
p <- 2; m <- 3
M <- matrix(0, p, m)
U <- toeplitz(p:1)
V <- toeplitz(m:1)
Vinv <- solve(V)
Uroot <- matrixsampling:::matrixroot(U)
Urootinv <- solve(Uroot)
Tsims <- rmatrixt(nsims, nu, M, U, V)
Bsims <- rmatrixbeta(nsims, p, (nu + p - 1)/2, m/2)
BIIsims <- rmatrixbetaII(nsims, p, m/2, (nu + p - 1)/2)
sims1 <- sims2 <- array(NA_real_, dim=c(p,p,nsims))
for(i in 1:nsims){
sims1[,,i] <- Uroot %*% solve(U + Tsims[,,i] %*% Vinv %*% t(Tsims[,,i])) %*% Uroot
sims2[,,i] <- Urootinv %*% Tsims[,,i] %*% Vinv %*% t(Tsims[,,i]) %*% Urootinv
}
layout(t(1:2))
curve(ecdf(apply(Bsims, 3, sum))(x), to=2)
curve(ecdf(apply(sims1, 3, sum))(x), add=TRUE, col="red")
curve(ecdf(apply(BIIsims, 3, sum))(x), to=2)
curve(ecdf(apply(sims2, 3, sum))(x), add=TRUE, col="red")
phiB <- Phi(Bsims)
phi1 <- Phi(sims1)
layout(t(c(1,2)))
plot(Re(phiB), type="o", pch=19)
lines(Re(phi1), col="blue", type="o", pch=19)
plot(Im(phiB), type="o", pch=19)
lines(Im(phi1), col="blue", type="o", pch=19)
phiBII <- Phi(VtoU(BIIsims))
phi2 <- Phi(VtoU(sims2))
layout(t(c(1,2)))
plot(Re(phiBII), type="o", pch=19)
lines(Re(phi2), col="blue", type="o", pch=19)
plot(Im(phiBII), type="o", pch=19)
lines(Im(phi2), col="blue", type="o", pch=19)
# Beta II by conditional approach
nsims <- 40000
p <- 6
a <- (p-1)/2+4; b <- (p-1)/2+4
p2 <- 4; p1 <- p-p2
V11 <- rmatrixbetaII(nsims, p1, a, b-p2/2)
V221 <- rmatrixbetaII(nsims, p2, a-p1/2, b)
V <- array(NA_real_, dim=c(p,p,nsims))
V22sims <- array(NA_real_, dim=c(p2,p2,nsims))
for(i in 1:nsims){
V21 <-
matrix(rmatrixt(1, 2*a+2*b-p+1, matrix(0, p2, p1),
diag(p2)+V221[,,i],
V11[,,i]*(diag(p1)+V11[,,i]), checkSymmetry = FALSE)[,,1], p2, p1)
V22 <- V22sims[,,i] <- V221[,,i] + V21%*%chol2inv(chol(V11[,,i]))%*%t(V21)
V[,,i] <- cbind(rbind(V11[,,i], V21), rbind(t(V21), V22))
}
detsims <- apply(V, 3, det)
betasims <- rUprime(nsims, p, 2*a, 2*b)
curve(ecdf(detsims)(x), from=quantile(detsims, 0.3), to=quantile(detsims, 0.7))
curve(ecdf(betasims)(x), add=TRUE, col="red")
#
detsims <- apply(V, 3, det)/apply(sweep(V, 1:2, diag(p), "+"), 3, det)
#detsims <- apply(array(apply(V, 3, function(x) x%*%solve(diag(p)+x)), dim=c(p,p,nsims)), 3, det)
betasims <- rU(nsims, p, 2*a, 2*b)
curve(ecdf(detsims)(x), from=quantile(detsims, 0.3), to=quantile(detsims, 0.7))
curve(ecdf(betasims)(x), add=TRUE, col="red") # ??????????
#
phiV1 <- Phi(VtoU(V))
phiV2 <- Phi(VtoU(rmatrixbetaII(nsims, p, a, b)))
layout(t(1:2))
plot(Re(phiV1), type="o", pch=19)
lines(Re(phiV2), col="blue", type="o", pch=19)
plot(Im(phiV1), type="o", pch=19)
lines(Im(phiV2), col="blue", type="o", pch=19)
V <- rmatrixbetaII(nsims, p, a, b)
V11 <- V[1:p1, 1:p1, ]
detsims <- apply(V11, 3, det)
betasims <- rUprime(nsims, p1, 2*a, 2*b-p2)
curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9))
curve(ecdf(betasims)(x), add=TRUE, col="red")
V221 <- array(NA_real_, dim=c(p2,p2,nsims))
for(i in 1:nsims){
V221[,,i] <- V[(p1+1):p,(p1+1):p,i] - V[(p1+1):p,1:p1,i] %*% chol2inv(chol(V[1:p1,1:p1,i])) %*% t(V[(p1+1):p,1:p1,i])
}
detsims <- apply(V221, 3, det)
betasims <- rUprime(nsims, p2, 2*a-p1, 2*b)
curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9))
curve(ecdf(betasims)(x), add=TRUE, col="red")
#
detsims <- apply(V11, 3, det)/apply(sweep(V11, 1:2, diag(p1), "+"), 3, det)
betasims <- rU(nsims, p1, 2*a, 2*b-p2)
curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9))
curve(ecdf(betasims)(x), add=TRUE, col="red")
#
detsims <- apply(V221, 3, det)/apply(sweep(V221, 1:2, diag(p2), "+"), 3, det)
betasims <- rU(nsims, p2, 2*a-p1, 2*b)
curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9))
curve(ecdf(betasims)(x), add=TRUE, col="red")
#
cor(V11[1,1,], V221[1,1,] + V221[2,1,])
cor(V11[1,1,], V221[1,1,]*V221[2,2,])
cor(V11[1,2,], V221[2,1,]^2+V221[2,2,])
Inverted t
# thm 5.2.3 Gupta & Nagar
nsims <- 40000
p <- 3; n1 <- 5; n2 <- 4
Sigma <- toeplitz(p:1)/100
X <- rmatrixnormal(nsims, matrix(0, p, n1), Sigma, diag(n1))
S2 <- rwishart(nsims, n2, Sigma)
Z <- array(NA_real_, dim=c(p,n1,nsims))
for(i in 1:nsims){
Z[,,i] <- matrixsampling:::invsqrtm(S2[,,i]+tcrossprod(X[,,i]))%*%X[,,i]
}
IT <- rmatrixit(nsims, n2-p+1, matrix(0, p, n1), diag(p), diag(n1))
Phi <- function(sims, z = seq(0.01, 4, length.out = 20)){
sapply(z,
function(z) mean(apply(sims, 3,
function(x) exp(-1i*sum(rep(z, p*n1)*c(t(x)))))))
}
PhiZ <- Phi(Z)
PhiIT <- Phi(IT)
layout(t(1:2))
plot(Re(PhiZ), type="o", pch=19)
lines(Re(PhiIT), type="o", pch=19, col="red")
plot(Im(PhiZ), type="o", pch=19) ## bizarre !
lines(Im(PhiIT), type="o", pch=19, col="red")
round(apply(Z, 1:2, var), 2)
round(apply(IT, 1:2, var), 2)
curve(ecdf(Z[1,1,]+Z[1,2,]+Z[3,3,]+Z[3,5,]^2)(x), from=-1.5, to=1.5)
curve(ecdf(IT[1,1,]+IT[1,2,]+IT[3,3,]+IT[3,5,]^2)(x), add=TRUE, col="red")
Z2 <- IT2 <- array(NA_real_, dim=c(p,p,nsims))
for(i in 1:nsims){
Z2[,,i] <- tcrossprod(Z[,,i])
IT2[,,i] <- tcrossprod(IT[,,i])
}
Phi <- function(sims, z = seq(0.01, 4, length.out = 20)){
sapply(z,
function(z) mean(apply(sims, 3,
function(x) exp(-1i*sum(rep(z, p*p)*c(t(x)))))))
}
PhiZ2 <- Phi(Z2)
PhiIT2 <- Phi(IT2)
layout(t(1:2))
plot(Re(PhiZ2), type="o", pch=19)
lines(Re(PhiIT2), type="o", pch=19, col="red")
plot(Im(PhiZ2), type="o", pch=19)
lines(Im(PhiIT2), type="o", pch=19, col="red")
# Beta I by conditional approach
nsims <- 40000
p <- 6
a <- (p-1)/2+6; b <- (p-1)/2+6
p2 <- 4; p1 <- p-p2
U11 <- rmatrixbeta(nsims, p1, a, b)
U221 <- rmatrixbeta(nsims, p2, a-p1/2, b)
U <- array(NA_real_, dim=c(p,p,nsims))
for(i in 1:nsims){
U21 <-
matrix(rmatrixit(1, 2*b-p+1, matrix(0, p2, p1),
diag(p2)-U221[,,i],
U11[,,i]*(diag(p1)-U11[,,i]), checkSymmetry = FALSE)[,,1], p2, p1)
U22 <- U221[,,i] + U21%*%chol2inv(chol(U11[,,i]))%*%t(U21)
U[,,i] <- cbind(rbind(U11[,,i], U21), rbind(t(U21), U22))
}
detsims <- apply(U, 3, det)
betasims <- rU(nsims, p, 2*a, 2*b)
curve(ecdf(detsims)(x), from=quantile(detsims, 0.3), to=quantile(detsims, 0.7))
curve(ecdf(betasims)(x), add=TRUE, col="red")
#
detsims <- apply(U, 3, det)/apply(-sweep(U, 1:2, diag(p), "-"), 3, det)
betasims <- rUprime(nsims, p, 2*a, 2*b)
curve(ecdf(detsims)(x), from=quantile(detsims, 0.3), to=quantile(detsims, 0.7))
curve(ecdf(betasims)(x), add=TRUE, col="red") # ??????????
phiU1 <- Phi(U, z=seq(1, 10, by=1))
phiU2 <- Phi(rmatrixbeta(nsims, p, a, b), z=seq(1, 10, by=1))
layout(t(1:2))
plot(Re(phiU1), type="o", pch=19)
lines(Re(phiU2), col="blue", type="o", pch=19)
plot(Im(phiU1), type="o", pch=19)
lines(Im(phiU2), col="blue", type="o", pch=19)
U <- rmatrixbeta(nsims, p, a, b)
U11 <- U[1:p1, 1:p1, ]
detsims <- apply(U11, 3, det)
betasims <- rU(nsims, p1, 2*a, 2*b)
curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9))
curve(ecdf(betasims)(x), add=TRUE, col="red")
U221 <- array(NA_real_, dim=c(p2,p2,nsims))
for(i in 1:nsims){
U221[,,i] <- U[(p1+1):p,(p1+1):p,i] - U[(p1+1):p,1:p1,i] %*% chol2inv(chol(U[1:p1,1:p1,i])) %*% t(U[(p1+1):p,1:p1,i])
}
detsims <- apply(U221, 3, det)
betasims <- rU(nsims, p2, 2*a-p1, 2*b)
curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9))
curve(ecdf(betasims)(x), add=TRUE, col="red")
#
detsims <- apply(U[1:p1, 1:p1, ], 3, det)/apply(-sweep(U[1:p1, 1:p1, ], 1:2, diag(p1), "-"), 3, det)
betasims <- rUprime(nsims, p1, 2*a, 2*b)
curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9))
curve(ecdf(betasims)(x), add=TRUE, col="red")
#
detsims <- apply(U221, 3, det)/apply(-sweep(U221, 1:2, diag(p2), "-"), 3, det)
betasims <- rUprime(nsims, p2, 2*a-p1, 2*b)
curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9))
curve(ecdf(betasims)(x), add=TRUE, col="red")
#
cor(U11[1,1,], U221[2,1,])
cor(U11[1,1,], U221[2,2,])
cor(U11[1,2,], U221[2,1,])
Confluent kind one
nsims <- 50000
nu1 <- 5; nu2 <- 8; theta <- 4; p <- 2
X1 <- rmatrixgamma(nsims, nu1, theta, p = p)
X2 <- rmatrixgamma(nsims, nu2, theta, p = p)
sims <- array(NA_real_, dim = c(p,p,nsims))
for(i in 1:nsims){
sims[,,i] <- (X1[,,i]+X2[,,i]) %*% chol2inv(chol(X1[,,i])) %*% (X1[,,i]+X2[,,i])
}
phisims <- Phi(sims, z=seq(1, 10, by=1)*500000)
CHsims <- rmatrixCHkind1(nsims, nu1+nu2, 2*nu1+nu2, 2*nu1+2*nu2, theta=theta, p=p)
phiCH <-
Phi(CHsims,
z=seq(1, 10, by=1)*500000)
layout(t(1:2))
plot(Re(phisims), type="o", pch=19)
lines(Re(phiCH), col="blue", type="o", pch=19)
plot(Im(phisims), type="o", pch=19)
lines(Im(phiCH), col="blue", type="o", pch=19)
mean(apply(sims, 3, det))
mean(apply(CHsims, 3, det))
stla/matrixsampling documentation built on Aug. 27, 2019, 3:03 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(collapse = TRUE, fig.width = 4, fig.height = 4, cache = TRUE) library(matrixsampling)
Wishart central
phiW <- function(nu, Sigma, z = seq(0.01, 4, length.out = 20)){ p <- nrow(Sigma) sapply(z, function(z) complexplus::Det(diag(p) - 2*1i*(z*diag(p)+matrix(z,p,p))%*%Sigma)^(-nu/2)) }
tr <- function(x) sum(diag(x)) Phi <- function(sims, z = seq(0.01, 4, length.out = 20)){ p <- nrow(sims[,,1]) sapply(z, function(z) mean(apply(sims, 3, function(x) exp(1i*tr((z*diag(p)+matrix(z,p,p))%*%x))))) }
nsims <- 30000 p <- 3 Sigma <- toeplitz(p:1)/100 # nu > p-1 nu <- p-1 + 3 W <- rwishart(nsims, nu, Sigma) phisims <- Phi(W) layout(t(1:2)) plot(Re(phiW(nu, Sigma)), type="o", pch=19) lines(Re(phisims), type="o", col="red", pch=19) plot(Im(phiW(nu, Sigma)), type="o", pch=19) lines(Im(phisims), type="o", col="red", pch=19)
# nu integer <p nu <- p-2 W <- rwishart(nsims, nu, Sigma) phisims <- Phi(W) layout(t(1:2)) plot(Re(phiW(nu, Sigma)), type="o", pch=19) lines(Re(phisims), type="o", col="red", pch=19) plot(Im(phiW(nu, Sigma)), type="o", pch=19) lines(Im(phisims), type="o", col="red", pch=19)
# check singular case Sigma <- cbind(c(3,2,0),c(2,2,0), c(0,0,0)) r <- 2 # rank Sigma SigmaPlus <- MASS::ginv(Sigma) nu <- 2 sims <- rwishart(nsims, nu, Sigma) sims2 <- numeric(nsims) for(i in 1:nsims){ y <- rnorm(3) sims2[i] <- (t(y) %*% SigmaPlus %*% y) / (t(y) %*% MASS::ginv(sims[,,i]) %*% y) } curve(ecdf(sims2)(x), to=10) curve(pchisq(x, nu-r+1), add=TRUE, col="red")
# check tcrossprod(matrix normal) Sigma <- toeplitz(2:1) nsims <- 50000 sims_normal <- rmatrixnormal(nsims, M = matrix(0, 2, 3), U = Sigma, V = diag(3)) sims_wishart <- rwishart(nsims, nu = 3, Sigma = Sigma) sims <- array(NA_real_, dim = c(2,2,nsims)) for(i in 1:nsims){ sims[,,i] <- tcrossprod(sims_normal[,,i]) } phisims <- Phi(sims) nu <- 3 layout(t(1:2)) plot(Re(phiW(nu, Sigma)), type="o", pch=19) lines(Re(phisims), type="o", col="red", pch=19) plot(Im(phiW(nu, Sigma)), type="o", pch=19) lines(Im(phisims), type="o", col="red", pch=19)
Wishart noncentral
phiW <- function(nu, Sigma, Theta, z = seq(0.01, 4, length.out = 20)){ p <- nrow(Sigma) sapply(z, function(z){ complexplus::Det(diag(p) - 2*1i*(z*diag(p)+matrix(z,p,p))%*%Sigma)^(-nu/2) * exp(1i*tr(solve(diag(p) - 2*1i*(z*diag(p)+matrix(z,p,p))%*%Sigma) %*% (z*diag(p)+matrix(z,p,p)) %*% Theta)) }) }
nsims <- 30000 p <- 3 Sigma <- toeplitz(p:1)/100 Theta <- matrix(0.1, p, p) # nu > 2*p-1 nu <- 2*p-1 + 3 W <- rwishart(nsims, nu, Sigma, Theta) phisims <- Phi(W) layout(t(1:2)) plot(Re(phiW(nu, Sigma, Theta)), type="o", pch=19) lines(Re(phisims), type="o", col="red", pch=19) plot(Im(phiW(nu, Sigma, Theta)), type="o", pch=19) lines(Im(phisims), type="o", col="red", pch=19)
# nu <= 2*p-1 integer nu <- p+1 W <- rwishart(nsims, nu, Sigma, Theta) phisims <- Phi(W) layout(t(1:2)) plot(Re(phiW(nu, Sigma, Theta)), type="o", pch=19) lines(Re(phisims), type="o", col="red", pch=19) plot(Im(phiW(nu, Sigma, Theta)), type="o", pch=19) lines(Im(phisims), type="o", col="red", pch=19)
# nu < 2*p-1 noninteger nu <- p+0.5 W <- rwishart(nsims, nu, Sigma, Theta) phisims <- Phi(W) layout(t(1:2)) plot(Re(phiW(nu, Sigma, Theta)), type="o", pch=19) lines(Re(phisims), type="o", col="red", pch=19) plot(Im(phiW(nu, Sigma, Theta)), type="o", pch=19) lines(Im(phisims), type="o", col="red", pch=19)
# nu = p-1 nu <- p-1 W <- rwishart(nsims, nu, Sigma, Theta) phisims <- Phi(W) layout(t(1:2)) plot(Re(phiW(nu, Sigma, Theta)), type="o", pch=19) lines(Re(phisims), type="o", col="red", pch=19) plot(Im(phiW(nu, Sigma, Theta)), type="o", pch=19) lines(Im(phisims), type="o", col="red", pch=19)
Beta central
rU <- function(nsims, p, e, h, theta=0){ betasims <- rbeta(nsims, e/2, h/2, ncp=theta) for(i in seq_len(p-1)){ betasims <- betasims * rbeta(nsims, (e-i)/2, h/2) } betasims } rUprime <- function(nsims, p, e, h){ betasims <- brr::rbeta2(nsims, e/2, h/2, scale=1) for(i in seq_len(p-1)){ betasims <- betasims * brr::rbeta2(nsims, (e-i)/2, (h-i)/2, scale=1) } betasims }
nsims <- 40000 # Beta I p <- 3 a <- 2; b <- 2 detsims <- apply(rmatrixbeta(nsims, p, a, b), 3, det) curve(ecdf(detsims)(x), to=quantile(detsims, 0.8)) betasims <- rU(nsims, p, 2*a, 2*b) curve(ecdf(betasims)(x), add=TRUE, col="red") # Beta II detsims <- apply(rmatrixbetaII(nsims, p, a, b), 3, det) curve(ecdf(detsims)(x), to=quantile(detsims, 0.8)) betasims <- rUprime(nsims, p, 2*a, 2*b) curve(ecdf(betasims)(x), add=TRUE, col="red")
Beta I noncentral rank one
theta <- 2 detsims <- apply(rmatrixbeta(nsims, p, a, b, Theta1 = diag(c(theta, rep(0,p-1)))), 3, det) curve(ecdf(detsims)(x), to=quantile(detsims, 0.8)) betasims <- rU(nsims, p, 2*a, 2*b, theta) curve(ecdf(betasims)(x), add=TRUE, col="red")
# check p=1 Theta <- 1 sims <- rmatrixbeta(nsims, 1, a, b, Theta1=Theta)[1,1,] curve(ecdf(sims)(x)) curve(pbeta(x, a, b, ncp=Theta), add=TRUE, col="red") # sims <- rmatrixbeta(nsims, 1, a, b, Theta2=Theta)[1,1,] curve(ecdf(sims)(x)) curve(1-pbeta(1-x, b, a, ncp=Theta), add=TRUE, col="blue")
Beta - the two definitions
Phi <- function(sims, z = seq(0.01, 4, length.out = 20)){ sapply(z, function(z) mean(apply(sims, 3, function(x) exp(1i*tr((z*diag(p)+matrix(z,p,p))%*%x))))) } VtoU <- function(Vsims){ array(apply(Vsims, 3, function(V){ chol2inv(chol(diag(p)+V))%*%V }), dim=dim(Vsims)) } p <- 5 a <- 3; b <- 3 Theta1 <- diag(p) Theta2 <- 3*diag(p) U1 <- rmatrixbeta(nsims, p, a, b, def=1, Theta1=Theta1, Theta2=Theta2) V1 <- rmatrixbetaII(nsims, p, a, b, def=1,Theta1=Theta1, Theta2=Theta2) U2 <- rmatrixbeta(nsims, p, a, b, def=2, Theta1=Theta1, Theta2=Theta2) V2 <- rmatrixbetaII(nsims, p, a, b, def=2, Theta1=Theta1, Theta2=Theta2) phiU1 <- Phi(U1) phiU2 <- Phi(U2) phiV1 <- Phi(VtoU(V1)) phiV2 <- Phi(VtoU(V2)) layout(t(1:2)) plot(Re(phiU1), type="o", pch=19) lines(Re(phiU2), col="red", type="o", pch=19) lines(Re(phiV1), col="green", type="o", pch=19) lines(Re(phiV2), col="blue", type="o", pch=19) plot(Im(phiU1), type="o", pch=19) lines(Im(phiU2), col="red", type="o", pch=19) lines(Im(phiV1), col="green", type="o", pch=19) lines(Im(phiV2), col="blue", type="o", pch=19)
# def 2 Theta1 <- diag(p) + matrix(1, p, p) Theta2 <- toeplitz(p:1) U2 <- rmatrixbeta(nsims, p, a, b, def=2, Theta1=Theta1, Theta2=Theta2) V2 <- rmatrixbetaII(nsims, p, a, b, def=2, Theta1=Theta1, Theta2=Theta2) phiU2 <- Phi(U2) phiV2 <- Phi(VtoU(V2)) layout(t(c(1,2))) plot(Re(phiU2), type="o", pch=19) lines(Re(phiV2), col="blue", type="o", pch=19) plot(Im(phiU2), type="o", pch=19) lines(Im(phiV2), col="blue", type="o", pch=19)
Student
nsims <- 30000 nu <- 6 p <- 2; m <- 3 M <- matrix(0, p, m) U <- toeplitz(p:1) V <- toeplitz(m:1) Vinv <- solve(V) Uroot <- matrixsampling:::matrixroot(U) Urootinv <- solve(Uroot) Tsims <- rmatrixt(nsims, nu, M, U, V) Bsims <- rmatrixbeta(nsims, p, (nu + p - 1)/2, m/2) BIIsims <- rmatrixbetaII(nsims, p, m/2, (nu + p - 1)/2) sims1 <- sims2 <- array(NA_real_, dim=c(p,p,nsims)) for(i in 1:nsims){ sims1[,,i] <- Uroot %*% solve(U + Tsims[,,i] %*% Vinv %*% t(Tsims[,,i])) %*% Uroot sims2[,,i] <- Urootinv %*% Tsims[,,i] %*% Vinv %*% t(Tsims[,,i]) %*% Urootinv } layout(t(1:2)) curve(ecdf(apply(Bsims, 3, sum))(x), to=2) curve(ecdf(apply(sims1, 3, sum))(x), add=TRUE, col="red") curve(ecdf(apply(BIIsims, 3, sum))(x), to=2) curve(ecdf(apply(sims2, 3, sum))(x), add=TRUE, col="red")
phiB <- Phi(Bsims) phi1 <- Phi(sims1) layout(t(c(1,2))) plot(Re(phiB), type="o", pch=19) lines(Re(phi1), col="blue", type="o", pch=19) plot(Im(phiB), type="o", pch=19) lines(Im(phi1), col="blue", type="o", pch=19)
phiBII <- Phi(VtoU(BIIsims)) phi2 <- Phi(VtoU(sims2)) layout(t(c(1,2))) plot(Re(phiBII), type="o", pch=19) lines(Re(phi2), col="blue", type="o", pch=19) plot(Im(phiBII), type="o", pch=19) lines(Im(phi2), col="blue", type="o", pch=19)
# Beta II by conditional approach nsims <- 40000 p <- 6 a <- (p-1)/2+4; b <- (p-1)/2+4 p2 <- 4; p1 <- p-p2 V11 <- rmatrixbetaII(nsims, p1, a, b-p2/2) V221 <- rmatrixbetaII(nsims, p2, a-p1/2, b) V <- array(NA_real_, dim=c(p,p,nsims)) V22sims <- array(NA_real_, dim=c(p2,p2,nsims)) for(i in 1:nsims){ V21 <- matrix(rmatrixt(1, 2*a+2*b-p+1, matrix(0, p2, p1), diag(p2)+V221[,,i], V11[,,i]*(diag(p1)+V11[,,i]), checkSymmetry = FALSE)[,,1], p2, p1) V22 <- V22sims[,,i] <- V221[,,i] + V21%*%chol2inv(chol(V11[,,i]))%*%t(V21) V[,,i] <- cbind(rbind(V11[,,i], V21), rbind(t(V21), V22)) } detsims <- apply(V, 3, det) betasims <- rUprime(nsims, p, 2*a, 2*b) curve(ecdf(detsims)(x), from=quantile(detsims, 0.3), to=quantile(detsims, 0.7)) curve(ecdf(betasims)(x), add=TRUE, col="red") # detsims <- apply(V, 3, det)/apply(sweep(V, 1:2, diag(p), "+"), 3, det) #detsims <- apply(array(apply(V, 3, function(x) x%*%solve(diag(p)+x)), dim=c(p,p,nsims)), 3, det) betasims <- rU(nsims, p, 2*a, 2*b) curve(ecdf(detsims)(x), from=quantile(detsims, 0.3), to=quantile(detsims, 0.7)) curve(ecdf(betasims)(x), add=TRUE, col="red") # ?????????? # phiV1 <- Phi(VtoU(V)) phiV2 <- Phi(VtoU(rmatrixbetaII(nsims, p, a, b))) layout(t(1:2)) plot(Re(phiV1), type="o", pch=19) lines(Re(phiV2), col="blue", type="o", pch=19) plot(Im(phiV1), type="o", pch=19) lines(Im(phiV2), col="blue", type="o", pch=19)
V <- rmatrixbetaII(nsims, p, a, b) V11 <- V[1:p1, 1:p1, ] detsims <- apply(V11, 3, det) betasims <- rUprime(nsims, p1, 2*a, 2*b-p2) curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9)) curve(ecdf(betasims)(x), add=TRUE, col="red") V221 <- array(NA_real_, dim=c(p2,p2,nsims)) for(i in 1:nsims){ V221[,,i] <- V[(p1+1):p,(p1+1):p,i] - V[(p1+1):p,1:p1,i] %*% chol2inv(chol(V[1:p1,1:p1,i])) %*% t(V[(p1+1):p,1:p1,i]) } detsims <- apply(V221, 3, det) betasims <- rUprime(nsims, p2, 2*a-p1, 2*b) curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9)) curve(ecdf(betasims)(x), add=TRUE, col="red") # detsims <- apply(V11, 3, det)/apply(sweep(V11, 1:2, diag(p1), "+"), 3, det) betasims <- rU(nsims, p1, 2*a, 2*b-p2) curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9)) curve(ecdf(betasims)(x), add=TRUE, col="red") # detsims <- apply(V221, 3, det)/apply(sweep(V221, 1:2, diag(p2), "+"), 3, det) betasims <- rU(nsims, p2, 2*a-p1, 2*b) curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9)) curve(ecdf(betasims)(x), add=TRUE, col="red") # cor(V11[1,1,], V221[1,1,] + V221[2,1,]) cor(V11[1,1,], V221[1,1,]*V221[2,2,]) cor(V11[1,2,], V221[2,1,]^2+V221[2,2,])
Inverted t
# thm 5.2.3 Gupta & Nagar nsims <- 40000 p <- 3; n1 <- 5; n2 <- 4 Sigma <- toeplitz(p:1)/100 X <- rmatrixnormal(nsims, matrix(0, p, n1), Sigma, diag(n1)) S2 <- rwishart(nsims, n2, Sigma) Z <- array(NA_real_, dim=c(p,n1,nsims)) for(i in 1:nsims){ Z[,,i] <- matrixsampling:::invsqrtm(S2[,,i]+tcrossprod(X[,,i]))%*%X[,,i] } IT <- rmatrixit(nsims, n2-p+1, matrix(0, p, n1), diag(p), diag(n1)) Phi <- function(sims, z = seq(0.01, 4, length.out = 20)){ sapply(z, function(z) mean(apply(sims, 3, function(x) exp(-1i*sum(rep(z, p*n1)*c(t(x))))))) } PhiZ <- Phi(Z) PhiIT <- Phi(IT) layout(t(1:2)) plot(Re(PhiZ), type="o", pch=19) lines(Re(PhiIT), type="o", pch=19, col="red") plot(Im(PhiZ), type="o", pch=19) ## bizarre ! lines(Im(PhiIT), type="o", pch=19, col="red")
round(apply(Z, 1:2, var), 2) round(apply(IT, 1:2, var), 2) curve(ecdf(Z[1,1,]+Z[1,2,]+Z[3,3,]+Z[3,5,]^2)(x), from=-1.5, to=1.5) curve(ecdf(IT[1,1,]+IT[1,2,]+IT[3,3,]+IT[3,5,]^2)(x), add=TRUE, col="red")
Z2 <- IT2 <- array(NA_real_, dim=c(p,p,nsims)) for(i in 1:nsims){ Z2[,,i] <- tcrossprod(Z[,,i]) IT2[,,i] <- tcrossprod(IT[,,i]) } Phi <- function(sims, z = seq(0.01, 4, length.out = 20)){ sapply(z, function(z) mean(apply(sims, 3, function(x) exp(-1i*sum(rep(z, p*p)*c(t(x))))))) } PhiZ2 <- Phi(Z2) PhiIT2 <- Phi(IT2) layout(t(1:2)) plot(Re(PhiZ2), type="o", pch=19) lines(Re(PhiIT2), type="o", pch=19, col="red") plot(Im(PhiZ2), type="o", pch=19) lines(Im(PhiIT2), type="o", pch=19, col="red")
# Beta I by conditional approach nsims <- 40000 p <- 6 a <- (p-1)/2+6; b <- (p-1)/2+6 p2 <- 4; p1 <- p-p2 U11 <- rmatrixbeta(nsims, p1, a, b) U221 <- rmatrixbeta(nsims, p2, a-p1/2, b) U <- array(NA_real_, dim=c(p,p,nsims)) for(i in 1:nsims){ U21 <- matrix(rmatrixit(1, 2*b-p+1, matrix(0, p2, p1), diag(p2)-U221[,,i], U11[,,i]*(diag(p1)-U11[,,i]), checkSymmetry = FALSE)[,,1], p2, p1) U22 <- U221[,,i] + U21%*%chol2inv(chol(U11[,,i]))%*%t(U21) U[,,i] <- cbind(rbind(U11[,,i], U21), rbind(t(U21), U22)) } detsims <- apply(U, 3, det) betasims <- rU(nsims, p, 2*a, 2*b) curve(ecdf(detsims)(x), from=quantile(detsims, 0.3), to=quantile(detsims, 0.7)) curve(ecdf(betasims)(x), add=TRUE, col="red") # detsims <- apply(U, 3, det)/apply(-sweep(U, 1:2, diag(p), "-"), 3, det) betasims <- rUprime(nsims, p, 2*a, 2*b) curve(ecdf(detsims)(x), from=quantile(detsims, 0.3), to=quantile(detsims, 0.7)) curve(ecdf(betasims)(x), add=TRUE, col="red") # ??????????
phiU1 <- Phi(U, z=seq(1, 10, by=1)) phiU2 <- Phi(rmatrixbeta(nsims, p, a, b), z=seq(1, 10, by=1)) layout(t(1:2)) plot(Re(phiU1), type="o", pch=19) lines(Re(phiU2), col="blue", type="o", pch=19) plot(Im(phiU1), type="o", pch=19) lines(Im(phiU2), col="blue", type="o", pch=19)
U <- rmatrixbeta(nsims, p, a, b) U11 <- U[1:p1, 1:p1, ] detsims <- apply(U11, 3, det) betasims <- rU(nsims, p1, 2*a, 2*b) curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9)) curve(ecdf(betasims)(x), add=TRUE, col="red") U221 <- array(NA_real_, dim=c(p2,p2,nsims)) for(i in 1:nsims){ U221[,,i] <- U[(p1+1):p,(p1+1):p,i] - U[(p1+1):p,1:p1,i] %*% chol2inv(chol(U[1:p1,1:p1,i])) %*% t(U[(p1+1):p,1:p1,i]) } detsims <- apply(U221, 3, det) betasims <- rU(nsims, p2, 2*a-p1, 2*b) curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9)) curve(ecdf(betasims)(x), add=TRUE, col="red") # detsims <- apply(U[1:p1, 1:p1, ], 3, det)/apply(-sweep(U[1:p1, 1:p1, ], 1:2, diag(p1), "-"), 3, det) betasims <- rUprime(nsims, p1, 2*a, 2*b) curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9)) curve(ecdf(betasims)(x), add=TRUE, col="red") # detsims <- apply(U221, 3, det)/apply(-sweep(U221, 1:2, diag(p2), "-"), 3, det) betasims <- rUprime(nsims, p2, 2*a-p1, 2*b) curve(ecdf(detsims)(x), from=quantile(detsims, 0.1), to=quantile(detsims, 0.9)) curve(ecdf(betasims)(x), add=TRUE, col="red") # cor(U11[1,1,], U221[2,1,]) cor(U11[1,1,], U221[2,2,]) cor(U11[1,2,], U221[2,1,])
Confluent kind one
nsims <- 50000 nu1 <- 5; nu2 <- 8; theta <- 4; p <- 2 X1 <- rmatrixgamma(nsims, nu1, theta, p = p) X2 <- rmatrixgamma(nsims, nu2, theta, p = p) sims <- array(NA_real_, dim = c(p,p,nsims)) for(i in 1:nsims){ sims[,,i] <- (X1[,,i]+X2[,,i]) %*% chol2inv(chol(X1[,,i])) %*% (X1[,,i]+X2[,,i]) } phisims <- Phi(sims, z=seq(1, 10, by=1)*500000) CHsims <- rmatrixCHkind1(nsims, nu1+nu2, 2*nu1+nu2, 2*nu1+2*nu2, theta=theta, p=p) phiCH <- Phi(CHsims, z=seq(1, 10, by=1)*500000) layout(t(1:2)) plot(Re(phisims), type="o", pch=19) lines(Re(phiCH), col="blue", type="o", pch=19) plot(Im(phisims), type="o", pch=19) lines(Im(phiCH), col="blue", type="o", pch=19) mean(apply(sims, 3, det)) mean(apply(CHsims, 3, det))
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.