05-function-reference.qmd

# Function Reference

## Execute Reservoir Frequency Analysis

### `rfa_simulate()`

```
_R_F_A _S_t_a_g_e-_F_r_e_q_u_e_n_c_y _S_i_m_u_l_a_t_i_o_n

_D_e_s_c_r_i_p_t_i_o_n:

     Performs reservoir stage frequency analysis using Monte Carlo
     simulation to develop stage-frequency relationships. Combines
     stratified flow-frequency sampling, flood seasonality, hydrograph
     scaling, and Modified Puls reservoir routing to estimate annual
     exceedance probabilities of reservoir stage.

_U_s_a_g_e:

     rfa_simulate(
       sim_type = "expected",
       bestfit_params,
       dist = "LP3",
       stage_ts,
       seasonality,
       hydrographs,
       resmodel,
       Nbins = 50,
       events_per_bin = 200,
       routing_dt = 1,
       Ncores = NULL,
       sim_name = NULL,
       results_dir = NULL
     )
     
_A_r_g_u_m_e_n_t_s:

sim_type: Character string specifying the simulation mode. One of
          ‘"median"’, ‘"expected"’, or ‘"full"’. Default is
          ‘"expected"’.

bestfit_params: Data frame or matrix of distribution parameters from
          RMC-BestFit MCMC output. Columns 1-3 are distribution
          parameters:

            * LP3: mean (log), standard deviation (log), skew (log)

            * GEV: location, scale, shape

          Column 4 is log-likelihood (used only in ‘"median"’ mode to
          identify the posterior mode). For ‘"expected"’ mode, must
          have at least ‘Nbins * events_per_bin’ rows. For ‘"full"’
          mode, each row is treated as an independent realization.

    dist: Character string specifying the frequency distribution.
          Either ‘"LP3"’ (Log-Pearson Type III, default) or ‘"GEV"’
          (Generalized Extreme Value).

stage_ts: Data frame of historical reservoir stage with columns ‘date’
          (character in M/D/YYYY format) and ‘stage’ (numeric, feet).

seasonality: Numeric vector of length 12 giving monthly flood
          occurrence probabilities (relative frequencies). Used to
          sample the month of each simulated event, which determines
          the antecedent pool elevation.

hydrographs: List of hydrograph data frames as returned by
          ‘hydrograph_setup’. Each element has columns ‘datetime’,
          ‘hour’, ‘inflow’, and ‘hydrograph_num’, with attributes
          ‘"obs_vol"’ (observed max n-day volume) and ‘"dt"’ (timestep
          in hours). The list must have a ‘"probs"’ attribute
          containing normalized sampling probabilities.

resmodel: Data frame with three columns: elevation (ft), storage
          (acre-ft), and discharge (cfs). Defines the reservoir model
          for Modified Puls routing.

   Nbins: Integer. Number of stratified sampling bins in EV1 space.
          Default is 50.

events_per_bin: Integer. Number of events sampled per bin. Default is
          200. Total simulations per realization = ‘Nbins *
          events_per_bin’.

routing_dt: Numeric. Routing timestep in hours. Passed to
          ‘scale_hydrograph’ to resample hydrographs before routing.
          Supported values are ‘0.25’ (15-min), ‘1’ (1-hour, default),
          ‘6’ (6-hour), and ‘24’ (24-hour).

  Ncores: Integer or ‘NULL’. Number of parallel workers for ‘"full"’
          mode. If ‘NULL’ (default), automatically selects based on
          available cores (capped at 16). Ignored for ‘"median"’ and
          ‘"expected"’ modes.

sim_name: Optional character string to label the simulation. Used in
          the output CSV filename. If ‘NULL’ (default), defaults to
          ‘"sim"’. Spaces are replaced with underscores in the
          filename.

results_dir: Optional path to the directory where results are saved. If
          ‘NULL’ (default), creates an ‘rfaR_results’ folder in the
          current working directory.

_D_e_t_a_i_l_s:

     Three simulation modes are available:

     ‘"median"’ Single realization using the most likely (posterior
          mode) parameter set from RMC-BestFit. Produces one
          stage-frequency curve with no uncertainty bounds. Useful for
          quick screening or debugging.

     ‘"expected"’ Single realization using expected value sampling,
          where each stratified flow sample is paired 1-to-1 with a
          different parameter set from the posterior distribution.
          Produces one stage-frequency curve representing the expected
          (mean) result integrated across parameter uncertainty.

     ‘"full"’ Nested Monte Carlo with full uncertainty quantification.
          The outer loop iterates over all parameter sets (one per row
          of ‘bestfit_params’), and each inner loop independently
          samples natural variability via stratified sampling,
          seasonality, starting pool, and hydrograph shape. Each
          realization produces an independent stage-frequency curve;
          curves are combined on a common stage grid to yield expected,
          median, and 5th/95th percentile confidence bounds. The outer
          loop is parallelized using the ‘future’ framework.

     For all modes, natural variability is captured through stratified
     sampling in with ‘Nbins’ bins and ‘events_per_bin’ events per bin,
     providing reliable coverage from common events down to
     approximately 1e-8 AEP.

     The stratified sampling approach divides the AEP range (default
     0.99 to 1e-8) into equal-width bins in EV1 (Gumbel reduced
     variate) space and samples events uniformly within each bin. This
     ensures adequate representation of rare flood events that would
     require orders of magnitude more samples under crude Monte Carlo
     sampling.

     Post-processing uses the law of total probability via
     ‘stage_frequency_curve’: for a grid of stage thresholds, the
     weighted exceedance fraction is computed within each bin and
     summed across bins.

     In ‘"full"’ mode, each parallel worker independently samples all
     random inputs (seasonality, starting pool, hydrograph shape, and
     stratified z-variates) in addition to using its own parameter set.
     This ensures full Monte Carlo independence across realizations.
     Each worker post-processes its results into a stage-frequency
     curve immediately, returning only the curve rather than the full
     peak stage matrix, keeping memory usage manageable for large
     numbers of realizations.

     Curves from individual realizations are combined by interpolating
     each onto a common stage grid (spanning the global min/max across
     all realizations) in log10(AEP) space, then computing summary
     statistics across realizations at each stage threshold.

     Results are automatically exported as a CSV file named
     ‘{sim_name}_{sim_type}_{MM_DD_YY_HHMM}.csv’ in the ‘results_dir’
     directory. For example, a median simulation named "John McGraw
     Dam" run on June 4, 2025 at 2:30 PM would produce
     ‘John_McGraw_Dam_median_06_04_25_1430.csv’.

_V_a_l_u_e:

     A list whose contents depend on ‘sim_type’:

     For ‘"median"’ and ‘"expected"’:

     stage_frequency Data frame with columns ‘stage’ and ‘AEP’ from
          ‘stage_frequency_curve’.

     peakStage Matrix of simulated peak stages (‘events_per_bin’ rows
          by ‘Nbins’ columns).

     peakFlow Matrix of simulated peak discharges (same dimensions).

     weights Stratified sampling bin weights.

     For ‘"full"’:

     stage_frequency Data frame with columns ‘stage’, ‘expected’,
          ‘median’, ‘lower_05’, and ‘upper_95’ on a common stage grid.

     all_curves List of per-realization stage-frequency data frames
          (each with ‘stage’ and ‘AEP’ columns).

     aep_matrix Matrix of interpolated AEP values (500 rows by
          ‘Nrealizations’ columns) on the common stage grid.

     common_stage Numeric vector of the common stage grid (length 500).

     Nrealizations Number of parameter set realizations processed.

     Nsims_per_realiz Number of simulations per realization.

_S_e_e _A_l_s_o:

     ‘flow_frequency_sampler’, ‘flow_frequency_sampler_expected’,
     ‘stratified_sampler’, ‘scale_hydrograph’, ‘mod_puls_routing’,
     ‘stage_frequency_curve’, ‘hydrograph_setup’

_E_x_a_m_p_l_e_s:

     ## Not run:
     
     # Not run: Expected and Median modes take ~15s; Full uncert
     # runs ~10,000 parameter realizations and takes 3-4 hours on a
     # typical machine.
     
     # --- Setup ---
     hydros <- hydrograph_setup(jmd_hydro_apr1999, jmd_hydro_jun1965,
                                jmd_hydro_may1955, jmd_hydro_pmf,
                                critical_duration = 2, routing_days = 10)
     
     # --- Expected only (default) ---
     results_exp <- rfa_simulate(
       sim_type       = "expected",
       bestfit_params = jmd_bf_parameter_sets,
       stage_ts       = jmd_wy1980_stage,
       seasonality    = jmd_seasonality$relative_frequency,
       hydrographs    = hydros,
       resmodel       = jmd_resmodel,
       sim_name       = "jmd"
     )
     
     # --- Median only ---
     results_med <- rfa_simulate(
       sim_type       = "median",
       bestfit_params = jmd_bf_parameter_sets,
       stage_ts       = jmd_wy1980_stage,
       seasonality    = jmd_seasonality$relative_frequency,
       hydrographs    = hydros,
       resmodel       = jmd_resmodel,
       sim_name       = "jmd"
     )
     
     # --- Full uncertainty (parallelized) ---
     results_full <- rfa_simulate(
       sim_type       = "full",
       bestfit_params = jmd_bf_parameter_sets,
       stage_ts       = jmd_wy1980_stage,
       seasonality    = jmd_seasonality$relative_frequency,
       hydrographs    = hydros,
       resmodel       = jmd_resmodel,
       Ncores         = 4,
       sim_name       = "jmd"
     )
     ## End(Not run)
     
```

