Dear R mailing list,
Greetings!!!
I'm doing a graphical test of suitability for the
aircondit data. I actually just followed the
exercise in Davison and Hinkley's book of Graphical
Test for model checking. The program for
testing if it comes from an exponential distribution
seems to be running okay but the program for
the gamma distribution distribution is running away
:-). I hope you can spare sometime to check
what I have written below. I am also in doubt about
the mle and functions I have written.
Thank you in advance for your help.
:-)
Jei
> ### Graphical test of suitability of the
exponential and gamma model for the aircondit data >
>
> library(boot)
>
> data(aircondit)
>
>
>
> # Testing if the aircondit data comes from an
exponential distribution>
> expqq.fun <- function(data, q) sort(data)/mean(data)
> exp.gen <- function(data, mle) rexp(length(data),
mle)> n <- nrow(aircondit)
> qq <- qexp((1:n)/(n+1))
> exp.boot <- boot(aircondit$hours, expqq.fun,R=999,
sim="parametric",
+ ran.gen=exp.gen, mle=1/mean(aircondit$hours),
q=qq)> env <- envelope(exp.boot, level=0.95)
> plot(qq, exp.boot$t0,xlab="Exponential quantiles",
+ ylab="Scaled order statistics",
xlim=c(0,max(qq)),
+
ylim=c(0,max(c(exp.boot$t0,env$overall[2,]))),pch=1)> lines(qq,env$overall[1,]); lines(qq,env$overall[2,])
> lines(qq,env$point[1,],lty=2);
lines(qq,env$point[2,],lty=2)>
>
>
> # Testing if the aircondit data comes from a gamma
distribution>
> gamqq.fun <- function(data, q) sort(data)/mean(data)
> gam.gen <- function(data, mle)
rgamma(length(data),shape=12, scale=108.0833/12)> n <- nrow(aircondit)
> p<-seq(.01,.99, by=.01)
> qq <- qgamma(p,shape=12,scale=108.0833/12)
> gam.boot <- boot(aircondit$hours, gamqq.fun,R=999,
sim="parametric",
+ ran.gen=gam.gen, mle=1/mean(aircondit$hours),
q=qq)> env <- envelope(gam.boot, level=0.95)
> plot(qq, gam.boot$t0,xlab="Gamma quantiles",
+ ylab="Scaled order statistics",
xlim=c(0,max(qq)),
+
ylim=c(0,max(c(gam.boot$t0,env$overall[2,]))),pch=1)
Error in xy.coords(x, y, xlabel, ylabel, log) :
x and y lengths differ> lines(qq,env$overall[1,]); lines(qq,env$overall[2,])
Error in xy.coords(x, y) : x and y lengths differ> lines(qq,env$point[1,],lty=2);
lines(qq,env$point[2,],lty=2)
Error in xy.coords(x, y) : x and y lengths differ>
