#Set up the working environment
rm(list = ls())
library(reshape2)
library(xtable)
library(ggplot2)
# RUN IF YOU HAVE NOT INSTALLED THE BICALC PACKAGE
# library(devtools)
# install_github("bceaton/bicalc")
library(bicalc)
load("GSD_data.RData") #load the dataset
load("GSD_simulations.RData") #load the first dataset in the list
#Analysis in the paper The following code chunks contain the commands used to create all the figures in the paper, using the data files provided.
#sample.data = sample(all1.data$D, size = k)
d.demo = sort(sample(demo.data$D, size = 200))
cfd.demo = seq_along(d.demo)/length(d.demo)
# calculate the binomial limits for specified values n and p
n = 200
p = 0.84
firstINT = QuantBD(n, p)
#the exact binomial confidence interval, based on measurements
l = firstINT$interval[1]
u = firstINT$interval[2]
dl = d.demo[l]
du = d.demo[u]
#the equal tail approximation applied to the complete dataset
le = firstINT$equaltail[1]
ue = firstINT$equaltail[2]
dle = as.numeric(approx(x=seq_along(d.demo), y = d.demo, xout = le )[2])
due = as.numeric(approx(x=seq_along(d.demo), y = d.demo, xout = ue )[2])
#the equal tail approximation applied to the binned data
d.binned = MakeCFD(d.demo)
dlb = 2^as.numeric(approx(x=d.binned$probs, y = log2(d.binned$size), xout = le/200 )[2])
dub = 2^as.numeric(approx(x=d.binned$probs, y = log2(d.binned$size), xout = ue/200 )[2])
dest = 2^as.numeric(approx(x=d.binned$probs, y = log2(d.binned$size), xout = 0.84 )[2])
plot(d.demo, cfd.demo,
type ="b",
pch = 21,
#cex = 0.75,
log = "x",
col = rgb(0,0,1,0.35),
ylim = c(0, 1),
xlim = c(0.5, 8),
xlab = "grain size (mm)",
ylab = "Proportion Finer"
)
lines(MakeCFD(d.demo),
col = rgb(1,0,0,0.75),
lwd = 1.5,
lty =1)
lines(c(0.4, dub), c(ue/200, ue/200), lty = 2)
lines( c(dub, dub), c(-0.1, ue/200), lty = 2)
lines(c(0.4, dlb), c(le/200, le/200), lty = 2)
lines( c(dlb, dlb), c(-0.1, le/200), lty = 2)
lines(c(0.4, dest), c(0.84, 0.84), lty = 1)
lines( c(dest, dest), c(-0.1, 0.84), lty = 1)
k = 150:190
pr = dbinom(k, n, p)
df = data.frame(k, pr)
# find upper and lower limits and output coverage probability.
qci = QuantBD(n, p)
lu_exact = qci$interval
p_c = qci$coverage
df$bar_col = ifelse(df$k >= lu_exact[1] & df$k < lu_exact[2], "in", "out")
# using ggplot since barplots in base graphics are a bit awkward
par(mar=c(4,4,1,1))
par(oma = c(0,0,0,0))
ggplot(df, aes(x = k, y = pr)) +
geom_col(aes(fill = bar_col), na.rm = T) +
scale_x_continuous(limits = c(150, 190)) +
scale_y_continuous(limits = c(0, 0.08)) +
scale_fill_manual(values = c("darkgrey", "lightgrey"), name = "") +
labs(x = expression(italic(k)),
y = expression(italic(P[r](k, 200, 0.84)))) +
theme_bw() +
theme(axis.text = element_text(size = 12)) +
theme(axis.title = element_text(size = 12)) +
geom_vline(xintercept = firstINT$equaltail+0.5, lty = 2)
par(mfrow = c(1,2))
mycolors = c(rgb(0,0,1,0.4), rgb(0,1,1,0.4), rgb(1,1,0,0.4), rgb(1,0,0,0.4))
linecolors = c(rgb(0,0,1,0.6), rgb(0,1,1,1), rgb(1,1,0,1), rgb(1,0,0,0.6))
plot(cfd1hr$size, cfd1hr$probs,
type ="l",
lty = 0,
log = "x",
col = linecolors[1],
ylim = c(0, 1),
xlim = c(0.5, 5),
xlab = "grain size (mm)",
ylab = "Proportion Finer"
)
polygon(polyhr1,
col=mycolors[1],
lty = 0)
polygon(polyhr2,
col=mycolors[2],
lty = 0)
polygon(polyhr3,
col=mycolors[3],
lty = 0)
polygon(polyhr4,
col=mycolors[4],
lty = 0)
lines(cfd1hr$size, cfd1hr$probs,
col = linecolors[1],
lty = 2)
lines(cfd2hr$size, cfd2hr$probs,
col = linecolors[2],
lty = 2)
lines(cfd3hr$size, cfd3hr$probs,
col = linecolors[3],
lty = 2)
lines(cfd4hr$size, cfd4hr$probs,
col = linecolors[4],
lty = 2)
#abline(h = c(0.05, 0.95), lty = 2)
legend("topleft",
inset = 0.02,
cex = 0.75,
bg = rgb(1,1,1),
#box.lty = 0,
legend = c("sample 1", "sample 2", "sample 3", "sample 4"),
lty = c(0,0,0,0,2),
pch = c(15,15, 15, 15, NA),
pt.cex = c(2,2,2,2,1),
col = c(mycolors,
rgb(0,0,0))
)
plot(cfd5hr$size, cfd5hr$probs,
type ="l",
lty = 0,
log = "x",
col = rgb(0,0,1),
ylim = c(0, 1),
xlim = c(0.5, 5),
xlab = "grain size (mm)",
ylab = "Proportion Finer"
)
polygon(polyhr5,
col=mycolors[1],
lty = 0)
polygon(polyhr6,
col=mycolors[2],
lty = 0)
polygon(polyhr7,
col=mycolors[3],
lty = 0)
polygon(polyhr8,
col=mycolors[4],
lty = 0)
lines(cfd5hr$size, cfd5hr$probs,
col =linecolors[1],
lty = 2)
lines(cfd6hr$size, cfd6hr$probs,
col =linecolors[2],
lty = 2)
lines(cfd7hr$size, cfd7hr$probs,
col = linecolors[3],
lty = 2)
lines(cfd8hr$size, cfd8hr$probs,
col = linecolors[4],
lty = 2)
#abline(h = c(0.05, 0.95), lty = 2)
legend("topleft",
inset = 0.02,
cex = 0.75,
bg = rgb(1,1,1),
#box.lty = 0,
legend = c("sample 5", "sample 6", "sample 7", "sample 8"),
lty = c(0,0,0,0,2),
pch = c(15,15, 15, 15, NA),
pt.cex = c(2,2,2,2,1),
col = mycolors
)
par(mfrow = c(1,3))
ind = 2
boxplot(results[,ind],
ylim = c(ci$lower[ind], ci$upper[ind])*c(0.8, 1.2),
xlab = "")
title(xlab="D16", line=1, cex.lab=1.2)
abline(h = ci$estimate[ind], lty = 2, col = rgb(0,0,0))
abline(h = c(ci$upper[ind], ci$lower[ind]), lty = 2, col = rgb(1,0,0))
abline(h = c(ci50$upper[ind], ci50$lower[ind]), lty = 2, col = rgb(0,0,1))
d16.in95ci =sum(results[,ind]>ci$lower[ind] & results[,ind]<ci$upper[ind])/length(results[,ind])
d16.in50ci =sum(results[,ind]>ci50$lower[ind] & results[,ind]<ci50$upper[ind])/length(results[,ind])
ind = 4
boxplot(results[,ind],
ylim = c(ci$lower[ind], ci$upper[ind])*c(0.8, 1.2),
xlab = "")
title(xlab="D50", line=1, cex.lab=1.2)
abline(h = ci$estimate[ind], lty = 2, col = rgb(0,0,0))
abline(h = c(ci$upper[ind], ci$lower[ind]), lty = 2, col = rgb(1,0,0))
abline(h = c(ci50$upper[ind], ci50$lower[ind]), lty = 2, col = rgb(0,0,1))
d50.in95ci =sum(results[,ind]>ci$lower[ind] & results[,ind]<ci$upper[ind])/length(results[,ind])
d50.in50ci =sum(results[,ind]>ci50$lower[ind] & results[,ind]<ci50$upper[ind])/length(results[,ind])
ind = 6
boxplot(results[,ind],
ylim = c(ci$lower[ind], ci$upper[ind])*c(0.8, 1.2),
xlab = "")
title(xlab="D84", line=1, cex.lab=1.2)
abline(h = ci$estimate[ind], lty = 2, col = rgb(0,0,0))
abline(h = c(ci$upper[ind], ci$lower[ind]), lty = 2, col = rgb(1,0,0))
abline(h = c(ci50$upper[ind], ci50$lower[ind]), lty = 2, col = rgb(0,0,1))
d84.in95ci =sum(results[,ind]>ci$lower[ind] & results[,ind]<ci$upper[ind])/length(results[,ind])
d84.in50ci =sum(results[,ind]>ci50$lower[ind] & results[,ind]<ci50$upper[ind])/length(results[,ind])
par(mfcol=c(1,2))
plot(k.surf$size, k.surf$probs,
type ="b",
log = "x",
pch = 19,
col = "blue",
ylim = c(0, 1),
xlim = c(2,200),
xlab = "grain size (mm)",
ylab = "Proportion Finer"
)
text(x = 180, y = 0.01, labels = "A")
lines(k.bulk$size, k.bulk$probs, lty = 2, col = rgb(1,0,0))
points(k.bulk$size, k.bulk$probs, pch = 1, col = rgb(1,0,0))
plot(k.surf$size, k.surf$probs,
type ="p",
pch = 19,
log = "x",
col = "blue",
ylim = c(0, 1),
xlim = c(2,200),
xlab = "grain size (mm)",
ylab = "Proportion Finer"
)
polygon(poly.surf,
col=rgb(0, 0, 1,0.3),
lty = 0)
#abline(h = c(0.05, 0.95), lty = 2)
abline(v = c(22.6), lty = 2, col = rgb(1,0,0))
#lines(k.bulk$size, k.bulk$probs, lty = 2, col = rgb(1,0,0))
points(k.bulk$size, k.bulk$probs, pch = 1, col = rgb(1,0,0))
text(x = 180, y = 0.01, labels = "B")
legend("topleft",
inset = 0.02,
cex = 0.75,
bg = rgb(1,1,1),
#box.lty = 0,
legend = c("pebble count", "bulk sample (> 4 mm)", "22.6 mm"),
lty = c(0,0,2),
pch = c(19,1, NA),
col = c(rgb(0, 0, 1), rgb(1,0,0), rgb(1,0,0))
)
par(mfcol=c(1,2))
plot(cfdpool$size, cfdpool$probs,
type ="p",
pch = 19,
log = "x",
col = "blue",
ylim = c(0, 1),
xlim = c(2,400),
xlab = "grain size (mm)",
ylab = "Proportion Finer"
)
polygon(polypool,
col=rgb(0, 0, 1,0.3),
lty = 0)
#abline(h = c(0.05, 0.95), lty = 2)
abline(v = c(22.6), lty = 2, col = rgb(1,0,0))
text(x = 360, y = 0.01, labels = "A")
#lines(cfdbar$size, cfdbar$probs, lty = 1, col = rgb(0.5,0,1))
points(cfdbar$size, cfdbar$probs, pch = 19, col = rgb(0.5,0,1))
polygon(polybar,
col=rgb(0.5, 0, 1,0.3),
lty = 0)
#lines(cfdrun$size, cfdrun$probs, lty = 1, col = rgb(1,0.25,0))
points(cfdrun$size, cfdrun$probs, pch = 19, col = rgb(1,0.25,0))
polygon(polyrun,
col=rgb(1, 0.25, 0,0.3),
lty = 0)
plot(cfdpool4$size, cfdpool4$probs,
type ="p",
pch = 19,
log = "x",
col = "blue",
ylim = c(0, 1),
xlim = c(2,400),
xlab = "grain size (mm)",
ylab = "Proportion Finer"
)
polygon(polypool4,
col=rgb(0, 0, 1,0.3),
lty = 0)
#abline(h = c(0.05, 0.95), lty = 2)
abline(v = c(22.6), lty = 2, col = rgb(1,0,0))
text(x = 360, y = 0.01, labels = "B")
#lines(cfdbar4$size, cfdbar4$probs, lty = 1, col = rgb(0.5,0,1))
points(cfdbar4$size, cfdbar4$probs, pch = 19, col = rgb(0.5,0,1))
polygon(polybar4,
col=rgb(0.5, 0, 1,0.3),
lty = 0)
#lines(cfdrun4$size, cfdrun4$probs, lty = 1, col = rgb(1,0.25,0))
points(cfdrun4$size, cfdrun4$probs, pch = 19, col = rgb(1,0.25,0))
polygon(polyrun4,
col=rgb(1, 0.25, 0,0.3),
lty = 0)
legend("topleft",
inset = 0.02,
cex = 0.75,
bg = rgb(1,1,1),
#box.lty = 0,
legend = c("pool", "exposed bar", "run/riffle", "22.6 mm"),
lty = c(0,0,0,2),
pch = c(19,19, 19, NA),
col = c(rgb(0, 0, 1), rgb(0.5, 0, 1), rgb(1, 0.25, 0), rgb(1,0,0))
)
par(mfcol = c(1,2))
plot(cfdopA$size, cfdopA$probs,
type ="p",
pch = 1,
log = "x",
col = "blue",
ylim = c(0, 1),
xlim = c(2,400),
xlab = "grain size (mm)",
ylab = "Proportion Finer"
)
lines(cfdopA$size, cfdopA$probs, lty = 1, col = rgb(0,0,1))
lines(cfdopB$size, cfdopB$probs, lty = 1, col = rgb(1,0,0))
points(cfdopB$size, cfdopB$probs, pch = 19, col = rgb(1,0,0))
text(x = 360, y = 0.01, labels = "A")
plot(cfdopA$size, cfdopA$probs,
type ="p",
pch = 1,
log = "x",
col = "blue",
ylim = c(0, 1),
xlim = c(2,400),
xlab = "grain size (mm)",
ylab = "Proportion Finer"
)
polygon(polyopA,
col=rgb(0, 0, 1,0.3),
lty = 0)
#abline(h = c(0.05, 0.95), lty = 2)
#abline(v = c(16), lty = 1, col = rgb(1,0,0))
#lines(cfdbar$size, cfdbar$probs, lty = 1, col = rgb(1,0,0))
points(cfdopB$size, cfdopB$probs, pch = 19, col = rgb(1,0,0))
polygon(polyopB,
col=rgb(1, 0, 0,0.3),
lty = 0)
text(x = 360, y = 0.01, labels = "B")
legend("topleft",
inset = 0.02,
cex = 0.75,
bg = rgb(1,1,1),
#box.lty = 0,
legend = c("Operator A", "Operator B"),
lty = c(0,0),
pch = c(1,19),
col = c(rgb(0, 0, 1),rgb(1,0,0))
)
par(mfcol = c(1,2))
plot(cfdheel$size, cfdheel$probs,
type ="p",
pch = 1,
log = "x",
col = "blue",
ylim = c(0, 1),
xlim = c(2,400),
xlab = "grain size (mm)",
ylab = "Proportion Finer"
)
lines(cfdheel$size, cfdheel$probs, lty = 1, col = rgb(0,0,1))
lines(cfdframe$size, cfdframe$probs, lty = 1, col = rgb(1,0,0))
points(cfdframe$size, cfdframe$probs, pch = 19, col = rgb(1,0,0))
text(x = 360, y = 0.01, labels = "A")
plot(cfdheel$size, cfdheel$probs,
type ="p",
pch = 1,
log = "x",
col = "blue",
ylim = c(0, 1),
xlim = c(2,400),
xlab = "grain size (mm)",
ylab = "Proportion Finer"
)
polygon(polyheel,
col=rgb(0, 0, 1,0.3),
lty = 0)
#abline(h = c(0.05, 0.95), lty = 2)
#abline(v = c(16), lty = 1, col = rgb(1,0,0))
#lines(cfdbar$size, cfdbar$probs, lty = 1, col = rgb(1,0,0))
points(cfdframe$size, cfdframe$probs, pch = 19, col = rgb(1,0,0))
polygon(polyframe,
col=rgb(1, 0, 0,0.3),
lty = 0)
text(x = 360, y = 0.01, labels = "B")
legend("topleft",
inset = 0.02,
cex = 0.75,
bg = rgb(1,1,1),
#box.lty = 0,
legend = c("Heel-toe", "Sampling frame"),
lty = c(0,0),
pch = c(1,19),
col = c(rgb(0, 0, 1),rgb(1,0,0))
)
CalcError = function(p, X, n= c(100,200,300,400,500,600,700, 800, 900, 1000)){
tmp = rbind(WolmanCI(X,n[1],probs = p),
WolmanCI(X,n[2], probs =p),
WolmanCI(X,n[3], probs =p),
WolmanCI(X,n[4], probs =p),
WolmanCI(X,n[5], probs =p),
WolmanCI(X,n[6], probs =p),
WolmanCI(X,n[7], probs =p),
WolmanCI(X,n[8], probs =p),
WolmanCI(X,n[9], probs =p),
WolmanCI(X,n[10], probs =p)
)
err = as.numeric(0.5*(tmp$upper - tmp$lower)/tmp$estimate)
return(err)
}
CalcSDlog = function(X){
tmp = WolmanCI(X, 100, p = c(0.5, 0.84))
SDlog = diff(log2(tmp$estimate))
return(SDlog)
}
srcs = c("C1", "C2", "C3", "C4", "C5", "C6", "C7", "C8", "cfdopA", "cfdopB", "cfdpool", "cfdpool4", "cfdrun", "cfdrun4", "cfdbar", "cfdbar4", "k.surf", "cfd1hr", "cfd2hr", "cfd3hr", "cfd4hr", "cfd5hr", "cfd6hr", "cfd6hr", "cfd7hr", "cfd8hr", "G309A", "G309B")
Err50 = matrix(data=NA, nrow = 10 , ncol = length(srcs)+1)
Err84 = matrix(data=NA, nrow = 10, ncol = length(srcs)+1)
SDlog = numeric(length = length(srcs))
Err50[,1] = seq(100,1000,100)
Err84[,1] = Err50[,1]
for(i in seq_along(srcs)) {
Err50[,i+1] = CalcError(0.5, get(srcs[i]))
Err84[,i+1] = CalcError(0.84, get(srcs[i]))
SDlog[i] = CalcSDlog(get(srcs[i]))
}
n100 = data.frame(med50 = NA)
n100$med50 = mean(( Err50[1,2:19]))
n100$min50 = min(( Err50[1,2:19]))
n100$max50 = max(( Err50[1,2:19]))
n100$med84 = mean(( Err84[1,2:19]))
n100$min84 = min(( Err84[1,2:19]))
n100$max84 = max(( Err84[1,2:19]))
n200 = data.frame(med50 = NA)
n200$med50 = mean(( Err50[2,2:19]))
n200$min50 = min(( Err50[2,2:19]))
n200$max50 = max(( Err50[2,2:19]))
n200$med84 = mean(( Err84[2,2:19]))
n200$min84 = min(( Err84[2,2:19]))
n200$max84 = max(( Err84[2,2:19]))
n500 = data.frame(med50 = NA)
n500$med50 = mean(( Err50[5,2:19]))
n500$min50 = min(( Err50[5,2:19]))
n500$max50 = max(( Err50[5,2:19]))
n500$med84 = mean(( Err84[5,2:19]))
n500$min84 = min(( Err84[5,2:19]))
n500$max84 = max(( Err84[5,2:19]))
par(mfcol = c(1,2))
par(mar=c(4,4,1,1))
par(oma = c(0,0,0,0))
matplot(Err50[,1], Err50[,2:19],
type = "l",
col = rgb(0,0,1, 0.3),
lty = 1,
lwd = 2,
ylim = c(0, 0.6),
xlim = c(100,1000),
xlab = "sample size",
ylab = expression(epsilon[50])
)
abline(v = c(200,500), lty = 1, lwd = 4, col = rgb(0,0,0,0.3))
text(x = 950, y = 0.0, labels = "A")
matplot(Err84[,1], Err84[,2:19],
type = "l",
col = rgb(0,0,1, 0.3),
lty = 1,
lwd = 2,
ylim = c(0, 0.6),
xlim = c(100,1000),
xlab = "sample size",
ylab = expression(epsilon[84])
)
abline(v = c(200,500), lty = 1, lwd = 4, col = rgb(0,0,0,0.3))
text(x = 950, y = 0.0, labels = "B")
par(mfcol = c(1,2))
plot.cols = character(length = NSim)
mycolors = colorRampPalette(c("blue", "cyan", "yellow", "red"))(6)
filt = SErr50$logSD < 0.75
plot.cols[filt] = mycolors[1]
filt = SErr50$logSD >= 0.75 & SErr50$logSD < 1.0
plot.cols[filt] = mycolors[2]
filt = SErr50$logSD >= 1 & SErr50$logSD < 1.25
plot.cols[filt] = mycolors[3]
filt = SErr50$logSD >= 1.25 & SErr50$logSD < 1.5
plot.cols[filt] = mycolors[4]
filt = SErr50$logSD >= 1.5 & SErr50$logSD < 1.75
plot.cols[filt] = mycolors[5]
filt = SErr50$logSD >= 1.75 & SErr50$logSD < 2.0
plot.cols[filt] = mycolors[6]
plot(SErr50$N, SErr50$err,
type = "p",
pch = 20,
cex = 0.75,
col = plot.cols,
ylim = c(0, 0.6),
xlim=(c(0,1000)),
xlab = "sample size",
ylab = expression(epsilon[50])
)
X = seq(0,1000,1)
Y1 = exp(powfit50$coefficients[1]+powfit50$coefficients[3]*0.5)*X^powfit50$coefficients[2]
Y2 = exp(powfit50$coefficients[1]+powfit50$coefficients[3]*2.0)*X^powfit50$coefficients[2]
lines(X,Y2, lwd = 2, col = "black" )
lines(X,Y1, lwd =2, col = "black", lty = 2 )
text(x = 950, y = 0.0, labels = "A")
legend("topright",
inset = 0.02,
cex = 1.0,
bg = rgb(1,1,1),
#box.lty = 0,
legend = c(expression(2.0~phi),
expression(0.5~phi),
expression((0.5-0.75~phi)),
expression((0.75-1.0~phi)),
expression((1.0-1.25~phi)),
expression((1.25-1.5~phi)),
expression((1.5-1.75~phi)),
expression((1.75-2.0~phi))
),
lty = c(1,2,0,0,0,0,0,0),
lwd = c(2,2,0,0,0,0,0,0),
pch = c(NA,NA,20,20,20,20,20,20),
col = c("black", "black", mycolors)
)
plot(SErr84$N, SErr84$err,
type = "p",
pch = 20,
cex = 0.75,
col = plot.cols,
ylim = c(0, 0.6),
xlim=(c(0,1000)),
xlab = "sample size",
ylab = expression(epsilon[84])
)
X = seq(0,1000,1)
Y3 = exp(powfit84$coefficients[1]+powfit84$coefficients[3]*0.5)*X^powfit84$coefficients[2]
Y4 = exp(powfit84$coefficients[1]+powfit84$coefficients[3]*2.0)*X^powfit84$coefficients[2]
lines(X,Y4, lwd = 2, col = "black" )
lines(X,Y3, lwd = 2, col = "black", lty = 2 )
text(x = 950, y = 0.0, labels = "B")
Coef = c("A","B")
CoefTable = data.frame(Coef)
CoefTable$SD0.50 = c(coeff50[1] + coeff50[3]*0.5,
coeff84[1] + coeff84[3]*0.5)
CoefTable$SD0.75 = c(coeff50[1] + coeff50[3]*0.75,
coeff84[1] + coeff84[3]*0.75)
CoefTable$SD1.00 = c(coeff50[1] + coeff50[3]*1.0,
coeff84[1] + coeff84[3]*1.0)
CoefTable$SD1.25 = c(coeff50[1] + coeff50[3]*1.25,
coeff84[1] + coeff84[3]*1.25)
CoefTable$SD1.50 = c(coeff50[1] + coeff50[3]*1.5,
coeff84[1] + coeff84[3]*1.5)
CoefTable$SD1.75 = c(coeff50[1] + coeff50[3]*1.75,
coeff84[1] + coeff84[3]*1.75)
CoefTable$SD2.00 = c(coeff50[1] + coeff50[3]*2.0,
coeff84[1] + coeff84[3]*2.00)
#<<fishtrap, echo=FALSE>>=
D50 = 0.055
Tstar = 0.045
S = 0.02
Tcrit = 9.81*1650*Tstar*D50
dcrit = Tcrit/(9810*S)
epsilon100 = exp(coeff50[1] + coeff50[3]*1.0)*100^coeff50[2]
D50low = 0.055*(1-epsilon100)
D50up = 0.055*(1+epsilon100)
Tcritlow = 9.81*1650*Tstar*D50low
dcritlow = Tcritlow/(9810*S)
Tcritup = 9.81*1650*Tstar*D50up
dcritup = Tcritup/(9810*S)
Q = 2.5
h = 0.23*Q^0.4
offset = h - dcrit
hfull = offset + 2*dcrit
Qfull = (hfull/0.23)^(1/0.4)
hlow = offset + 2*dcritlow
Qlow = (hlow/0.23)^(1/0.4)
hup = offset + 2*dcritup
Qup = (hup/0.23)^(1/0.4)
# Specify Fishtrap Creek peak flows for 1972 - 2011
Q = c(14, 3.77, 8.72, 12, 4.11, 6.03, 4.33, 3.31, 3.78, 9.64,
9.55, 9.27, 4.13, 11.7, 3.6, 5.85, 3.98, 3.88, 10.5, 6.01, 9.67,
10.8, 15.3, 10.1, 13.5, 5.05, 3.13, 8.68, 4.2, 6.3, 10.1, 8.8, 8.02,
8.08, 8.02, 7.9, 9.73)
# Set up x axis tick positions and labels
Ttick = c(1.001,1.01,1.1,1.5,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,25,30,35,40,45,50,60,70,80,90,100)
# Fit a line by method of moments, along with 95% confidence intervals
KTtick = -(sqrt(6)/pi)*(0.5772 + log(log(Ttick/(Ttick-1))))
QTtick = mean(Q) + KTtick*sd(Q)
RPlow = as.numeric(approx(x = QTtick, y = Ttick, xout = Qlow)[2])
RPup = as.numeric(approx(x = QTtick, y = Ttick, xout = Qup)[2])
epsilon500 = exp(coeff50[1] + coeff50[3]*1.0)*500^coeff50[2]
D50low2 = 0.055*(1-epsilon500)
D50up2 = 0.055*(1+epsilon500)
Tcritlow2 = 9.81*1650*Tstar*D50low2
dcritlow2 = Tcritlow2/(9810*S)
Tcritup2 = 9.81*1650*Tstar*D50up2
dcritup2 = Tcritup2/(9810*S)
hlow2 = offset + 2*dcritlow2
Qlow2 = (hlow2/0.23)^(1/0.4)
hup2 = offset + 2*dcritup2
Qup2 = (hup2/0.23)^(1/0.4)
RPlow2 = as.numeric(approx(x = QTtick, y = Ttick, xout = Qlow2)[2])
RPup2 = as.numeric(approx(x = QTtick, y = Ttick, xout = Qup2)[2])