## Primary RFA Modules

### `stratified_sampler()`

```
_S_t_r_a_t_i_f_i_e_d _S_a_m_p_l_e_r

_D_e_s_c_r_i_p_t_i_o_n:

     Creates stratified bins in standard normal space for use in
     stratified sampling. Stratification improves sampling efficiency
     by ensuring adequate coverage of rare events.

_U_s_a_g_e:

     stratified_sampler(
       minAEP = 1e-08,
       maxAEP = 0.99,
       dist = "EV1",
       Nbins = NULL,
       Mevents = NULL,
       verbose = FALSE
     )
     
_A_r_g_u_m_e_n_t_s:

  minAEP: Minimum annual exceedance probability. Default is ‘1E-8’.

  maxAEP: Maximum annual exceedance probability. Default is ‘0.99’.

    dist: Character. Probability space for stratification.
          Case-insensitive. Default is ‘"EV1"’. Invalid values trigger
          a warning and default to ‘"EV1"’.

          ‘"EV1"’ Extreme Value Type I (Gumbel) space. Recommended for
              flood frequency analysis. Allocates more bins to rare
              events through the transformation ‘-log(-log(1-AEP))’,
              improving tail estimation efficiency.

          ‘"Normal"’ Standard normal (z-score) space. Uniform bins in
              z-space. Use when the underlying phenomenon is normally
              distributed.

          ‘"Uniform"’ Uniform probability space. Equal probability
              width per bin. Generally inefficient for rare event
              estimation.

   Nbins: Number of stratified bins. Default is ‘20’.

 Mevents: Number of events per bin. Default is ‘500’.

 verbose: Logical. If ‘TRUE’, prints a completion message summarizing
          the stratification. Default is ‘FALSE’.

_D_e_t_a_i_l_s:

     The function divides the probability space into bins using the
     specified transformation. EV1 transformation is recommended for
     heavy-tailed distributions common in flood frequency analysis, as
     it naturally allocates more sampling effort to rare events
     critical for dam safety assessments.

_V_a_l_u_e:

     A list containing:

     normOrd Vector of standard normal ordinates spanning all bins

     Weights Vector of probability weights for each bin

     Zlower Vector of lower bounds for each bin (standard normal)

     Zupper Vector of upper bounds for each bin (standard normal)

     Nbins Number of bins

     Mevents Number of events per bin

_E_x_a_m_p_l_e_s:

     # Default stratification
     strat <- stratified_sampler()
     
     # Custom bins and events
     strat <- stratified_sampler(minAEP = 1E-6,
                                 maxAEP = 0.5,
                                 dist = "EV1",
                                 Nbins = 10,
                                 Mevents = 100)
     
     # With verbose message
     strat <- stratified_sampler(verbose = TRUE)
     
```

