Nothing
# computes the CDF of the normal midrange
# input: the quantile q in R and the number of means size >= 2.
pNMRc <- function(q, size, np = 32)
{
if (q == -Inf) return(0.0) else if (q == Inf) return(1.0)
xx <- GaussLegendre_aux(np)
x <- xx$nodes
w <- xx$weights
# interval 1: (-infty; q - 8]
y <- (q - 8) + (1 + x) / (x - 1) # from -1, 1 to -infty, q - 8
aux <- pnorm(2*q-y)-pnorm(y)
aux[aux <= 0] <- .Machine$double.eps^19
fy <- log(dnorm(y)) + (size-1) * log(aux)
fy <- exp(fy)
fy <- size * (2 / (x - 1)^2) * fy
I <- sum(fy * w)
# interval 2: [q-8; q]
a <- q - 8
b <- q
y <- (b - a) / 2 * x + (a + b) / 2 # from -1, 1 to (q - 8), q
aux <- pnorm(2*q-y)-pnorm(y)
aux[aux <= 0] <- .Machine$double.eps^19
fy <- log(dnorm(y)) + (size-1) * log(aux)
fy <- exp(fy)
fy <- (b - a) / 2 * size * fy
I <- I + sum(fy * w)
return(I)
}
# auxiliar function 2
pMR_aux2c <- function(s, q, size, df, np = 32)
{
Ii <- pNMRc(s * q, size, np)
fx <- df/2 * log(df) - lgamma(df/2) - (df/2-1) * log(2) +
(df - 1) * log(s) - df * s * s / 2
fx <- exp(fx)
I <- fx * Ii
return(I)
}
# computes the CDF of the externally studentized normal midrange
# input: the quantile q in R and the number of means size >= 2,
# degrees of freedom df > 0, and number of points of the
# gaussian quadrature (np)
pMRc <- function(q, size, df, np = 32)
{
if (is.nan(q) == TRUE) return(NaN)
if (is.nan(size) == TRUE) return(NaN)
if (is.nan(df) == TRUE) return(NaN)
if (size<=1) return(NaN)
if (df <= 0) return(NaN)
if (q == Inf) return(1)
if (q == -Inf) return(0)
if (df == Inf) return(pNMRc(q, size, np))
xx <- GaussLegendre_aux(np)
x <- xx$nodes
w <- xx$weights
# xx <- GaussLegendre_aux(np)
# x <- xx$nodes
# w <- xx$weights
# Integration interval from 0 to 1
y <- -0.5 * x + 0.5 # from [-1; 1] to [0; 1]
y <- matrix(y, np, 1)
fy <- 0.5 * apply(y, 1, pMR_aux2c, q, size, df, np)
I <- sum(fy * w)
# integration interval from 1 to +infty
y <- 1+(1+x)/(1-x) # from [-1; 1] to [1; +infty)
y <- matrix(y, np, 1)
fy.aux <- apply(y, 1, pMR_aux2c, q, size, df, np)
fy <- log(fy.aux)+log(2)-2*log(1-x)
fy <- exp(fy)
I <- I+sum(fy * w)
return(I)
}
# computes the PDF of the normal midrange
# input: the quantile q in R and the number of means size >= 2,
# degrees of freedom df > 0 and the number of points of the
# gaussian quadrature (np). We use two divisions of the
# integration limits. The first division was between Inf and q-8,
# and the second division was between q-8 and q.
dNMRc <- function(q, size, np = 32)
{
xx <- GaussLegendre_aux(np)
x <- xx$nodes
w <- xx$weights
# interval 1: (-infty; q - 8]
y <- (q - 8) + (1 + x) / (x - 1) # from [-1; 1] to (-infty; q - 8]
aux <- pnorm(2 * q - y) - pnorm(y)
aux[aux <= 0] <- .Machine$double.eps^19
fy <- log(dnorm(y)) + (size - 2) * log(aux) + log(dnorm(2 * q - y)) + log(2) +
log(size) + log(size-1)
fy <- exp(fy)
fy <- (2 / (x - 1)^2) * fy
I <- sum(fy * w)
# interval 2: [q-8; q]
a <- q - 8
b <- q
y <- (b - a) / 2 * x + (a + b) / 2 #from [-1; 1] to [(q - 8); q]
aux <- pnorm(2 * q - y) - pnorm(y)
aux[aux <= 0] <- .Machine$double.eps^19
fy <- log(dnorm(y)) + (size - 2) * log(aux) + log(dnorm(2 * q - y)) + log(2) +
log(size) + log(size-1)
fy <- exp(fy)
fy <- (b - a) / 2 * fy
I <- I + sum(fy * w)
return(I)
}
# auxiliar function 3
pMR_aux3c <- function(s, q, size, df, np = 32)
{
Ii <- dNMRc(s * q, size, np)
fx <- df/2 * log(df) - lgamma(df/2) - (df/2-1) * log(2) +
(df - 1) * log(s) - df * s * s / 2
fx <- exp(fx) * s
I <- fx * Ii
return(I)
}
# computes the PDF of the externally studentized normal midrange
# input: the quantile q in R and the number of means size >= 2,
# degrees of freedom df > 0 and the number of points
# of the gaussian quadrature (np).
# We use two divisions of the integration limits.
# The first division was between Inf and q-8,
# and the second division was between q-8 and q.
dMRc <- function(q, size, df, np = 32)
{
if (is.nan(q) == TRUE) return(NaN)
if (is.nan(size) == TRUE) return(NaN)
if (is.nan(df) == TRUE) return(NaN)
if (size <= 1) return(NaN)
if (df <= 0) return(NaN)
if ((q == -Inf) | (q == Inf)) return(0)
if (df == Inf) return(dNMRc(q, size, np))
xx <- GaussLegendre_aux(np)
x <- xx$nodes
w <- xx$weights
# integration limit between 0 and 1
y <- -0.5 * x + 0.5 # changing the interval [-1,1] to [0,1]
y <- matrix(y, np, 1)
fy <- 0.5 * apply(y, 1, pMR_aux3c, q, size, df, np)
I <- sum(fy * w)
# integration limit [1,Infty] - transformation y <- -ln(x), 0 < x < exp(-1)
b <- exp(-1)
a <- 0.0
x <- (b - a) / 2.0 * x + (b + a) / 2.0 # changing the interval [-1,1] to [0,exp(0,-1)]
y <- -log(x) # changing the interval [0,exp(-1)] to [1,infty]
y <- matrix(y, np, 1)
fy <- (b - a) / 2.0 * apply(y, 1, pMR_aux3c, q, size, df, np) / x
I <- I + sum(fy * w)
return(I)
}
# computes the quantiles of the normal midrange
# inputs: the cumulative probability (0 < p < 1),
# the number of means (size >= 2) n and the number
# of points of the gaussian quadrature (np)
qNMRc <- function(p, size, np = 32, eps = 1.0e-13, maxit = 5000)
{
if (p < 0.5) q0 <- -0.5 else
if (p > 0.5) q0 <- 0.5 else q0 <- 0 # arbitrary initial value
found <- FALSE
it <- 0
while ((found == FALSE) & (it <= maxit))
{
q1 <- q0 - (pNMRc(q0, size, np) - p) / dNMRc(q0, size, np)
if (abs(q1-q0) <= eps) found = TRUE
it <- it + 1
q0 <- q1
}
return(q1)
}
# computes the quantiles of the externally studentized normal midrange
# inputs: the cumulative probability (0 < p < 1),
# the number of means (size >= 2) size, degrees of freedom df > 0
# and number of points of the gaussian quadrature (np)
qMRc <- function(p, size, df, np = 32, eps = 1e-13, maxit = 5000)
{
if (is.nan(p) == TRUE) return(NaN)
if (is.nan(size) == TRUE) return(NaN)
if (is.nan(df) == TRUE) return(NaN)
if (size <= 1) return(NaN)
if (df <= 0) return(NaN)
if (p == 1) return(Inf)
if (p == 0) return(-Inf)
if (p < 0) return(NaN) else if (p > 1) return(NaN)
if (df == Inf) return(qNMRc(p, size, np, eps, maxit))
if (p < 0.5) q0 <- -0.5 else
if (p > 0.5) q0 <- 0.5 else q0 <- 0 # arbitrary initial value
found <- FALSE
it <- 0
while ((found == FALSE) & (it <= maxit))
{
q1 <- q0 - (pMRc(q0, size, df, np) - p) / dMRc(q0, size, df, np)
if (abs(q1-q0) <= eps) found = TRUE
it <- it + 1
q0 <- q1
}
return(q1)
}
# computes the vector of random numbers of the normal midrange or
# externally studentized normal midrange.
# Inputs - n: size of the vector to be simulated, size: sample size
# df in the interval [0, Inf], is the degrees of freedom.
rMRc <- function(n, size, df = Inf)
{
if (is.nan(n) == TRUE) return(NaN)
if (is.nan(size) == TRUE) return(NaN)
if (is.nan(df) == TRUE) return(NaN)
if (n < 1) return(NaN)
if (size <= 1) return(NaN)
if (df <= 0) return(NaN)
if (df == Inf) X <- (matrix(rnorm(n * size), n, size)) else
X <- (matrix(rnorm(n * size), n, size)) / (rchisq(n, df) / df)^0.5
midrange <- function(x)
{
return((max(x)+ min(x)) / 2)
}
res <- apply(X, 1, midrange)
return(res)
}
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.