Nothing
randomGVARmodel <- function(
Nvar,
probKappaEdge = 0.1,
probKappaPositive = 0.5,
probBetaEdge = 0.1,
probBetaPositive = 0.5,
maxtry = 10,
kappaConstant = 1.1
){
try <- 0
repeat{
kappaRange = c(0.5,1)
## Approach from
# Yin, J., & Li, H. (2011). A sparse conditional gaussian graphical model for analysis of genetical genomics data. The annals of applied statistics, 5(4), 2630.
trueKappa <- matrix(0,Nvar,Nvar)
trueKappa[upper.tri(trueKappa)] <- sample(c(0,1),sum(upper.tri(trueKappa)),TRUE,prob=c(1-probKappaEdge,probKappaEdge))
trueKappa[upper.tri(trueKappa)] <- trueKappa[upper.tri(trueKappa)] * sample(c(1,-1),sum(upper.tri(trueKappa)),TRUE,prob=c(1-probKappaPositive,probKappaPositive)) *
runif(sum(upper.tri(trueKappa)), min(kappaRange),max(kappaRange))
trueKappa[lower.tri(trueKappa)] <- t(trueKappa)[lower.tri(trueKappa)]
diag(trueKappa) <- kappaConstant * rowSums(abs(trueKappa))
diag(trueKappa) <- ifelse(diag(trueKappa)==0,1,diag(trueKappa))
trueKappa <- trueKappa/diag(trueKappa)[row(trueKappa)]
trueKappa <- (trueKappa + t(trueKappa)) / 2
Vmin <- min(abs(trueKappa[trueKappa!=0]))
trueBeta <- matrix(sample(c(0,1),Nvar^2,TRUE,prob=c(1-probBetaEdge,probBetaEdge)),Nvar,Nvar)
trueBeta <- trueBeta * sample(c(-1,1),Nvar^2,TRUE,prob=c(1-probBetaPositive,probBetaPositive)) *
runif(Nvar^2, Vmin,1)
diag(trueBeta) <- Vmin
evK <- eigen(trueKappa)$values
evB <- eigen(trueBeta)$values
while(any(Re(evB)^2 + Im(evB)^2 > 1)){
trueBeta <- 0.95*trueBeta
evB <- eigen(trueBeta)$values
}
if (all(evK > 0) & all(Re(evB)^2 + Im(evB)^2 < 1)){
break
}
try <- try + 1
if (try > maxtry){
stop("Maximum number of tries reached.")
}
}
Res <- list(
kappa = trueKappa,
beta = trueBeta,
PCC = computePCC(trueKappa),
PDC = computePDC(trueBeta,trueKappa)
)
class(Res) <- "gVARmodel"
return(Res)
}
graphicalVARsim <- function(
nTime, # Number of time points
beta, # if dim is 2xnVar, assume changescores
kappa,
mean = rep(0,ncol(kappa)),
# sd = rep(1, ncol(kappa)),
init = mean,
warmup = 100,
lbound = rep(-Inf, ncol(kappa)),
ubound = rep(Inf, ncol(kappa))
# skewed = rep(0,ncol(kappa)),
# WN = FALSE
){
stopifnot(!missing(beta))
stopifnot(!missing(kappa))
Nvar <- ncol(kappa)
init <- rep(init, length = Nvar)
totTime <- nTime + warmup
Data <- t(matrix(init, Nvar, totTime))
Sigma <- solve(kappa)
# lbound <- (lbound - mean) / sd
# ubound <- (ubound - mean) / sd
#
# if (sum(skewed==rep(0,ncol(kappa)))==ncol(kappa)){
for (t in 2:totTime){
Data[t,] <- mean + t(beta %*% (Data[t-1,]-mean)) + mvtnorm::rmvnorm(1, rep(0,Nvar), Sigma)
Data[t,] <- ifelse(Data[t,] < lbound, lbound, Data[t,] )
Data[t,] <- ifelse(Data[t,] > ubound, ubound, Data[t,] )
}
# }else{
# for (t in 2:totTime){#Needed to round Omega to avoid error "not symmetrical"
# Data[t,] <- t(beta %*% Data[t-1,]) + sn::rmsn(n=1,dim=1*Nvar, mu=rep(0,Nvar),Omega=round(Sigma,digits = 5),
# alpha=skewed)[1,]
# Data[t,] <- ifelse(Data[t,] < lbound, lbound, Data[t,] )
# Data[t,] <- ifelse(Data[t,] > ubound, ubound, Data[t,] )
# }
# }
# if (WN){
# for (t in 2:totTime){
# Data[t,] <- Data[t,] + rmvnorm(1, rep(0,Nvar), Sigma)
# }
# }
return(Data[-seq_len(warmup), ,drop=FALSE])
}
#
# graphicalVARsim <- function(
# nTime, # Number of time points
# beta, # if dim is 2xnVar, assume changescores
# kappa,
# scaledMatrices = TRUE,
# mean = rep(0,ncol(kappa)),
# sd = rep(1, ncol(kappa)),
# init = mean,
# intercepts = 0,
# warmup = 100,
# pi =rep(0,ncol(kappa)),
# nudgeMean = rep(0.1,ncol(kappa)),
# nudgeSD = rep(0.1,ncol(kappa)),
# lbound = rep(-Inf, ncol(kappa)),
# ubound = rep(Inf, ncol(kappa))
# ){
#
#
# stopifnot(!missing(beta))
# stopifnot(!missing(kappa))
#
# Nvar <- ncol(kappa)
# init <- rep(init, length = Nvar)
# intercepts <- rep(intercepts, length = Nvar)
#
#
# # Temp solulution, add change scores:
# if (ncol(beta)==Nvar){
# beta <- cbind(beta,matrix(0,Nvar,Nvar))
# }
# if (ncol(beta)!=Nvar*2) stop("Beta must contain only lag-1 and changescore effects.")
#
# totTime <- nTime + warmup
#
# Data <- t(matrix(init, Nvar, totTime))
#
# Sigma <- solve(kappa)
#
# lbound <- (lbound - mean) / sd
# ubound <- (ubound - mean) / sd
#
# for (t in 3:totTime){
# residMean <- ifelse(runif(Nvar) < pi, rnorm(Nvar,nudgeMean,nudgeSD), 0)
#
# Data[t,] <- t(intercepts + beta[,1:Nvar] %*% Data[t-1,] + beta[,Nvar + (1:Nvar)] %*% (Data[t-1,]-Data[t-2,])) + rmvnorm(1, residMean, Sigma)
# Data[t,] <- ifelse(Data[t,] < lbound, lbound, Data[t,] )
# Data[t,] <- ifelse(Data[t,] > ubound, ubound, Data[t,] )
# }
#
# for (t in 1:totTime){
# Data[t,] <- Data[t,]*sd + mean
# }
#
# return(Data[-seq_len(warmup), ,drop=FALSE])
# }
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.