Flux (mass per unit time) and Load (mass) are computed from timeseries estimates of analyte concentration and discharge. Concentration timeseries are computed from regression between surrogate and user-provided analyte. Uncertainty in calculations is estimated using bootstrap resampling. Code for the processing of acoustic backscatter from horizontally profiling acoustic Doppler current profilers is provided. All methods detailed in Livsey et al (2020) <doi:10.1007/s12237-020-00734-z>, Livsey et al (2023) <doi:10.1029/2022WR033982>, and references therein.

Fits a multivariate value-added model (VAM), see Broatch, Green, and Karl (2018) <doi:10.32614/RJ-2018-033> and Broatch and Lohr (2012) <doi:10.3102/1076998610396900>, with normally distributed test scores and a binary outcome indicator. A pseudo-likelihood approach, Wolfinger (1993) <doi:10.1080/00949659308811554>, is used for the estimation of this joint generalized linear mixed model.  The inner loop of the pseudo-likelihood routine (estimation of a linear mixed model) occurs in the framework of the EM algorithm presented by Karl, Yang, and Lohr (2013) <DOI:10.1016/j.csda.2012.10.004>. This material is based upon work supported by the National Science Foundation under grants DRL-1336027 and DRL-1336265.

Calculate the probability density functions (PDFs) for two threshold evidence 
    accumulation models (EAMs). These are defined using the following Stochastic 
    Differential Equation (SDE), dx(t) = v(x(t),t)*dt+D(x(t),t)*dW, where x(t) is
    the accumulated evidence at time t, v(x(t),t) is the drift rate, D(x(t),t) is 
    the noise scale, and W is the standard Wiener process. The boundary conditions 
    of this process are the upper and lower decision thresholds, represented by b_u(t) 
    and b_l(t), respectively. Upper threshold b_u(t) > 0, while lower threshold b_l(t) < 0. 
    The initial condition of this process x(0) = z where b_l(t) < z < b_u(t). We 
    represent this as the relative start point w = z/(b_u(0)-b_l(0)), defined as 
    a ratio of the initial threshold location. This package generates the PDF using
    the same approach as the 'python' package it is based upon, 'PyBEAM' by Murrow and Holmes
    (2023) <doi:10.3758/s13428-023-02162-w>. First, it converts the SDE model into the 
    forwards Fokker-Planck equation dp(x,t)/dt = d(v(x,t)*p(x,t))/dt-0.5*d^2(D(x,t)^2*p(x,t))/dx^2,
    then solves this equation using the Crank-Nicolson method to determine p(x,t). Finally, 
    it calculates the flux at the decision thresholds, f_i(t) = 0.5*d(D(x,t)^2*p(x,t))/dx 
    evaluated at x = b_i(t), where i is the relevant decision threshold, either upper (i = u) 
    or lower (i = l). The flux at each thresholds f_i(t) is the PDF for each threshold, 
    specifically its PDF. We discuss further details of this approach in this package 
    and 'PyBEAM' publications. Additionally, one can calculate the cumulative distribution 
    functions of and sampling from the EAMs.