Nothing
# generates random sample from multivariate normal distribution
# with means mu and covariance sigma positive definite
# input: the sample size N, the parameters mu and sigma
rmultvariate <- function (N = 1, mu, Sigma, tol = 1e-06)
{
p <- length(mu)
if (!all(dim(Sigma) == c(p, p)))
stop("incompatible arguments")
Sd <- eigen(Sigma, symmetric = TRUE)
eva <- Sd$values
if (!all(eva >= -tol * abs(eva[1L])))
stop("'Sigma' is not positive definite")
X <- matrix(rnorm(p * N), N)
X <- drop(mu) + Sd$vectors %*% diag(sqrt(pmax(eva, 0)), p) %*% t(X)
if (N == 1)
drop(X) else t(X)
}
# Function to calculate nodes and weights of Gauss-Legendre
# quadrature, where n is the number of points. It uses the
# Golub-Welsh algorithm
GaussLegendre <- function(n)
{
n <- as.integer(n)
if (n < 0)
stop("Must be a non-negative number of nodes!")
if (n == 0)
return(list(x = numeric(0), w = numeric(0)))
i <- 1:n
j <- 1:(n-1)
mu0 <- 2
b <- j / (4 * j^2 - 1)^0.5
A <- rep(0, n * n)
A[(n + 1) * (j - 1) + 2] <- b
A[(n + 1) * j] <- b
dim(A) <- c(n, n)
sd <- eigen(A, symmetric = TRUE)
w <- rev(as.vector(sd$vectors[1, ]))
w <- mu0 * w^2
x <- rev(sd$values)
return(list(nodes = x, weights = w))
}
# Function to compute cumulative distribution function of the noncentral
# unilateral Dunnett's test statistic with infinity degrees of freedom nu
# given the vector of quantile q, the vector of correlation r, the vector
# of the noncentrality parameter delta, the vector of degrees of freedom
# nu, the number of points n of the Gauss-Legendre quadrature and the
# list x of nodes and weights of the quadrature. It returns cumulative
# probabilities.
pNDUD <- function(q, r, delta, n = 32, x)
{
k <- length(r)
# transformando de -1 a 1 para -infty a +infty, t = x/(1-x^2), x em -1 a 1
y <- matrix(x$nodes/(1-(x$nodes)^2), n, 1)
pArgYi <- function(y, r, delta)
{
ti <- pnorm((r^0.5 * y + q - delta)/(1-r)^0.5 , log.p=TRUE)
ti <- exp(sum(ti)+dnorm(y, log=TRUE))
return(ti)
}
I <- apply(y, 1, pArgYi, r, delta)
I <- I * ((1+x$nodes^2)/(1-(x$nodes)^2)^2)
I <- sum(I * x$weights)
return(I)
}
# Function to compute cumulative distribution function of the noncentral
# bilateral Dunnett's test statistic with infinity degrees of freedom nu
# given the vector of quantile q, the vector of correlation r, the vector
# of the noncentrality parameter delta, the vector of degrees of freedom
# nu, the number of points n of the Gauss-Legendre quadrature and the
# list x of nodes and weights of the quadrature. It returns cumulative
# probabilities.
pNDBD <- function(q, r, delta, n = 32, x)
{
k <- length(r)
# transform from [-1; 1] to (-infty, +infty), t = x/(1-x^2), x in [-1, 1]
y <- matrix(x$nodes/(1-(x$nodes)^2), n, 1)
pArgYi <- function(y, r, delta)
{
ti <- pnorm((r^0.5 * y + q - delta)/(1-r)^0.5 , log.p=FALSE)
ti <- ti - pnorm((r^0.5 * y - q - delta)/(1-r)^0.5, log.p=FALSE)
if (any(ti <= 0))
{
ti[ti <= 0] <- -743.7469
ti[ti > 0] <- log(ti[ti > 0])
} else ti <- log(ti)
ti <- exp(sum(ti)+dnorm(y, log=TRUE))
return(ti)
}
I <- apply(y, 1, pArgYi, r, delta)
I <- I * ((1+x$nodes^2)/(1-x$nodes^2)^2)
I <- sum(I * x$weights)
return(I)
}
# Function to compute the probability density function of the noncentral
# unilateral Dunnett's test statistic with infinity degrees of freedom nu
# given the vector of quantile q, the vector of correlation r, the vector
# of the noncentrality parameter delta, the vector of degrees of freedom
# nu, the number of points n of the Gauss-Legendre quadrature and the
# list x of nodes and weights of the quadrature. It returns density
# values.
dNDUD <- function(q, r, delta, n = 32, x)
{
k <- length(r)
# transform from [-1; 1] to (-infty, +infty), t = x/(1-x^2), x in [-1, 1]
y <- matrix(x$nodes/(1-(x$nodes)^2), n, 1)
rgdgi <- cbind(r, delta, (1:k))
loopg <- function(rd, ym, q)
{
tg <- pnorm((rd[1]^0.5 * ym + q - rd[2]) / (1 - rd[1])^0.5, log.p=TRUE)
return(tg)
}
loopi <- function(rd, ym, q)
{
ti <- dnorm((rd[1]^0.5 * ym + q - rd[2]) / (1 - rd[1])^0.5, log=TRUE) - 0.5 * log(1-rd[1])
return(ti)
}
# loop of y's
loopy <- function(ym, q, rd)
{
Phy <- apply(rd, 1, loopg, ym, q) # loopg
phy <- apply(rd, 1, loopi, ym, q) # loopi
sPhy <- sum(Phy) # sum of Phi's differences
sTi <- sum(exp(phy + sPhy - Phy)) * dnorm(ym)
return(sTi)
}
I <- apply(y, 1, loopy, q, rgdgi)
I <- I * ((1+x$nodes^2)/(1-x$nodes^2)^2)
I <- sum(I * x$weights)
return(I)
}
# Function to compute the probability density function of the noncentral
# bilateral Dunnett's test statistic with infinity degrees of freedom nu
# given the vector of quantile q, the vector of correlation r,the vector
# of the noncentrality parameter delta, the vector of degrees of freedom
# nu, the number of points n of the Gauss-Legendre quadrature and the
# list x of nodes and weights of the quadrature. It returns density
# probabilities.
dNDBD <- function(q, r, delta, n = 32, x)
{
k <- length(r)
# transform from [-1; 1] to (-infty, +infty), t = x/(1-x^2), x in [-1, 1]
y <- matrix(x$nodes/(1-(x$nodes)^2), n, 1)
rgdgi <- cbind(r, delta, (1:k))
loopg <- function(rd, ym, q)
{
tg <- pnorm((rd[1]^0.5 * ym + q - rd[2]) / (1 - rd[1])^0.5) -
pnorm((rd[1]^0.5 * ym - q - rd[2]) / (1 - rd[1])^0.5)
if (tg <= 0) tg <- -743.7469 else tg <- log(tg)
return(tg)
}
loopi <- function(rd, ym, q)
{
ti <- dnorm((rd[1]^0.5 * ym + q - rd[2]) / (1 - rd[1])^0.5) +
dnorm((rd[1]^0.5 * ym - q - rd[2]) / (1 - rd[1])^0.5)
if (ti <= 0) ti <- -743.7469 else ti <- log(ti)
aux <- ti - 0.5 * log(1 - rd[1])
return(aux)
}
# loop of y's
loopy <- function(ym, q, rd)
{
Phy <- apply(rd, 1, loopg, ym, q) # loopg
phy <- apply(rd, 1, loopi, ym, q) # loopi
sPhy <- sum(Phy) # sum of Phi's differences
sTi <- sum(exp(phy + sPhy - Phy)) * dnorm(ym)
return(sTi)
}
I <- apply(y, 1, loopy, q, rgdgi)
I <- I * ((1+x$nodes^2)/(1-x$nodes^2)^2)
I <- sum(I * x$weights)
return(I)
}
# Function to compute quantiles from the noncentral unilateral
# Dunnett's test statistic distribution with infinity degrees of freedom
# nu given the vector of quantile q, the vector of correlation r, the
# vector of the noncentrality parameter delta, the vector of degrees of
# freedom nu, the number of points n of the Gauss-Legendre quadrature
# and the list x of nodes and weights of the quadrature. It returns
# quantiles.
qNDUD <- function(p, r, delta, n = 32, x)
{
k <- length(r)
alpha <- 1 - p
alpha_k <- 1 - (1-alpha)^(1/k)
q0 <- qnorm(1 - alpha_k) # initial value
maxIt <- 5000
tol <- 1e-13
conv <- FALSE
it <- 0
while ((conv == FALSE) & (it <= maxIt))
{
q1 <- q0 - (pNDUD(q0, r, delta, n, x) - p) / dNDUD(q0, r, delta, n, x)
if (abs(q1-q0) <= tol) conv <- TRUE
q0 <- q1
it <- it + 1
}
return(q1)
}
# Function to compute quantiles from the noncentral bilateral
# Dunnett's test statistic distribution with infinity degrees of freedom
# nu given the vector of quantile q, the vector of correlation r, the
# vector of the noncentrality parameter delta, the vector of degrees of
# freedom nu, the number of points n of the Gauss-Legendre quadrature
# and the list x of nodes and weights of the quadrature. It returns
# quantiles.
qNDBD <- function(p, r, delta, n = 32, x)
{
k <- length(r)
alpha <- 1 - p
alpha_k <- 1 - (1-alpha/2)^(1/k)
q0 <- qnorm(1 - alpha_k) # initial value
if (q0 < 0) q0 <- -q0
maxIt <- 5000
tol <- 1e-13
conv <- FALSE
it <- 0
while ((conv == FALSE) & (it <= maxIt))
{
q1 <- q0 - (pNDBD(q0, r, delta, n, x) - p) / dNDBD(q0, r, delta, n, x)
if (abs(q1-q0) <= tol) conv <- TRUE
q0 <- q1
it <- it + 1
}
return(q1)
}
# Function to compute the Cumulative distribution of the noncentral
# unilateral Dunnett's test statistic with finite degrees of freedom
# nu given the vector of quantile q, the vector of correlation r, the
# vector of the noncentrality parameter delta, the vector of degrees of
# freedom nu, the number of points n of the Gauss-Legendre quadrature
# and the list x of nodes and weights of the quadrature. It returns
# cumulative probabilities.
pNDUDF <- function(q, r, nu, delta, n = 32, xx)
{
k <- length(r)
x <- xx$nodes
w <- xx$weights
# computing integral in [0,1]
y <- 0.5 * x + 0.5 # from [-1, 1] to [0, 1]
fx <- nu/2 * log(nu) - lgamma(nu/2) - (nu/2-1) * log(2) + (nu-1) * log(y) - nu * y * y / 2
fx <- exp(fx)
y <- matrix(y, n, 1) * q
fy <- 0.5 * apply(y, 1, pNDUD, r, delta, n, xx) * fx
I <- sum(fy * w)
# computing integral in [1, infty), using the transformation
# y <- -ln(x), 0 < x < exp(-1)
b <- exp(-1)
a <- 0
x <- (b - a) / 2 * x + (b + a) / 2 # from [-1,1] to [0, exp(-1)]
y <- -log(x) # from [0, exp(-1)] to [1, +infty)
fx <- nu/2 * log(nu) - lgamma(nu/2) - (nu/2-1) * log(2) + (nu-1) * log(y) - nu * y * y / 2
fx <- exp(fx)
y <- matrix(y, n, 1) * q
fy <- (b - a) / 2 * apply(y, 1, pNDUD, r, delta, n, xx) / x * fx
I <- I + sum(fy * w)
return(I)
}
# Function to compute the Cumulative distribution of the noncentral
# bilateral Dunnett's test statistic with finite degrees of freedom
# nu given the vector of quantile q, the vector of correlation r, the
# vector of the noncentrality parameter delta, the vector of degrees of
# freedom nu, the number of points n of the Gauss-Legendre quadrature
# and the list x of nodes and weights of the quadrature. It returns
# cumulative probabilities.
pNDBDF <- function(q, r, nu, delta, n = 32, xx)
{
k <- length(r)
x <- xx$nodes
w <- xx$weights
# computing integral in [0,1]
y <- 0.5 * x + 0.5 # from [-1,1] to [0, 1]
fx <- nu/2 * log(nu) - lgamma(nu/2) - (nu/2-1) * log(2) + (nu-1) * log(y) - nu * y * y / 2
fx <- exp(fx)
y <- matrix(y, n, 1) * q
fy <- 0.5 * apply(y, 1, pNDBD, r, delta, n, xx) * fx
I <- sum(fy * w)
# computing integral in [1, infty), using the transformation
# y <- -ln(x), 0 < x < exp(-1)
b <- exp(-1)
a <- 0
x <- (b - a) / 2 * x + (b + a) / 2 # from [-1,1] to [0, exp(-1)]
y <- -log(x) # from [0, exp(-1)] to [1, +infty)
fx <- nu/2 * log(nu) - lgamma(nu/2) - (nu/2-1) * log(2) + (nu-1) * log(y) - nu * y * y / 2
fx <- exp(fx)
y <- matrix(y, n, 1) * q
fy <- (b - a) / 2 * apply(y, 1, pNDBD, r, delta, n, xx) / x * fx
I <- I + sum(fy * w)
return(I)
}
# Function to compute the probability density function of the noncentral
# unilateral Dunnett's test statistic with finite degrees of freedom
# nu given the vector of quantile q, the vector of correlation r, the
# vector of the noncentrality parameter delta, the vector of degrees of
# freedom nu, the number of points n of the Gauss-Legendre quadrature
# and the list x of nodes and weights of the quadrature. It returns
# density values.
dNDUDF <- function(q, r, nu, delta, n = 32, xx)
{
k <- length(r)
x <- xx$nodes
w <- xx$weights
# computing integral in [0,1]
y <- 0.5 * x + 0.5 # from [-1,1] to [0, 1]
fx <- nu/2 * log(nu) - lgamma(nu/2) - (nu/2-1) * log(2) +
(nu-1) * log(y) - nu * y * y / 2 + log(y) # + jacobian
fx <- exp(fx)
y <- matrix(y, n, 1) * q
fy <- 0.5 * apply(y, 1, dNDUD, r, delta, n, xx) * fx
I <- sum(fy * w)
# computing integral in [1, infty), using the transformation
# y <- -ln(x), 0 < x < exp(-1)
b <- exp(-1)
a <- 0
x <- (b - a) / 2 * x + (b + a) / 2 # from [-1,1] to [0, exp(-1)]
y <- -log(x) # from [0, exp(-1)] to [1, +infty)
fx <- nu/2 * log(nu) - lgamma(nu/2) - (nu/2-1) * log(2) +
(nu-1) * log(y) - nu * y * y / 2 + log(y) # + jacobian
fx <- exp(fx)
y <- matrix(y, n, 1) * q
fy <- (b - a) / 2 * apply(y, 1, dNDUD, r, delta, n, xx) / x * fx
I <- I + sum(fy * w)
return(I)
}
# Function to compute the probability density function of the noncentral
# bilateral Dunnett's test statistic with finite degrees of freedom
# nu given the vector of quantile q, the vector of correlation r, the
# vector of the noncentrality parameter delta, the vector of degrees of
# freedom nu, the number of points n of the Gauss-Legendre quadrature
# and the list x of nodes and weights of the quadrature. It returns
# density values.
dNDBDF <- function(q, r, nu, delta, n = 32, xx)
{
k <- length(r)
x <- xx$nodes
w <- xx$weights
# computing integral in [0,1]
y <- 0.5 * x + 0.5 # from [-1,1] to [0, 1]
fx <- nu/2 * log(nu) - lgamma(nu/2) - (nu/2-1) * log(2) +
(nu-1) * log(y) - nu * y * y / 2 + log(y) # + jacobian
fx <- exp(fx)
y <- matrix(y, n, 1) * q
fy <- 0.5 * apply(y, 1, dNDBD, r, delta, n, xx) * fx
I <- sum(fy * w)
# computing integral in [1, infty), using the transformation
# y <- -ln(x), 0 < x < exp(-1)
b <- exp(-1)
a <- 0
x <- (b - a) / 2 * x + (b + a) / 2 # from [-1,1] to [0, exp(-1)]
y <- -log(x) # from [0, exp(-1)] to [1, +infty)
fx <- nu/2 * log(nu) - lgamma(nu/2) - (nu/2-1) * log(2) +
(nu-1) * log(y) - nu * y * y / 2 + log(y) # + jacobian
fx <- exp(fx)
y <- matrix(y, n, 1) * q
fy <- (b - a) / 2 * apply(y, 1, dNDBD, r, delta, n, xx) / x * fx
I <- I + sum(fy * w)
return(I)
}
# Function to compute quantiles from the noncentral
# unilateral Dunnett's test statistic with finite degrees of freedom
# nu given the vector of quantile q, the vector of correlation r, the
# vector of the noncentrality parameter delta, the vector of degrees of
# freedom nu, the number of points n of the Gauss-Legendre quadrature
# and the list x of nodes and weights of the quadrature. It returns
# quantiles.
qNDUDF <- function(p, r, nu, delta, n = 32, x)
{
k <- length(r)
alpha <- 1 - p
alpha_k <- 1 - (1-alpha)^(1/k)
q0 <- qt(1 - alpha_k, nu) + mean(delta)# initial value
p0 <- pNDUDF(q0, r, nu, delta, n, x)
if (abs(p0-p) > 0.1) q0 <- qNDUD(p, r, delta, n, x)
maxIt <- 5000
tol <- 1e-11
conv <- FALSE
it <- 0
while ((conv == FALSE) & (it <= maxIt))
{
q1 <- q0 - (pNDUDF(q0, r, nu, delta, n, x) - p) / dNDUDF(q0, r, nu, delta, n, x)
if (abs(q1-q0) <= tol) conv <- TRUE
if (q1 < 0) q1 <- -q1
q0 <- q1
it <- it + 1
}
return(q1)
}
# Function to compute quantiles from the noncentral
# bilateral Dunnett's test statistic with finite degrees of freedom
# nu given the vector of quantile q, the vector of correlation r, the
# vector of the noncentrality parameter delta, the vector of degrees of
# freedom nu, the number of points n of the Gauss-Legendre quadrature
# and the list x of nodes and weights of the quadrature. It returns
# quantiles.
qNDBDF <- function(p, r, nu, delta, n = 32, x)
{
k <- length(r)
alpha <- 1 - p
alphak <- 1 - (1-alpha/2)^(1/k)
q0 <- qt(1 - alphak, nu) + mean(delta) # valor
if (q0 < 0) q0 <- -q0
p0 <- pNDBDF(q0, r, nu, delta, n, x)
if (abs(p0-p) > 0.1) q0 <- qNDBD(p, r, delta, n, x)
maxIt <- 5000
tol <- 1e-11
conv <- FALSE
it <- 0
while ((conv == FALSE) & (it <= maxIt))
{
p0 <- pNDBDF(q0, r, nu, delta, n, x)
q1 <- q0 - (p0 - p) / dNDBDF(q0, r, nu, delta, n, x)
if (abs(q1-q0) <= abs(q0) * tol) conv <- TRUE
if (q1 < 0) q1 <- -q1
q0 <- q1
it <- it + 1
}
return(q1)
}
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.