Standard deviation
sdev.r=function(x)
{sqrt(var(x,na.rm=T))}
Skew
skew.r=function(x)
{
m=mean(x,na.rm=T)
x=x-m
s=sqrt(var(x,na.rm=T))
skew=mean(x*x*x,na.rm=T)/(s*s*s)
return(skew)
}
Sequent peak
seqpeak.r=function(q,r=0.7*mean(q))
{
# r is used as input so that for simulations that may have
different mean
# the same release may be used. r may be a vector to impose
a release pattern
rmq=r-q # net change in each time period. Also this
takes care of vector wrap around
# if r and q are of different lengths
n=length(q)
k=rep(0,n)
k[1]=max(rmq[1],0) # first time k[0] is 0
for (j in 2:n)
{
k[j]=max(k[j-1]+rmq[j],0)
}
s=max(k)
return(s)
}
Sequent peak use example (code fragment)
# This requires flows in summable units
# Use cfs.days as the unit and ignore leap years
days=c(31,30,31,31,28,31,30,31,30,31,30,31)
# This is oct, nov, dec, jan … sep for a water year
# Again dropping first 6 months
flowcfd=flow[7:1074]*days # again taking advantage of automatic
wrap around
y=c((1:9)/10,0.95)
s=rep(0,length(y))
for (j in 1:length(y)){
s[j]=seqpeak.r(flowcfd,y[j]*mean(flowcfd))
}
par("mar"=c(5.1,4.1,7.1,4.1)) # increase margin space for
extra axes
plot(s,y*mean(flowcfd),type="l",xlab="Storage cfs.day",ylab="Yield
cfs.day")
labs=c((0:50)/10)
axis(4,at=labs*mean(flowcfd),labels=labs)
mtext("Yield%",side=4,line=2)
axis(3,at=labs*mean(flowcfd)*12,labels=labs)
mtext("% Meanannualrunoff",side=3,line=2)
title("Deterministic Storage Yield Plot")
Boxplots (code fragment)
bmat=matrix(0,m,12) # matrix with m rows and 12 columns,1
for each month
for (i in 1:m){
fmatsim=matrix(flowsim[i,],ncol=12,byrow=T)
for(j in 1:12){
bmat[i,j]=mean(fmatsim[,j])
}
}
mno=c(10:12,1:9)
boxplot(split(bmat,col(bmat)),names=mno,ylab="cfs") # The split
function separates
# columns into a list that is needed by boxplot. names gives
the x axis labels
lines(1:12,fmave)
title("Mean Monthly Flow Boxplot")
Histograms (code fragment)
mno=c(10,11,12,1:9)
par(mfrow=c(4,3))
for (i in 1: 12){
hist(fmat[,i],main=paste("Month ",mno[i]),xlab="cfs")
box()
}
Normal
probability plot
normplot.r=function(x,labvals=c(0.1,0.2,0.5,0.8,0.9,0.95,0.99,0.995),b=0)
{
xs=sort(x)
n=length(xs)
i=1:n
p=(i-b)/(n+1-2*b) # General
plotting position, See Chow et al p396
q=qnorm(p) # Normal quantile
plot(q,xs,xaxt="n")
qv=qnorm(labvals)
axis(1,at=qv,label=labvals)
}
Gumbel probability plot
gumbelplot.r=function(x,labvals=c(0.1,0.2,0.5,0.8,0.9,0.95,0.99,0.995),b=0)
{
xs=sort(x)
n=length(xs)
i=1:n
p=(i-b)/(n+1-2*b) #
General plotting position, See Chow et al p396
q=-log(-log(p)) # Gumbel quantile
plot(q,xs,xaxt="n")
qv=-log(-log(labvals))
axis(1,at=qv,label=labvals)
}
# Generating an ARMA time series
armagen=function(phi=c(),theta=c(),siga=1,n=10000,wn=100)
{
p=length(phi) # The AR order of the model is the length of the first parameter vector input
q=length(theta) # The MA order of the model is the length of the first parameter vector input
nt=n+wn # the total number of time steps simulated
# wn is the number of warm up time steps discarded to avoid initialization issues
z=rep(0,nt)
a=rnorm(nt,sd=siga)
v=a
# Evaluate the moving average process
i=0
for(th in theta)
{
i=i+1
v[(1+i):nt]=v[(1+i):nt]-th*a[1:(nt-i)]
}
# Evaluate the AR process only
for order p and greater to avoid difficulty initializing
z=v
if(p > 0){
for(i in (p+1):nt){
for(j in
1:p)z[i]=z[i]+phi[j]*z[i-j]
}
}
return(z[(wn+1):nt])
}