### `flow_frequency_sampler()`

```
_F_l_o_w _F_r_e_q_u_e_n_c_y _S_a_m_p_l_e_r

_D_e_s_c_r_i_p_t_i_o_n:

     Generates a stratified matrix of flow values from a single set of
     frequency distribution parameters using stratified Monte Carlo
     sampling. Used internally by ‘rfa_simulate()’ for median-only and
     full uncertainty modes.

_U_s_a_g_e:

     flow_frequency_sampler(
       bestfit_params,
       freq_dist = "LP3",
       strat_dist = "ev1",
       Nbin = NULL,
       Mevent = NULL
     )
     
_A_r_g_u_m_e_n_t_s:

bestfit_params: Numeric vector of length 3 containing distribution
          parameters. For LP3: ‘c(mean_log, sd_log, skew_log)’. For
          GEV: ‘c(location, scale, shape)’.

freq_dist: Character. Distribution type. Either ‘"LP3"’ (default) or
          ‘"GEV"’.

strat_dist: Character. Probability space for stratification bins.
          Passed to ‘stratified_sampler()’. One of ‘"ev1"’ (default),
          ‘"normal"’, or ‘"uniform"’. See ‘stratified_sampler()’ for
          details.

    Nbin: Integer. Number of stratified bins. Default is ‘50’.

  Mevent: Integer. Number of events per bin. Default is ‘200’.

_V_a_l_u_e:

     A list containing:

     flow Matrix of sampled flow values [Mevent x Nbin]

     nbins Number of stratified bins

     mevents Number of events per bin

     weights Probability weights for each bin from
          ‘stratified_sampler()’

_S_e_e _A_l_s_o:

     ‘stratified_sampler()’, ‘qp3()’, ‘rfa_simulate()’

_E_x_a_m_p_l_e_s:

     # Single parameter set (posterior mode)
     params <- c(4.85, 0.39, -0.15)
     result <- flow_frequency_sampler(params, freq_dist = "LP3",
                                      Nbin = 20, Mevent = 500)
     # Dimensions of result
     dim(result$flow)  # 500 x 20
     
     # Distribution of Sampled Flows
     hist(result$flow)
     
```

### `flow_frequency_sampler_expected()`

