#Chapter 10 countywinds functions
sWeib<-function(x,a,b,lambda=1)
{
 #returns rate of events exceeding x, for the Weibull marked Poisson process
 #returns vector of survival probabilities for Weibull distribution when lambda=1
 #a,b are location and shape vectors
 #lambda is the event frequency
 x<-(abs(x)+x)/2 #x is now 0 or positive
 lambda*exp(-(x/b)^a) 
}
fWeibMax<-function(x,a,b,lambda,t)
{
#returns distribution function of maximum over some interval t with rate lambda
exp(-t*sWeib(x,a,b,lambda))
}
rlWeibPois<-function(n,a,b,lambda,t=1)
{
 #returns a matrix or array of return levels for the Weibull 
 #marked Poisson process.
 #n an integer vector, n>0
 #t is a vector of interval lengths
 #Given a return level, rl, 
 #Divide up the time series by intervals of length t.
 #Any interval with any event of size rl or larger is a success.
 #Then n-1 is the average number of failures before the first success.
 #Also n-1 is the average number of failures between two successes.
 #Here we test n convert it to an integer. 
 if(any(n<=1)) stop("n must be greater than 1")
 n<-ceiling(n) #n is integer number of intervals. 
 np<-1/log(n/(n-1)) #close to n-.5
 #we want to handle matrices as well,
 #we have to collect attributes, convert to vectors.
 dims<-dim(lambda)
 dnames<-dimnames(lambda)
 b<-as.vector(b)
 a<-as.vector(a)
 lambda<-as.vector(lambda)+rep(0,length(a))+rep(0,length(b)) #replication of lambda
 out<-b*(log ((t*lambda) %o% np) )^(1/a) #generate a matrix note use of outer operator
 if(!is.null(dims))
 { #add attributes dims and dimnames
  dim(out)<-c(dims,length(np))
  if(is.null(dnames))
    dimnames(out)<-list(x=1:dims[1],y=1:dims[2],intervals=n)
  else
    dimnames(out)<-c(dnames,list(intervals=n))
 }
 else
 {
  colnames(out)<-n
 } 
 out
}

require(MASS)
sampleParameters<-function(R,fitc,fitw,X=list(lambda=1,a=1,b=1),n=c(2,5,10,20,50,100,200,500,1000),linkinv=list(lambda=exp,a=exp,b=exp),t=1)
{
 #R number of replicates
 #X is a list with three data frames of covariates named lambda, a, b 
 #Each data frame of X must include intercept for all terms as first column.
 #fitc fit count, Poisson count model output from glm 
 #fitw fit winds, Weibull wind model output from gamlss
 #n vector of return periods
 #linkinv a list of link terms by parameter
 #returns a  length(R)+1 by length(n) matrix of return levels, suitable for plotting
 #first row contains MLE
 #We have to have new data frames or matrices but we allow vectors. Note that the intercept is implied.
 
 X<-lapply(X,function(x) if(is.vector(x)) matrix(x,nrow=1))
 lambdacov<-vcov(fitc) 
 lambdacoef<-coef(fitc) 
 vcovfitw<-vcov(fitw,type="all") 
 thetalambda<-rbind(lambdacoef,mvrnorm(R,mu=lambdacoef,Sigma=lambdacov))
 thetaba<-rbind(vcovfitw$coef,mvrnorm(R,mu=vcovfitw$coef,Sigma=vcovfitw$vcov)) 
 lambda<-linkinv$lambda(thetalambda%*%t(X$lambda)) #
 dimnames(lambda)<-list(R=0:R,X=1:ncol(lambda)) #for call to return levels
 pb<-length(coef(fitw,"mu")) #get number of covariates for b
 a<-linkinv$a(thetaba[,-(1:pb),drop=F]%*%t(X$a)) #a parameters
 b<-linkinv$b(thetaba[,1:pb,drop=F]%*%t(X$b)) #b parameters
 #Note b parameters before a parameters
 #Now we have all the data to pass to the call rlWeibPois
 return(rlWeibPois(n=n,a=a,b=b,lambda=lambda,t=t))
}