Nothing
# testing all of: plot.lp.evm* and plot.rl.evm* where * is opt, sim and boot
# first with no covariates
n <- 100
for(Family in list(gpd,gev)){
set.seed(20130513)
pst <- function(msg) texmex:::texmexPst(msg,Family=Family)
u <- switch(Family$name,GPD=14,GEV=-Inf)
data <- switch(Family$name,GPD=rain,GEV=portpirie$SeaLevel)
fit.opt <- evm(data,th=u,family=Family)
fit.sim <- evm(data,th=u,method="sim",trace=100000,family=Family)
fit.boot <- evmBoot(fit.opt,R=20,trace=30)
M <- seq(20,1000,len=20)
p.opt <- predict(fit.opt,M=M,ci=TRUE)
p.sim <- predict(fit.sim,M=M,ci=TRUE)
p.boot <- predict(fit.boot,M=M,ci=TRUE)
par(mfrow=c(3,3))
plot(p.opt,main=paste(Family$name,"\nMLE"))
plot(p.sim,main=paste(Family$name,"\nMCMC"))
plot(p.boot,main=paste(Family$name,"\nBootstrap"))
p.lp.opt <- predict(fit.opt,type="lp",ci=TRUE)
p.lp.sim <- predict(fit.sim,type="lp",ci=TRUE)
p.lp.boot <- predict(fit.boot,type="lp",ci=TRUE)
expect_error(plot(p.lp.opt), info=pst("plot.lp.evmOpt: fail if no covariates"))
expect_error(plot(p.lp.sim), info=pst("plot.lp.evmSim: fail if no covariates"))
expect_error(plot(p.lp.boot), info=pst("plot.lp.evmBoot: fail if no covariates"))
# now with covariates
n <- 1000
M <- 1000
mu <- 1
X <- data.frame(a = rnorm(n),b = runif(n,-0.1,0.1))
param <- switch(Family$name,GPD=cbind(X[,1],X[,2]),GEV=cbind(mu,X[,1],X[,2]))
th <- switch(Family$name,GPD=0,GEV=-Inf)
X$Y <- Family$rng(n,param,list(threshold=th))
start <- switch(Family$name,GPD=c(0,1,0,1),GEV=c(mu,0,1,0,1))
fit.opt <- evm(Y,data=X,phi=~a,xi=~b, th=th,family=Family,start=start)
fit.sim <- evm(Y,data=X,phi=~a,xi=~b, th=th,family=Family,method="sim",trace=100000,start=start)
o <- options(warn=-1)
fit.boot <- evmBoot(fit.opt,R=20,trace=30)
options(o)
nx <- 3
M <- seq(5,1000,len=20)
newX <- data.frame(a=rnorm(nx),b=runif(nx,-0.1,0.1))
p.opt <- predict(fit.opt,M=M,newdata=newX,ci=TRUE)
p.sim <- predict(fit.sim,M=M,newdata=newX,ci=TRUE)
p.boot <- predict(fit.boot,M=M,newdata=newX,ci=TRUE)
p.lp.opt <- predict(fit.opt,type="lp",newdata=newX,ci=TRUE)
p.lp.sim <- predict(fit.sim,type="lp",newdata=newX,ci=TRUE)
p.lp.boot <- predict(fit.boot,type="lp",newdata=newX,ci=TRUE)
par(mfrow=c(3,3))
plot(p.opt,sameAxes=FALSE,main=paste(Family$name,"MLE\ndifferent axes"))
plot(p.sim,sameAxes=FALSE,main=paste(Family$name,"MCMC\ndifferent axes"))
plot(p.boot,sameAxes=FALSE,main=paste(Family$name,"Bootstrap\ndifferent axes"))
plot(p.opt,sameAxes=TRUE,main=paste(Family$name,"MLE\nsame axes"))
plot(p.sim,sameAxes=TRUE,main=paste(Family$name,"MCMC\nsame axes"))
plot(p.boot,sameAxes=TRUE,main=paste(Family$name,"Bootstrap\nsame axes"))
par(mfrow=c(4,4))
plot(p.lp.opt,main=paste(Family$name,"MLE"))
plot(p.lp.sim,main=paste(Family$name,"MCMC\nmean"),type="mean")
plot(p.lp.sim,main=paste(Family$name,"MCMC\nmedian"),type="median")
plot(p.lp.boot,main=paste(Family$name,"Bootstrap"))
# single covariate only:
param <- switch(Family$name,GPD=cbind(X[1,1],X[,2]),GEV=cbind(mu,X[1,1],X[,2]))
X$Y <- Family$rng(n,param,list(threshold=th))
fit.opt <- evm(Y,data=X,xi=~b,th=th,family=Family)
fit.sim <- evm(Y,data=X,xi=~b,th=th,family=Family,method="sim",trace=100000)
o <- options(warn=-1)
fit.boot <- evmBoot(fit.opt,R=20,trace=30)
options(o)
p.opt <- predict(fit.opt,M=M,newdata=newX,ci=TRUE)
p.sim <- predict(fit.sim,M=M,newdata=newX,ci=TRUE)
p.boot <- predict(fit.boot,M=M,newdata=newX,ci=TRUE)
p.lp.opt <- predict(fit.opt,type="lp",newdata=newX,ci=TRUE)
p.lp.sim <- predict(fit.sim,type="lp",newdata=newX,ci=TRUE)
p.lp.boot <- predict(fit.boot,type="lp",newdata=newX,ci=TRUE)
par(mfrow=c(3,3))
plot(p.opt,sameAxes=FALSE,main=paste(Family$name,"MLE"))
plot(p.sim,sameAxes=FALSE,main=paste(Family$name,"MCMC"))
plot(p.boot,sameAxes=FALSE,main=paste(Family$name,"Bootstrap"))
par(mfrow=c(3,1))
plot(p.lp.opt,main=paste(Family$name,"MLE"),polycol="cyan")
plot(p.lp.sim,main=paste(Family$name,"MCMC"),polycol="cyan")
plot(p.lp.boot,main=paste(Family$name,"Bootstrap"),polycol="cyan")
}
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.