```
_F_l_o_w _F_r_e_q_u_e_n_c_y _S_a_m_p_l_e_r (_E_x_p_e_c_t_e_d _O_n_l_y)

_D_e_s_c_r_i_p_t_i_o_n:

     Generates a stratified matrix of flow values by pairing each
     z-ordinate with a different parameter set from the posterior
     distribution. This collapses the nested Monte Carlo structure into
     a single pass, simultaneously sampling natural variability (via
     stratified z-ordinates) and knowledge uncertainty (via varying
     parameters). Used internally by ‘rfa_simulate()’ for expected-only
     mode.

_U_s_a_g_e:

     flow_frequency_sampler_expected(
       bestfit_params,
       freq_dist = "LP3",
       strat_dist = "ev1",
       Nbin = NULL,
       Mevent = NULL
     )
     
_A_r_g_u_m_e_n_t_s:

bestfit_params: Data frame of distribution parameters from RMC-BestFit.
          Must have ‘Nbin * Mevent’ rows (one parameter set per
          z-ordinate). For LP3: columns are mean (log), sd (log), skew
          (log). For GEV: columns are location, scale, shape.

freq_dist: Character. Distribution type. Either ‘"LP3"’ (default) or
          ‘"GEV"’.

strat_dist: Character. Probability space for stratification bins.
          Passed to ‘stratified_sampler()’. One of ‘"ev1"’ (default),
          ‘"normal"’, or ‘"uniform"’. See ‘stratified_sampler()’ for
          details.

    Nbin: Integer. Number of stratified bins. Default is ‘50’.

  Mevent: Integer. Number of events per bin. Default is ‘200’.

_V_a_l_u_e:

     A list containing:

     flow Matrix of sampled flow values [Mevent x Nbin]

     nbins Number of stratified bins

     mevents Number of events per bin

     weights Probability weights for each bin from
          ‘stratified_sampler()’

_S_e_e _A_l_s_o:

     ‘stratified_sampler()’, ‘qp3()’, ‘flow_frequency_sampler()’,
     ‘rfa_simulate()’

_E_x_a_m_p_l_e_s:

     # Using a pre-loaded parameter set (all 10,000 parameter sets)
     result <- flow_frequency_sampler_expected(jmd_bf_parameter_sets,
                                               Nbin = 20, Mevent = 500)
     dim(result$flow)  # 500 x 20
     
     # Using bootstrapped parameter samples
     jmd_samples <- bootstrap_vfc(
       c(jmd_vfc_parameters$mean_log,
         jmd_vfc_parameters$sd_log,
         jmd_vfc_parameters$skew_log),
       dist = "LP3",
       ERL  = jmd_vfc_parameters$erl)
     
     jmd_result <- flow_frequency_sampler_expected(
       jmd_samples$params,
       freq_dist = "LP3",
       Nbin      = 20,
       Mevent    = 500)
     
```

### `hydrograph_setup()`

```
_H_y_d_r_o_g_r_a_p_h _S_e_t_u_p _f_o_r _R_F_A _S_i_m_u_l_a_t_i_o_n

_D_e_s_c_r_i_p_t_i_o_n:

     Prepares hydrograph data frames copied from RMC-RFA for use in RFA
     simulation. Converts date/time columns and adds sequential hour
     and hydrograph ID columns.

_U_s_a_g_e:

     hydrograph_setup(
       ...,
       critical_duration = NULL,
       routing_days = NULL,
       weights = NULL
     )
     
_A_r_g_u_m_e_n_t_s:

     ...: Data frame representing an input hydrograph with columns:
          Ord, Date, Time, Flow (cfs). Copied directly from RMC-RFA.
          Date & Time should be ‘class() = "character"’

critical_duration: Critical duration in days.

routing_days: Desired length of routing simulation in days.

 weights: Optional numeric vector of sampling weights for each
          hydrograph. Must be the same length as the number of input
          hydrographs. Weights are normalized to probabilities
          internally. If NULL (default), all hydrographs are weighted
          equally.

_V_a_l_u_e:

     A list of formatted hydrograph data frames, each containing:

     datetime Date-time (POSIXct)

     hour Hours from start of event

     inflow Inflow (cfs)

     hydrograph_num Hydrograph ID number for sampling

     obs_vol Max n-day inflow volume (stored as an attribute of each
          dataframe)

     dt Hydrograph timestep (delta time, dt) in hours (stored as an
          attribute of each dataframe)

     The returned list also has a ‘probs’ attribute containing the
     normalized sampling probabilities derived from ‘weights’.

_E_x_a_m_p_l_e_s:

     # Setup with equal weights (default)
     hydros <- hydrograph_setup(jmd_hydro_apr1999,
                                jmd_hydro_jun1965,
                                jmd_hydro_pmf,
                                critical_duration = 2,
                                routing_days = 10)
     
     # Setup with custom weights (PMF 3x more likely to be sampled)
     hydros <- hydrograph_setup(jmd_hydro_apr1999,
                                jmd_hydro_jun1965,
                                jmd_hydro_pmf,
                                critical_duration = 2,
                                routing_days = 10,
                                weights = c(1, 1, 3))
     
     # view normalized probabilities
     attr(hydros, "probs")
     
```

### `scale_hydrograph()`

