Nothing
set.seed(20160832)
# bivariate normal curve at different levels of probability: check value of curve at points;
# different margins
# higher dimensional input
# 2-d bvn no model fitting; check against theoretical values
rho <- 0.5
Mean <- c(3,12)
Sigma <- matrix(c(1,rho,rho,1),ncol=2)
n <- 10000
require(mvtnorm)
Sample <- rmvnorm(n,sigma=Sigma,mean=Mean)
colnames(Sample) <- c("x","y")
# check theoretical probabilities at calculated points on curve match input probabilities
p <- seq(0.5,0.95,by=0.05)
c1 <- lapply(p, function(p)JointExceedanceCurve(Sample,p))
c2 <- lapply(c1,function(o){
X <- cbind(o[[1]],o[[2]])
sapply(1:length(o[[1]]), function(i)pmvnorm(lower=as.numeric(X[i,]),
upper=rep(Inf, 2), mean=Mean, corr=Sigma))
})
res <- sapply(c2,mean)
expect_equal(res, target = p, tolerance=0.01,
info="JointExceedanceCurve: theoretical probabilities at calculated points on curve match input probabilities")
# check points on diagonal match theoretical quantile
tq <- sapply(p,function(p)qmvnorm(p, sigma = Sigma,mean=Mean, tail="upper")$quantile)
sq <- sapply(c1,function(o){
D <- (o[[1]]-o[[2]])^2
Min <- which.min(D)
o[[1]][Min]
})
expect_equal(sq, target = tq, tolerance=0.01,
info="JointExceedanceCurve: calculated points on diagonal match theoretical quantile")
# check that curve calculated properly at user specified points
p <- seq(0.05,0.01,by=-0.01)
x <- quantile(Sample[,1],seq(0.7,0.95,by=0.05))
names(x) <- NULL
y <- lapply(p, function(p)JointExceedanceCurve(Sample,p,x=x)[[2]])
z <- sapply(1:length(p), function(i)JointExceedanceCurve(Sample[,2:1],p[i],x=y[[i]])[[2]])
for(i in 1:length(p)){
expect_equal(z[,i],x,tol=0.01,info="JointExceedanceCurve: user specified points of x for curve calculation")
}
# check for estimation from fitted models from texmex
m <- mex(winter,mqu=0.7,dqu=0.7,which="NO")
m2 <- predict(m,nsim=n,pqu=0.9)
m3 <- predict(m,nsim=n,pqu=0.92)
j2a <- JointExceedanceCurve(m2,0.0005, which=c(1,3)) # columns of predict matrix, which have been re-ordered from original
j2b <- JointExceedanceCurve(m2,0.0005, which=c("O3","NO"))
j2c <- JointExceedanceCurve(m2,0.0005, which=c(3,4)) # columns of predict matrix, which have been re-ordered from original
j2d <- JointExceedanceCurve(m2,0.0005, which=c("NO","SO2"))
expect_equal(j2a,j2b,info="JointExceedanceCurve: specifying which pair to return by column number or name")
expect_equal(j2c,j2d,info="JointExceedanceCurve: specifying which pair to return by column number or name")
O3vals <- c(15,18,21,24)
j2e <- JointExceedanceCurve(m2,0.0005, which=c("O3","NO"),x=O3vals)
j2f <- JointExceedanceCurve(m2,0.0005, which=c("NO","O3"),x=j2e$NO)
expect_equal(j2f$O3,O3vals,tol=0.01,info="JointExceedanceCurve: user specified values of points at which to calc curve")
NO2vals <- 55:65 # these need to be chosen with care so that the two curves are actually estimable from the two different importance samples!
p <- 0.02
j2g <- JointExceedanceCurve(m2,p, which=c("NO2","NO"),x=NO2vals)
j2h <- JointExceedanceCurve(m3,p, which=c("NO2","NO"),x=NO2vals)
expect_equal(attributes(j2g)$ExceedanceProb,p,info="JointExceedanceCurve: attribute ExceedanceProb")
expect_equal(attributes(j2h)$ExceedanceProb,p,info="JointExceedanceCurve: attribute ExceedanceProb")
expect_equal(j2g,j2h,tol=0.05,info="JointExceedanceCurve: curves estimated from different importance samples get same answers")
# mexMCMC
m <- mexAll(winter,mqu=0.7,dqu=rep(0.7,5))
m4 <- mexMonteCarlo(nSample=5000,mexList=m)
NOvals <- seq(370,450,by=10)
p <- 0.01
j3a <- JointExceedanceCurve(m2,p, which=c("NO","NO2"),x=NOvals)
j3b <- JointExceedanceCurve(m4,p, which=c("NO","NO2"),x=NOvals)
expect_equal(attributes(j3a)$ExceedanceProb,p,info="JointExceedanceCurve: for mexMCMC object attribute ExceedanceProb")
expect_equal(attributes(j3a)$ExceedanceProb,p,info="JointExceedanceCurve: for mexMCMC object attribute ExceedanceProb")
expect_equal(j3a,j3b,tol=0.05,
info="JointExceedanceCurve: curves estimated from mexMCMC object and predict.mex object give same answers (up to fitted model accuracy and sampling variation")
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.