#Functions for working with county level wind data 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 collect attributes and 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)) #replicate 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 need 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=FALSE] %*% t(X$a)) #a parameters b = linkinv$b(thetaba[, 1:pb, drop=FALSE] %*% 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)) }