```
_S_c_a_l_e _H_y_d_r_o_g_r_a_p_h

_D_e_s_c_r_i_p_t_i_o_n:

     Scales inflow hydrograph given a sampled inflow volume from the
     volume-frequency curve. Scale factor is defined by the sampled
     inflow volume and the maximum volume of the corresponding
     duration. If the native hydrograph timestep differs from the
     target routing timestep, resampling is applied prior to scaling:
     finer-to-coarser uses block averaging; coarser-to-finer uses
     linear interpolation.

_U_s_a_g_e:

     scale_hydrograph(
       hydrograph_shape,
       observed_volume,
       sampled_volume,
       routing_dt = 1
     )
     
_A_r_g_u_m_e_n_t_s:

hydrograph_shape: Data frame with two columns: time (hours) and inflow
          (cfs). Must be in that order.

observed_volume: Maximum n-day inflow volume from input hydrograph
          shape.

sampled_volume: Sampled inflow volume from volume-frequency curve.

routing_dt: Target routing timestep in hours. Supported values are
          ‘0.25’ (15-min), ‘1’ (1-hour, default), ‘6’ (6-hour), and
          ‘24’ (24-hour).

_V_a_l_u_e:

     A data frame with two columns: ‘time_hrs’ and ‘inflow_cfs’, at the
     target routing timestep, with time starting at hour 0.

_E_x_a_m_p_l_e_s:

     # Example hydrograph. Requires pre-processing
     hydro_example <- hydrograph_setup(jmd_hydro_jun1965_15min, critical_duration = 2, routing_days = 10)
     hydrograph_shape <- hydro_example[[1]][, 2:3]
     
     # Default 1-hour routing timestep
     scaled <- scale_hydrograph(hydrograph_shape,
                                observed_volume = 50000,
                                sampled_volume = 55000)
     
     # 15-min input hydrograph, 1-hour routing timestep
     scaled <- scale_hydrograph(hydrograph_shape,
                                observed_volume = 50000,
                                sampled_volume = 55000,
                                routing_dt = 1)
     
```

### `mod_puls_routing()`

```
_M_o_d_i_f_i_e_d _P_u_l_s _R_e_s_e_r_v_o_i_r _R_o_u_t_i_n_g

_D_e_s_c_r_i_p_t_i_o_n:

     Performs Modified Puls (level pool) routing of inflow hydrograph
     given a defined reservoir geometry (Stage (ft), Storage (ac-ft),
     Discharge (cfs)).

_U_s_a_g_e:

     mod_puls_routing(resmodel_df, inflow_df, initial_elev, full_results = FALSE)
     
_A_r_g_u_m_e_n_t_s:

resmodel_df: Data frame with three columns: elevation/stage (ft),
          storage (acre-feet), and discharge (cfs). Must be in that
          order.

inflow_df: Data frame with two columns: time (hours) and inflow (cfs).
          Must be in that order.

initial_elev: Starting water surface elevation in feet.

full_results: Logical. If ‘FALSE’ (default), returns only peak stage
          and discharge. If ‘TRUE’, returns the complete routing
          result.

_V_a_l_u_e:

     If ‘full_results = FALSE’, a named numeric vector with
     ‘peak_stage_ft’ and ‘peak_discharge_cfs’. If ‘full_results =
     TRUE’, a data frame with columns: ‘time_hr’, ‘inflow_cfs’,
     ‘elevation_ft’, ‘storage_acft’, and ‘outflow_cfs’.

_R_e_f_e_r_e_n_c_e_s:

     Chow, V.T. (1959). Open-Channel Hydraulics. McGraw-Hill.

_E_x_a_m_p_l_e_s:

     # Example hydrograph. Requires pre-processing
     hydro_example <- hydrograph_setup(jmd_hydro_jun1965_15min,
                                        critical_duration = 2,
                                        routing_days = 10)
     hydrograph_shape <- hydro_example[[1]][, 2:3]
     scaled_hydrograph <- scale_hydrograph(hydrograph_shape,
                                           observed_volume = 50000,
                                           sampled_volume = 55000)
     
     # Peak values only
     mod_puls_routing(jmd_resmodel, scaled_hydrograph, initial_elev = 3830)
     
     # Full routing table
     jmd_full_routing <- mod_puls_routing(jmd_resmodel,
                                          scaled_hydrograph,
                                          initial_elev = 3830,
                                          full_results = TRUE)
     head(jmd_full_routing)
     
```

### `stage_frequency_curve()`

```
_C_o_m_p_u_t_e _S_t_a_g_e-_F_r_e_q_u_e_n_c_y _C_u_r_v_e

_D_e_s_c_r_i_p_t_i_o_n:

     Converts a matrix of routed peak stages into a stage-frequency
     curve using the law of total probability across stratified bins.
     For each candidate stage, the conditional exceedance probability
     within each bin is weighted by the bin probability and summed to
     produce the unconditional annual exceedance probability (AEP).

_U_s_a_g_e:

     stage_frequency_curve(peakStage, weights, stage_bins = 1000)
     
_A_r_g_u_m_e_n_t_s:

peakStage: Matrix of peak stages [Mevents x Nbins] from routing.

 weights: Numeric vector of bin weights from ‘stratified_sampler()’.
          Must have length equal to ‘ncol(peakStage)’.

stage_bins: Numeric value of exceedance stages used to compute
          stage-frequency in each bin. Default is 1000.

_V_a_l_u_e:

     A data frame with columns:

     stage Evaluated stage values

     AEP Annual exceedance probability at each stage

_S_e_e _A_l_s_o:

     ‘stratified_sampler()’, ‘flow_frequency_sampler()’,
     ‘rfa_simulate()’

```

