Installing bicalc

If you have not installed the bicalc package yet, there is a copy of it availble from GitHub. You will need to install the devtools package, if you have not already done so.

#  install devtools (uncomment the line below to install devtools)
# install.packages("devtools")

#  use devtools to install the bicalc package from BC Eaton's GitHub repository (uncomment to run)
# devtools::install_github("bceaton/bicalc")
 
 #  turn the bicalc library on
 library(bicalc)

Introduction to bicalc

This section of the markdown file introduces the main three functions in bicalc.

Generating some data

Let’s assume that you have recorded 200 individual b-axis measurements from stones on the surface of a gravel bed. We could enter those values in a variable called diameters (do so if you wish), or we could simply generate some randomly.

diameters = 2^rnorm(200, mean = 5, sd = 1)   #generate some log-normally distributed random numbers
head(diameters)   #print the first few data points to the console
## [1] 39.970116 26.258692 38.827837 24.039340 62.185430  7.483985

Now that we have some data to work with, we can use it to generate a cumulative frequency distribution of grain sizes. the function MakeCFD is used for this purpose.

  test.cfd = MakeCFD(diameters)
  head(test.cfd)   #print the CFD to the console
##        size probs
## 1  4.000000  0.00
## 2  5.656854  0.00
## 3  8.000000  0.03
## 4 11.313708  0.07
## 5 16.000000  0.16
## 6 22.627417  0.26

You can enter your own cumulative frequency data, but please ensure that it has the same format and variable names as the output generated by the code shown above.

Estimating the confidence bounds

Once we have the cumulative frequency distribution and the number of observations upon which it is based, we can calculate the confidence intervals about the estimate of any given percentile in the distribution. The function we use for this is WolmanCI, named in honour of the Reds Wolman, who first described the commonly used method of bed surface sampling used by many people today. Below, we use the function to estimate the \(D_{50}\) of the surface, and the 95% confidence bounds for that estimate.

D50.estimate = WolmanCI(cfd = test.cfd, n = 200, probs = c(0.5))
D50.estimate   #print the results for D50 to the console
##   probs estimate    lower    upper
## 1   0.5 33.84993 29.86333 38.36764

We can also do with for multiple percentiles of interest as follows.

D.estimates = WolmanCI(cfd = test.cfd, n = 200, probs = c(0.25, 0.5, 0.75))
D.estimates   #print the results for D25, D50 and D75 to the console
##   probs estimate    lower    upper
## 1  0.25 21.85664 17.66845 24.52301
## 2  0.50 33.84993 29.86333 38.36764
## 3  0.75 54.92579 48.23208 61.61601

Plotting the confidence intervals

let’s begin by first generating a conventional grain size distribution plot, then adding a polygon representing the confidence interval for the entire distribution. The polygon that includes the upper and lower estimates for each percentile in the distribution is created using PolyCI.

#  create the standard plot using base graphics
  plot(test.cfd,  
     type ="b",
     pch = 21,
     log = "x",
     col = rgb(0,0,1,0.35),
     ylim = c(0, 1),
     xlab = "grain size (mm)",
     ylab = "Proportion Finer"
     )

#  create the confidence interval polygon
test.poly = PolyCI(cfd = test.cfd, n = 200, alpha = 0.05)

#  add it to the basic plot
polygon(test.poly, col = rgb(0,0,1,0.2), lty = 0)