nn = 2000
mydat = data.frame(x = seq(0, 1, length=nn))
mydat = transform(mydat, mu = loge(-0+1.5*x+0.2*x^2, inverse=TRUE))
mydat = transform(mydat, mu = loge(0-sin(8*x), inverse=TRUE))
mydat = transform(mydat, y = rexp(nn, rate=1/mu))
(fit = vgam(y ~ s(x,df=5), amlexponential(w=c(0.001,0.1,0.5,5,60)),
mydat, trace=TRUE))
fit@extra
# These plots are against the sqrt scale (to increase clarity)
par(mfrow=c(1,2))
# Quantile plot
with(mydat, plot(x, sqrt(y), col="blue", las=1, main=
paste(paste(round(fit@extra$percentile, dig=1), collapse=", "),
"percentile-expectile curves")))
with(mydat, matlines(x, sqrt(fitted(fit)), lwd=2, col="blue", lty=1))
# Compare the fitted expectiles with the quantiles
with(mydat, plot(x, sqrt(y), col="blue", las=1, main=
paste(paste(round(fit@extra$percentile, dig=1), collapse=", "),
"percentile curves are red")))
with(mydat, matlines(x, sqrt(fitted(fit)), lwd=2, col="blue", lty=1))
for(ii in fit@extra$percentile)
with(mydat, matlines(x, sqrt(qexp(p=ii/100, rate=1/mu)), col="red"))
Run the code above in your browser using DataLab