## Rejection Sampling

### `pmf_stage_lognormal()`

```
_P_a_r_a_m_e_t_e_r_i_z_e _a _T_h_r_e_e-_P_a_r_a_m_e_t_e_r _L_o_g_n_o_r_m_a_l _P_M_F _f_o_r _S_t_a_g_e

_D_e_s_c_r_i_p_t_i_o_n:

     Computes the parameters of a three-parameter lognormal
     distribution for use as a probabilistic maximum stage (PMF) in
     rejection sampling. The shift parameter defines the lower or upper
     bound of the distribution (typically, the lower value). The
     function supports two modes: (1) supplying a best estimate and
     assuming sigma, or (2) supplying a best estimate, low, and high to
     solve for sigma numerically.

_U_s_a_g_e:

     pmf_stage_lognormal(
       pmf_shift,
       pmf_mean,
       pmf_sigma = NULL,
       pmf_low = NULL,
       pmf_high = NULL
     )
     
_A_r_g_u_m_e_n_t_s:

pmf_shift: Numeric. Assumed lower-bound of the lognormal distribution.
          This will define the shift value. Must be less than all of
          ‘pmf_mean’, ‘pmf_low’, and ‘pmf_high’.

pmf_mean: Numeric. Assumed mean of the shifted PMF stage distribution.

pmf_sigma: Numeric. Assumed standard deviation on the log scale.
          Required if ‘pmf_low’ and ‘pmf_high’ are not supplied.
          Defaults to ‘NULL’. Suggested value = 0.5.

 pmf_low: Numeric. Low estimate of PMF stage (assumed 5th percentile).
          Optional. If supplied, ‘pmf_high’ must also be supplied.

pmf_high: Numeric. High estimate of PMF stage (assumed 95th
          percentile). Optional. If supplied, ‘pmf_low’ must also be
          supplied.

_V_a_l_u_e:

     A named list with the following elements:

     pmf_shift Hard lower bound of the distribution.

     pmf_mean Assumed mean of shifted PMF stage distribution.

     pmf_sigma Standard deviation on the log scale.

     pmf_mu Derived location parameter on the log scale.

     pmf_p05 5th percentile of the PMF stage distribution.

     pmf_p95 95th percentile of the PMF stage distribution.

_E_x_a_m_p_l_e_s:

     # Mode 1: assume sigma
     pmf <- pmf_stage_lognormal(pmf_shift = 239, pmf_mean = 241.9, pmf_sigma = 0.5)
     
     # Mode 2: solve for sigma from low/high estimates
     pmf <- pmf_stage_lognormal(pmf_shift = 239, pmf_mean = 241.9,
                                 pmf_low = 239, pmf_high = 245)
     
```

### `rejection_sampling_stage()`

```
_R_e_j_e_c_t_i_o_n _S_a_m_p_l_i_n_g _o_f _S_t_a_g_e _B_o_u_n_d_e_d _b_y _a _P_r_o_b_a_b_i_l_i_s_t_i_c _M_a_x_i_m_u_m _S_t_a_g_e

_D_e_s_c_r_i_p_t_i_o_n:

     Draws ‘n_samples’ stage values from a stage-frequency curve,
     rejecting any draws that exceed a probabilistic maximum stage
     (PMF) sampled from a three-parameter lognormal distribution.
     Accepted samples are ranked and assigned Weibull plotting
     positions for use in frequency analysis.

_U_s_a_g_e:

     rejection_sampling_stage(
       pmf_stage_LN,
       stage_freq_df = NULL,
       aep = NULL,
       stage = NULL,
       n_samples = 1e+07
     )
     
_A_r_g_u_m_e_n_t_s:

pmf_stage_LN: A named list returned by ‘pmf_stage_lognormal’ containing
          the lognormal PMF parameters.

stage_freq_df: A data frame with AEP in column 1 and stage in column 2.
          Mutually exclusive with ‘aep’ and ‘stage’.

     aep: Numeric vector of annual exceedance probabilities. Must be
          supplied with ‘stage’. Mutually exclusive with
          ‘stage_freq_df’.

   stage: Numeric vector of stages corresponding to ‘aep’. Must be
          supplied with ‘aep’. Mutually exclusive with ‘stage_freq_df’.

n_samples: Integer. Number of samples to draw. Defaults to ‘1e7’.

_V_a_l_u_e:

     A data frame of thinned accepted samples with columns for z-score
     and stage, suitable for plotting a probabilistically bounded
     stage-frequency curve. Output is produced by ‘thin_samples()’.

_S_e_e _A_l_s_o:

     ‘pmf_stage_lognormal’

_E_x_a_m_p_l_e_s:

     ## Not run:
     
     pmf_ln <- pmf_stage_lognormal(pmf_shift = 239, pmf_mean = 241.9, pmf_sigma = 0.5)
     
     result <- rejection_sampling_stage(
       pmf_stage_LN  = pmf_ln,
       stage_freq_df = my_stage_freq_df,
       n_samples     = 1e7
     )
     ## End(Not run)
     
```

## Additional Utilities

