Nothing
set.seed(20101111)
# large sample estimation, testing Normality of estimators with and without covariates in model
myTest <- function(Mu, Phi, Xi, ..., info){
ests <- lapply(1:nreps, function(i) evm(data=data.frame(y=x[,i],x=covariate), y=y,family=glo,...))
par <- sapply(1:nreps,function(i)ests[[i]]$par)
se <- sapply(1:nreps,function(i)ests[[i]]$se)
true <- cbind(Mu,Phi,Xi)
p.vals <- apply((t(par)-true)/t(se),2,function(x)shapiro.test(x)$p.value)
expect_true(min(p.vals) >= alpha, info=info)
}
#*************************************************************
# Set up simulation test info
ntests <- 21 #total number of parameters estimated in tests below:
nreps <- 20
nsim <- 10000
alpha <- 0.05/ntests # bonferroni
p <- matrix(runif(6*nreps, -1, 1),ncol=6)
p[,3] <- p[,3]
p[,5] <- p[,5]/3
p[,6] <- p[,6]/4
covariate <- seq(0.1,0.9,length=nsim)
#*************************************************************
# Test stationary model
x <- sapply(1:nreps,function(i)rglo(nsim,mu=p[i,1],sigma=exp(p[i,3]), xi=p[i,5]))
myTest(Mu=p[,1], Phi=p[,3], Xi=p[,5], info="glo: no covariates, p-value")
#*************************************************************
# Test covariate in mu only
x <- sapply(1:nreps,function(i)rglo(nsim,mu=p[i,1] + p[i,2]*covariate,sigma=exp(p[i,3]), xi=p[i,5]))
myTest(Mu=p[,1:2],Phi=p[,3], Xi=p[,5], mu=~x, info="glo: covariates in mu only, p-value")
#*************************************************************
# Test covariate in sigma only
x <- sapply(1:nreps,function(i)rglo(nsim,mu=p[i,1],sigma=exp(p[i,3]+p[i,4]*covariate), xi=p[i,5]))
myTest(Mu=p[,1],Phi=p[,3:4], Xi=p[,5], phi=~x,info="glo: covariates in phi only, p-value")
#*************************************************************
# Test covariate in xi only
x <- sapply(1:nreps,function(i)rglo(nsim,mu=p[i,1],sigma=exp(p[i,3]), xi=p[i,5]+p[i,6]*covariate))
myTest(Mu=p[,1],Phi=p[,3], Xi=p[,5:6], xi=~x, info="glo: covariates in xi only, p-value")
#*************************************************************
# Test covariate in all variables
x <- sapply(1:nreps,function(i)rglo(nsim,mu=p[i,1]+ p[i,2]*covariate,sigma=exp(p[i,3]+p[i,4]*covariate), xi=p[i,5]+p[i,6]*covariate))
myTest(Mu=p[,1:2],Phi=p[,3:4], Xi=p[,5:6], mu=~x,phi=~x,xi=~x, info="glo: covariates in all parameters, p-value")
##############################################################################
## pglo
set.seed(20101111)
probabilities <- runif(10)
xi.values <- c(0, seq(-0.3, 0.3, length.out=10))
core.sanity.test <- function(xi) {
randoms <- sort(rglo(10, 0, 1, xi))
expect_true(all(diff(pglo(randoms, 0, 1, xi)) >= 0),
info="pglo: ascending")
expect_true(all(diff(pglo(randoms, 0,1,xi,lower.tail=FALSE)) <= 0),
info="pglo: descending")
expect_equal(log(pglo(randoms, 0, 1,xi)),
pglo(randoms, 0, 1, xi, log.p=TRUE),
info="pglo: log.p (1)")
expect_equal(log(pglo(randoms, 0, 1, xi, lower.tail=FALSE)),
pglo(randoms, 0, 1, xi, lower.tail=FALSE, log.p=TRUE),
info="pglo: log.p (2)")
mu <- runif(1, -5, 5)
sigma <- rexp(1)
expect_equal(pglo(randoms, 0, 1,xi),
pglo(mu + sigma * randoms, mu, sigma, xi),
info="pglo: shift and scale")
} # Close core.sanity.tests
qglo.comparison <- function(xi) {
quantiles <- qglo(probabilities, 0, 1, xi)
my.probs <- pglo(quantiles, 0, 1, xi)
expect_equal(my.probs, probabilities, info="pglo:straighttest")
quantiles <- qglo(probabilities, 0, 1, xi, lower.tail=FALSE)
my.probs <- pglo(quantiles, 0, 1, xi, lower.tail=FALSE)
expect_equal(my.probs, probabilities, info="pglo:lowertail")
my.probs.2 <- pglo(quantiles, 0, 1, xi, lower.tail=TRUE)
expect_equal(probabilities+my.probs.2,
rep(1, length(probabilities)),
info="pglo: tail flip")
} # Close qglo.comparison
lapply(xi.values, core.sanity.test)
##lapply(xi.values, qglo.comparison)
##############################################################################
## qglo
## get the probabilities that we'll use and sort them
## into ascending order for safekeeping
set.seed(20101111)
probabilities <- sort(runif(10))
core.sanity.test <- function(xi) {
base.quantiles <- qglo(probabilities, 0, 1, xi)
## check that the values are ascending
expect_true(all(diff(base.quantiles)>=0), info="qglo:ascendingquantiles")
## and check that we're descending correctly for upper tail
bq2 <- qglo(probabilities, 0, 1, xi, lower.tail=FALSE)
expect_true(all(diff(bq2)<=0), info="qglo:descendingquantiles")
## does lower.tail work
expect_equal(base.quantiles,
qglo(1 - probabilities, 0, 1, xi, lower.tail=FALSE),
info="qglo: lower.tail works correctly")
## does log.p work?
expect_equal(base.quantiles,
qglo(log(probabilities), 0, 1, xi, log.p=TRUE),
info="qglo: log.p works")
## check shift and scale property
sigma <- rexp(1)
mu <- runif(1, -5, 5)
shifted <- mu + sigma * base.quantiles
expect_equal(shifted, qglo(probabilities, mu, sigma, xi),
info="qglo: shift and scale")
} # Close core.sanity.test
##lapply(c(0, seq(-5, 5, length.out=10)), core.sanity.test)
## known values
expect_equal(0, qglo(0.5, 0, 1, 0),
info="qglo: median match at zero xi")
xi <- seq(-2, 2, length=10)
expect_equal(qglo(0.5, 0, 1, xi),
rep(0,10),
info="qglo: median match at nonzero xi")
##############################################################################
## dglo
local.dglo <- function(x,loc,scale,shape){ # function is not vectorised
if(shape == 0){
Ex <- exp(-(x-loc)/scale)
out <- Ex/(scale* (1+Ex)^2)
} else {
Ox <- 1+shape/scale * (x-loc)
out <- Ox^{-1/shape - 1} / (scale * (1+ Ox^(-1/shape))^2)
}
out
}
myTest <- function(mu, sig, xi, info){
myd <- sapply(1:nreps, function(i) dglo(x[,i], mu[i],sig[i], xi[i]))
ed <- sapply(1:nreps, function(i) local.dglo(x[,i], loc=mu[i], scale=sig[i], shape=xi[i]))
expect_equal(ed, myd, info=info)
}
set.seed(20101111)
#*************************************************************
# Test dglo. Note that local.dglo is NOT vectorized.
nreps <- 100
nsim <- 1000
p <- matrix(runif(2*nreps, -1, 1),ncol=2)
p[, 1] <- p[, 1] + 1
mu <- rep(runif(nreps,0,1))
x <- sapply(1:nreps,
function(i)rglo(nsim,mu=mu[i],sigma=p[i,1], xi=p[i,2]))
myTest(mu=mu,sig=p[,1], xi=p[,2], info="dglo: random xi")
#*************************************************************
# Test dglo when some or all of xi == 0
p[sample(1:nreps,nreps/2),2] <- 0
x <- sapply(1:nreps,function(i)rglo(nsim,mu=mu[i],sigma=p[i,1],xi=p[i,2]))
myTest(mu=mu,sig=p[,1], xi=p[,2], info="dglo: some zero xi")
p[,2] <- 0
x <- sapply(1:nreps,function(i)rglo(nsim,mu=mu[i],sigma=p[i,1],xi=p[i,2]))
myTest(mu=mu,sig=p[,1], xi=p[,2], info="dglo: all zero xi")
#*************************************************************
# Test vectorization of dglo.
mu <- rnorm(nsim)
sig <- runif(nsim, 0, 2)
xi <- runif(nsim)
x <- rglo(nsim, mu, sig, xi)
myd <- dglo(x, mu, sig, xi)
ed <- sapply(1:nsim, function(i) local.dglo(x[i], loc=mu[i], scale=sig[i], shape=xi[i]))
expect_equal(ed, myd, info="dglo:vectorisation")
#*************************************************************
# test log.d argument
ld <- dglo(x,mu,sig,xi,log.d=TRUE)
expect_equal(myd, exp(ld), info="dglo:logdensity")
##############################################################################
## rglo
## so, how do we test an RNG...
num.simple <- 1000
num.quantile <- 1e6
xi.values <- c(0, seq(-5, 5, length.out=10))
test.quantiles <- c(0.25, 0.5, 0.75)
core.sanity.test <- function(xi) {
seed <- as.integer(runif(1, -1, 1)*(2**30))
set.seed(seed)
samples <- rglo(num.simple, 0, 1, xi)
expect_equal(length(samples), num.simple,
info = "rglo: output of correct length")
if (xi > 0) {
expect_true(all(samples>=-1/xi), info = "rglo:lowerboundcheck")
} else if (xi < 0) {
expect_true(all(samples<=-1/xi), info = "rglo:upperboundcheck")
}
## scale and shift property
sigma <- rexp(1)
mu <- runif(1, -5, 5)
shifted <- mu + sigma * samples
set.seed(seed)
expect_equal(shifted, rglo(num.simple, mu, sigma, xi), info = "rglo: scale and shift")
} # Close core.sanity.test
quantile.test <- function(xi) {
## here are the sampled quantiles
quantiles <- quantile(pglo(rglo(num.quantile, 0, 1, xi),
0, 1, xi),
probs=test.quantiles,
names=FALSE,na.rm=TRUE)
## this is a bit crude, but hey...
expect_equal(test.quantiles, quantiles, tolerance=0.02, info="rglo: quantile test")
} # Close quantile.test
lapply(xi.values, core.sanity.test)
lapply(xi.values, quantile.test)
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.