### `bootstrap_vfc()`

```
_B_o_o_t_s_t_r_a_p _o_p_t_i_o_n _f_o_r _V_F_C (_R_F_A-_s_t_y_l_e _w/ _E_R_L)

_D_e_s_c_r_i_p_t_i_o_n:

     Generates matrix of parameters using posterior mode/mean and ERL

_U_s_a_g_e:

     bootstrap_vfc(bestfit_postmode, dist = "LP3", ERL = 150, Nboots = 10000)
     
_A_r_g_u_m_e_n_t_s:

bestfit_postmode: Vector of distribution parameters from RMC-BestFit.
          For LP3: columns are mean (log), sd (log), skew (log). For
          GEV: columns are location, scale, shape.

    dist: Distribution type. Either ‘"LP3"’ (default) or ‘"GEV"’.

     ERL: psuedo effective record length for bootstrapping

  Nboots: number of bootstraps

_V_a_l_u_e:

     A list containing:

     params Matrix of bootstrapped parameters

     dist Distribution selected (LP3 or GEV)

     postmode Posterior mode parameters provided to bootstrap

     ERL pseudo effective record length provided to bootstrap

_S_e_e _A_l_s_o:

     ‘flow_frequency_sampler()’

_E_x_a_m_p_l_e_s:

     # Sample using JMD VFC parameters
     jmd_samples <- bootstrap_vfc(c(jmd_vfc_parameters$mean_log,
                                     jmd_vfc_parameters$sd_log,
                                     jmd_vfc_parameters$skew_log),
                                     dist = "LP3",
                                     ERL = jmd_vfc_parameters$erl)
     
     # GEV Example with parent distribution parameters
     gev_example <- c(3.0, 1.0, -0.1)
     gev_samples <- bootstrap_vfc(gev_example,
                                   dist = "GEV",
                                   ERL = 200,
                                   Nboots = 5000)
     hist(gev_samples$params[,1])
     
     # LP3 Example with parent distribution parameters
     lp3_example <- c(3.5, 0.22, 0.1)
     lp3_samples <- bootstrap_vfc(lp3_example,
                                   dist = "LP3",
                                   ERL = 300,
                                   Nboots = 1000)
     hist(lp3_samples$params[,3])
     
```

### `qp3()`

```
_P_e_a_r_s_o_n _T_y_p_e _I_I_I _I_n_v_e_r_s_e _C_D_F (_Q_u_a_n_t_i_l_e _F_u_n_c_t_i_o_n)

_D_e_s_c_r_i_p_t_i_o_n:

     Computes quantiles from a Pearson Type III distribution given
     probabilities and distribution parameters (mean, standard
     deviation, skewness).

_U_s_a_g_e:

     qp3(p, mu, sigma, gamma)
     
_A_r_g_u_m_e_n_t_s:

       p: Vector of probabilities (between 0 and 1).

      mu: Mean of the distribution.

   sigma: Standard deviation of the distribution.

   gamma: Skewness coefficient of the distribution.

_D_e_t_a_i_l_s:

     NOTE: This is the non-vectorized version of the function (it will
     be used in a loop) Vectorized version is part of future
     development.

     The Pearson Type III distribution is parameterized by mean (‘mu’),
     standard deviation (‘sigma’), and skewness (‘gamma’). These are
     converted internally to location, scale, and shape parameters.

     When skewness is near zero (‘abs(gamma) < 1E-3’), the normal
     distribution is used as an approximation.

_V_a_l_u_e:

     Vector of quantiles corresponding to the input probabilities.

_E_x_a_m_p_l_e_s:

     # Single quantile
     qp3(0.99, mu = 10, sigma = 2, gamma = 0.5)
     
     # Multiple quantiles
     qp3(c(0.5, 0.9, 0.99), mu = 10, sigma = 2, gamma = 0.5)
     
     # Use with JMD volume-frequency parameters
     qp3(0.99,
         mu = jmd_vfc_parameters$mean_log,
         sigma = jmd_vfc_parameters$sd_log,
         gamma = jmd_vfc_parameters$skew_log)
     
```

### `aep2stage()`

```
_I_n_t_e_r_p_o_l_a_t_e _S_t_a_g_e _f_r_o_m _A_E_P

_D_e_s_c_r_i_p_t_i_o_n:

     Given a stage-frequency curve, interpolates stage values at target
     AEPs. Interpolation is performed in standard normal (z-variate)
     space for better behavior across orders of magnitude in AEP.

_U_s_a_g_e:

     aep2stage(aep, stage, interp_aep)
     
_A_r_g_u_m_e_n_t_s:

     aep: Numeric vector of annual exceedance probabilities from the
          known curve.

   stage: Numeric vector of stage values corresponding to ‘aep’.

interp_aep: Numeric vector of target AEPs to interpolate at.

_V_a_l_u_e:

     Numeric vector of interpolated stage values (same length as
     ‘interp_aep’).

_S_e_e _A_l_s_o:

     ‘stage2aep’ for the inverse operation.

_E_x_a_m_p_l_e_s:

     # Get stage at the 1.12% and 0.175% AEP
     aep2stage(jmd_rfa_expected$AEP, jmd_rfa_expected$Expected, c(0.0112,1.75E-3))
     
```

### `stage2aep()`

```
_I_n_t_e_r_p_o_l_a_t_e _A_E_P _f_r_o_m _S_t_a_g_e

_D_e_s_c_r_i_p_t_i_o_n:

     Given a stage-frequency curve, interpolates AEP values at target
     stages. Interpolation is performed in standard normal (z-variate)
     space for better behavior across orders of magnitude in AEP.

_U_s_a_g_e:

     stage2aep(aep, stage, interp_stage)
     
_A_r_g_u_m_e_n_t_s:

     aep: Numeric vector of annual exceedance probabilities from the
          known curve.

   stage: Numeric vector of stage values corresponding to ‘aep’.

interp_stage: Numeric vector of target stage values to interpolate at.

_D_e_t_a_i_l_s:

     Stage-frequency curves typically flatten at very rare AEPs once
     the reservoir reaches the reservoir model discharge capacity. This
     produces tied stage values. ‘approx()’ handles these by averaging
     the corresponding AEPs, which may emit a "collapsing to unique 'x'
     values" warning. This is expected behavior.

_V_a_l_u_e:

     Numeric vector of interpolated AEP values (same length as
     ‘interp_stage’).

_S_e_e _A_l_s_o:

     ‘aep2stage’ for the inverse operation.

_E_x_a_m_p_l_e_s:

     # Get AEP at stage 3862.0-ft and 3872-ft
     stage2aep(jmd_rfa_expected$AEP, jmd_rfa_expected$Expected, c(3862.0, 3872))
     
```

### `interpolate_aep_matrix()`

```
_I_n_t_e_r_p_o_l_a_t_e _A_E_P_s _a_t _T_a_r_g_e_t _S_t_a_g_e_s _A_c_r_o_s_s _R_e_a_l_i_z_a_t_i_o_n_s

_D_e_s_c_r_i_p_t_i_o_n:

     For each realization's stage-frequency curve, interpolates the AEP
     at every target stage value. Interpolation is performed in
     standard normal (z-variate) space.

_U_s_a_g_e:

     interpolate_aep_matrix(results_list, target_stages)
     
_A_r_g_u_m_e_n_t_s:

results_list: List of realization data frames, each with columns
          ‘stage’ and ‘AEP’.

target_stages: Numeric vector of stage values to interpolate at.

_V_a_l_u_e:

     Matrix of AEP values with rows = target stages and columns =
     realizations.

```

### `interpolate_stage_matrix()`

```
_I_n_t_e_r_p_o_l_a_t_e _S_t_a_g_e_s _a_t _T_a_r_g_e_t _A_E_P_s _A_c_r_o_s_s _R_e_a_l_i_z_a_t_i_o_n_s

_D_e_s_c_r_i_p_t_i_o_n:

     For each realization's stage-frequency curve, interpolates the
     stage at every target AEP value. Interpolation is performed in
     standard normal (z-variate) space.

_U_s_a_g_e:

     interpolate_stage_matrix(results_list, target_aeps)
     
_A_r_g_u_m_e_n_t_s:

results_list: List of realization data frames, each with columns
          ‘stage’ and ‘AEP’.

target_aeps: Numeric vector of AEP values to interpolate at.

_V_a_l_u_e:

     Matrix of stage values with rows = target AEPs and columns =
     realizations.

```

### `power_function()`

```
_P_o_w_e_r _F_u_n_c_t_i_o_n

_D_e_s_c_r_i_p_t_i_o_n:

     Raises 10 to the 1st through the nth power

_U_s_a_g_e:

     power_function(x)
     
_A_r_g_u_m_e_n_t_s:

       x: Highest power of 10 to use

_V_a_l_u_e:

     Vector of 10 raised sequentially from the 1st to the nth power.

```

### `theme_rfar_conceptual()`

```
_r_f_a_R _g_g_p_l_o_t_2 _T_h_e_m_e _f_o_r _C_o_n_c_e_p_t_u_a_l _P_l_o_t_s

_D_e_s_c_r_i_p_t_i_o_n:

     A custom ‘ggplot2’ theme used in the rfaR conceptual realization
     vignette. Based on ‘theme_bw’ with tightened text sizes, italic
     titles, and the legend suppressed. Intended for figures that
     emphasize curves over annotation density.

_U_s_a_g_e:

     theme_rfar_conceptual()
     
_V_a_l_u_e:

     A ‘ggplot2’ theme object that can be added to a ‘ggplot’ via ‘+’.

_S_e_e _A_l_s_o:

     ‘theme_bw’

_E_x_a_m_p_l_e_s:

     library(ggplot2)
     ggplot(jmd_rfa_expected, aes(x = AEP, y = Expected)) +
       geom_line() +
       scale_x_continuous(transform = c("log10", "reverse")) +
       labs(title = "JMD Expected Stage-Frequency",
            subtitle = "Demonstrating theme_rfar_conceptual()",
            x = "AEP", y = "Stage (ft)") +
       theme_rfar_conceptual()
     
```