Package 'nsRFA'

Non-supervised Regional Frequency Analysis

Description

The estimation of hydrological variables in ungauged basins is a very important topic for many purposes, from research to engineering applications and water management (see the PUB project, Sivapalan et al., 2003). Regardless of the method used to perform such estimation, the underlying idea is to transfer the hydrological information from gauged to ungauged sites. When observations of the same variable at different measuring sites are available and many data samples are used together as source of information, the methods are called regional methods. The well known regional frequency analysis (e.g. Hosking and Wallis, 1997), where the interest is in the assessment of the frequency of hydrological events, belong to this class of methods. In literature, the main studied variable is the flood peak and the most used regional approach is the index-flood method of Dalrymple (1960), in which it is implicitly assumed that the flood frequency distribution for different sites belonging to a homogeneous region is the same except for a site-specific scale factor, the index-flood (see Hosking and Wallis, 1997, for details). Hence, the estimation of the flood frequency distribution for an ungauged site can be divided into two parts: estimation of the index-flood (more in general, the index-value) through linear/non-linear relations with climatic and basin descriptors; estimation of the adimentional flood frequency distribution, the growth curve, assigning the ungauged basin to one homogeneous region.

nsRFA is a collection of statistical tools for objective (non-supervised) applications of the Regional Frequency Analysis methods in hydrology. This does not mean that Regional Frequency Analysis should be non-supervised. These tools are addressed to experts, to help their expert judgement. The functions in nsRFA allow the application of the index-flood method in the following points:

-	regionalization of the index-value;
-	formation of homogeneous regions for the growth curves;
-	fit of a distribution function to the empirical growth curve of each region;

Regionalization of the index-value

The index-value can be either the sample mean (e.g. Hosking and Wallis, 1997) or the sample median (e.g. Robson and Reed, 1999) or another scale parameter. Many methodological approaches are available for the index-value estimation, and their differences can be related to the amount of information available. Excluding direct methods, that use information provided by flow data available at the station of interest, regional estimation methods require ancillary hydrological and physical information. Those methods can be divided in two classes: the multiregressive approach and the hydrological simulation approach. For both of them, the best estimator is the one that optimizes some criterion, such as the minimum error, the minimum variance or the maximum efficiency. Due to its simplicity, the most frequently used method is the multiregressive approach (see e.g. Kottegoda & Rosso, 1998; Viglione et al., 2007a), that relates the index-flow to catchment characteristics, such as climatic indices, geologic and morphologic parameters, land cover type, etc., through linear or non-linear equations.

R provides the functions lm and nls for linear and non-linear regressions (package stats). With the package nsRFA, a tool to select the best linear regressions given a set of candidate descriptors, bestlm, is provided. In REGRDIAGNOSTICS several functions are provided to analyze the output of lm, such as: the coefficient of determination (classic and adjusted); the Student t significance test; the variance inflation factor (VIF); the root mean squared error (RMSE); the mean absolute error (MAE); the prediction intervals; predicted values by a jackknife (cross-validation) procedure. The functions in DIAGNOSTICS provide more general diagnostics of model results (that can be also non-linear), comparing estimated values with observed values.

More details are provided in vignettes:

nsRFA_ex01	How to use the package nsRFA (example 1):
	Regional frequency analysis of the annual flows in Piemonte
	and Valle d'Aosta

that can be accessed via vignette("nsRFA_ex01", package="nsRFA").

Formation of homogeneous regions for the growth curves

Different techniques exist, for example those that lead to the formation of fixed regions through cluster analysis (Hosking and Wallis, 1997, Viglione, 2007), or those based on the method of the region of influence (ROI, Burn, 1990). The regional procedure can be divided into two parts: the formation of regions and the assignment of an ungauged site to one of them. Concerning the first point, the sites are grouped according to their similarity in terms of those basin descriptors that are assumed to explain the shape of the growth curve. This shape is usually quantified in a parametric way. For instance, the coefficient of variation (CV) or the L-CV of the curve can be used for this purpose. The package nsRFA provide the functions in moments and Lmoments to calculate sample moments and L-moments. Subsequently, the selected parameter is related with basin descriptors through a linear or a more complex model. A regression analysis is performed with different combination of descriptors, and descriptors that are strongly related with the parameter are used to group sites in regions. The same tools used for the regionalization of the index value, i.e. bestlm, REGRDIAGNOSTICS and DIAGNOSTICS, can be used if the parametric method is chosen.

nsRFA also provide a non-parametric approach that considers the dimensionless growth curve as a whole (see, Viglione et al., 2006; Viglione, 2007). The multiregressive approach can still be used if we reason in terms of (dis)similarity between pairs of basins in the following way: (1) for each couple of stations, a dissimilarity index between non-dimensional curves is calculated using a quantitatively predefined metric, for example using the Anderson-Darling statistic (A2), and organising the distances in a matrix with AD.dist; (2) for each basin descriptor, the absolute value (or another metric) of the difference between its measure in two basins is used as distance between them, using dist of the package stats to obtain distance matrices; (4) a multiregressive approach (bestlm, lm) is applied considering the matrices as variables and the basin descriptors associated to the best regression are chosen; (5) the significance of the linear relationship between distance matrices is assessed through the Mantel test with mantel.lm.

In the suitable-descriptor's space, stations with similar descriptor values can be grouped into disjoint regions through a cluster analysis (using functions in traceWminim) or the ROI method can be used adapting a region to the ungauged basin (roi). In both cases, the homogeneity of the regions can be assessed with the functions in HOMTESTS, where the Hosking and Wallis heterogeneity measures (HW.tests, see Hosking and Wallis, 1997) and the Anderson-Darling homogeneity test (ADbootstrap.test, see Viglione et al., 2007b) are provided.

More details are provided in vignettes:

nsRFA_ex01	How to use the package nsRFA (example 1):
	Regional frequency analysis of the annual flows in Piemonte
	and Valle d'Aosta
nsRFA_ex02	How to use the package nsRFA (example 2):
	Region-Of-Influence approach, some FEH examples

that can be accessed via vignette("nsRFA_ex01", package="nsRFA").

Fit of a distribution function to the empirical growth curve of each region

Once an homogeneous region is defined, the empirical growth curves can be pooled together and a probability distribution can be fitted to the pooled sample. The choice of the best distribution can be assisted by a Model Selection Criteria with MSClaio2008 (see, Laio et al., 2008). The parameters of the selected distribution can be estimated using the method of moments (moment_estimation), L-moments (par.GEV, par.genpar, par.gamma, ...) or maximum-likelihood (MLlaio2004). Goodness-of-fit tests are also available: the Anderson-Darling goodness of fit test with GOFlaio2004 (Laio. 2004), and Monte-Carlo based tests with GOFmontecarlo. Confidence intervals for the fitted distribution can be calculated with a Markov Chain Monte Carlo algorithm, using BayesianMCMC.

More details are provided in vignettes:

nsRFA_ex01	How to use the package nsRFA (example 1):
	Regional frequency analysis of the annual flows in Piemonte
	and Valle d'Aosta
MSClaio2008	Model selection techniques for the frequency analysis
	of hydrological extremes: the MSClaio2008 R function

that can be accessed via vignette("nsRFA_ex01", package="nsRFA").

Other functions

varLmoments provides distribution-free unbiased estimators of the variances and covariances of sample L-moments, as described in Elamir and Seheult (2004).

More details are provided in vignettes:

Fig1ElamirSeheult	Figure 1 in Elamir and Seheult (2004)

Details

Package:	nsRFA
Version:	0.7

The package provides several tools for Regional Frequency Analysis of hydrological variables. The first version dates to 2006 and was developed in Turin at the Politecnico by Alberto Viglione.

For a complete list of the functions, use library(help="nsRFA").

Main changes in version 0.7

0.7-17:	removal of old Fortran code and therefore of functions bestlm (and therefore the vignette nsRFA_ex01) and HW.original;
0.7-12:	refinement of function BayesianMCMC allowing several threshold and new function BayesianMCMCreg;
0.7-1:	refinement of function BayesianMCMC;
0.7-0:	plotting position for historical information in DISTPLOTS;

Main changes in version 0.6

0.6-9:	new vignette Fig11GriffisStedinger;
0.6-8:	exponential and Gumbel distributions added in GOFmontecarlo;
0.6-6:	some plotting position/probability plots have been added in DISTPLOTS;
0.6-4:	refinement of function BayesianMCMC;
0.6-2:	new vignette nsRFA_ex02;
0.6-2:	refinement of function BayesianMCMC;
0.6-0:	new vignette nsRFA_ex01;
0.6-0:	new function bestlm;
0.6-0:	the plotting position/probability plots in DISTPLOTS have been reshaped;
0.6-0:	this list of changes has been added;

Author(s)

Alberto Viglione

Maintainer: Alberto Viglione <[email protected]>

References

All the manual references are listed here:

Beirlant, J., Goegebeur, Y., Teugels, J., Segers, J., 2004. Statistics of Extremes: Theory and Applications. John Wiley and Sons Inc., 490 pp.

Burn, D.H., 1990. Evaluation of regional flood frequency analysis with a region of influence approach. Water Resources Research 26(10), 2257-2265.

Castellarin, A., Burn, D.H., Brath, A., 2001. Assessing the effectiveness of hydrological similarity measures for flood frequency analysis. Journal of Hydrology 241, 270-285.

Chowdhury, J.U., Stedinger, J.R., Lu, L.H., Jul. 1991. Goodness-of-fit tests for regional generalized extreme value flood distributions. Water Resources Research 27(7), 1765-1776.

D'Agostino, R., Stephens, M., 1986. Goodness-of-fit techniques. Vol.68 of Statistics: textBooks and monographs. Department of Statistics, Southern Methodist University, Dallas, Texas.

Dalrymple, T., 1960. Flood frequency analyses. Vol. 1543-A of Water Supply Paper. U.S. Geological Survey, Reston, Va.

Durbin, J., Knott, M., 1971. Components of Cramer-von Mises statistics. London School of Economics and Political Science, pp. 290-307.

El Adlouni, S., Bob\'ee, B., Ouarda, T.B.M.J., 2008. On the tails of extreme event distributions in hydrology. Journal of Hydrology 355(1-4), 16-33.

Elamir, E. A.H., Seheult, A.H., 2004. Exact variance structure of sample l-moments. Journal of Statistical Planning and Inference 124, 337-359.

Everitt, B., 1974. Cluster Analysis. Social Science Research Council. Halsted Press, New York.

Fill, H., Stedinger, J., 1998. Using regional regression within index flood procedures and an empirical bayesian estimator. Journal of Hydrology 210(1-4), 128-145.

Greenwood, J., Landwehr, J., Matalas, N., Wallis, J., 1979. Probability weighted moments: Definition and relation to parameters of several distributions expressible in inverse form. Water Resources Research 15, 1049-1054.

Hosking, J., 1986. The theory of probability weighted moments. Tech. Rep. RC12210, IBM Research, Yorktown Heights, NY.

Hosking, J., 1990. L-moments: analysis and estimation of distributions using linear combinations of order statistics. J. Royal Statistical Soc. 52, 105-124.

Hosking, J., Wallis, J., 1993. Some statistics useful in regional frequency analysis. Water Resources Research 29(2), 271-281.

Hosking, J., Wallis, J., 1997. Regional Frequency Analysis: An Approach Based on L-Moments. Cambridge University Press.

Institute of Hydrology, 1999. Flood Estimation Handbook, Institute of Hydrology, Oxford.

Kendall, M., Stuart, A., 1961-1979. The Advanced Theory of Statistics. Charles Griffin & Company Limited, London.

Kottegoda, N.T., Rosso, R., 1997. Statistics, Probability, and Reliability for Civil and Environmental Engineers, international Edition. McGraw-Hill Companies.

Laio, F., 2004. Cramer-von Mises and Anderson-Darling goodness of fit tests for extreme value distributions with unknown parameters, Water Resour. Res., 40, W09308, doi:10.1029/2004WR003204.

Laio, F., Tamea, S., 2007. Verification tools for probabilistic forecast of continuous hydrological variables. Hydrology and Earth System Sciences 11, 1267-1277.

Laio, F., Di Baldassarre G., Montanari A., 2008. Model selection techniques for the frequency analysis of hydrological extremes, Water Resour. Res., Under Revision.

Regione Piemonte, 2004. Piano di tutela delle acque. Direzione Pianificazione Risorse Idriche.

Robson, A., Reed, D., 1999. Statistical procedures for flood frequency estimation. In: Flood Estimation HandBook. Vol.~3. Institute of Hydrology Crowmarsh Gifford, Wallingford, Oxfordshire.

Sankarasubramanian, A., Srinivasan, K., 1999. Investigation and comparison of sampling properties of l-moments and conventional moments. Journal of Hydrology 218, 13-34.

Scholz, F., Stephens, M., 1987. K-sample Anderson-Darling tests. Journal of American Statistical Association, 82(399), pp. 918-924.

Sivapalan, M., Takeuchi, K., Franks, S.W., Gupta, V.K., Karambiri, H., Lakshmi, V., Liang, X., McDonnell, J.J., Mendiondo, E.M., O'Connell, P.E., Oki, T., Pomeroy, J.W, Schertzer, D., Uhlenbrook, S., Zehe, E., 2003. IAHS decade on Predictions in Ungauged Basins (PUB), 2003-2012: Shaping an exciting future for the hydrological sciences, Hydrological Sciences - Journal - des Sciences Hydrologiques, 48(6), 857-880.

Smouse, P.E., Long, J.C., Sokal, R.R., 1986. Multiple regression and correlation extensions of the Mantel test of matrix correspondence. Systematic Zoology, 35(4), 627-632.

Stedinger, J., Lu, L., 1995. Appraisal of regional and index flood quantile estimators. Stochastic Hydrology and Hydraulics 9(1), 49-75.

Stedinger, J.R., Vogel, R.M. and Foufoula-Georgiou, E., 1993. Frequency analysis of extreme events. In David R. Maidment, editor, Hand-Book of Hydrology, chapter 18. McGraw-Hill Companies, international edition.

Viglione, A., Claps, P., Laio, F., 2006. Utilizzo di criteri di prossimit\'a nell'analisi regionale del deflusso annuo, XXX Convegno di Idraulica e Costruzioni Idrauliche - IDRA 2006, Roma, 10-15 Settembre 2006.

Viglione, A., 2007. Metodi statistici non-supervised per la stima di grandezze idrologiche in siti non strumentati, PhD thesis at the Politecnico of Turin, February 2007.

Viglione, A., Claps, P., Laio, F., 2007a. Mean annual runoff estimation in North-Western Italy, In: G. La Loggia (Ed.) Water resources assessment and management under water scarcity scenarios, CDSU Publ. Milano.

Viglione, A., Laio, F., Claps, P., 2007b. A comparison of homogeneity tests for regional frequency analysis. Water Resources Research 43~(3).

Vogel, R., Fennessey, N., 1993. L moment diagrams should replace product moment diagrams. Water Resources Research 29(6), 1745-1752.

Vogel, R., Wilson, I., 1996. Probability distribution of annual maximum, mean, and minimum streamflows in the united states. Journal of Hydrologic Engineering 1(2), 69-76.

Ward, J., 1963. Hierarchical grouping to optimize an objective function, Journal of the American Statistical Association, 58, pp. 236-244.

---

Anderson-Darling distance matrix for growth curves

Description

Distance matrix for growth curves. Every element of the matrix is the Anderson-Darling statistic calculated between two series.

Usage

AD.dist (x, cod, index=2)

Arguments

x	vector representing data from many samples defined with cod
cod	array that defines the data subdivision among sites
index	if index=1 samples are divided by their average value; if index=2 (default) samples are divided by their median value

Details

The Anderson-Darling statistic used here is the one defined in https://en.wikipedia.org/wiki/Anderson-Darling_test as $A^2$.

Value

AD.dist returns the distance matrix between growth-curves built with the Anderson-Darling statistic.

Note

For information on the package and the Author, and for all the references, see nsRFA.

See Also

traceWminim, roi.

Examples

data(hydroSIMN)

annualflows
summary(annualflows)
x <- annualflows["dato"][,]
cod <- annualflows["cod"][,]

# Ad.dist
AD.dist(x,cod)             # it takes some time

---

Data-sample

Description

Systematic flood data and historical flood data for the Ard\'eche region (France) as described in: \ Naulet, R. (2002). Utilisation de l'information des crues historiques pour une meilleure pr\'ed\'etermination du risque d'inondation. Application au bassin de l'Ard\‘eche \'a Vallon Pont-d’Arc et St-Martin d'Ard\‘eche. PhD thesis at CEMAGREF, Lyon, and at the Universit\’e du Qu\'ebec, pp. 322. \ and \ Nguyen, C.C. (2012). Am\'elioration des approches Bay\'esiennes MCMC pour l'analyse r\'egionale des crues (Improvement of BayesianMCMC approaches for regional flood frequency analyses). PhD thesis at the University of Nantes, pp. 192.

Usage

data(Ardechedata)

Format

Ardeche_areas, areas (km2) of the gauged and ungauged catchments in the Ard\'eche region (France); Ardeche_ungauged_extremes, flood peaks (m3/s) reconstructed in ungauged catchments and number of years for which the peak was not exceeded; Beauvene_cont, sistematic flood peaks (m3/s) recorded at one station; Chambonas_cont, sistematic flood peaks (m3/s) recorded at one station; SaintLaurent_cont, sistematic flood peaks (m3/s) recorded at one station; SaintMartin_cont, sistematic flood peaks (m3/s) recorded at one station; SaintMartin_hist, values for historical peaks (m3/s) for one station and for flood perception thresholds (m3/s) non exceeded in the periods indicated; Vogue_cont, sistematic flood peaks (m3/s) recorded at one station.

Examples

data(Ardechedata)
SaintMartin_cont
SaintMartin_hist

---

Bayesian MCMC frequency analysis

Description

Bayesian Markov Chain Monte Carlo algorithm for flood frequency analysis with historical and other information. The user can choose between a local and a regional analysis.

Usage

BayesianMCMC (xcont, xhist=NA, infhist=NA, suphist=NA, 
               nbans=NA, seuil=NA, nbpas=1000, nbchaines=3, 
               confint=c(0.05, 0.95), dist="GEV",
               apriori=function(...){1}, 
               parameters0=NA, varparameters0=NA)
 BayesianMCMCcont (x, nbpas=NA)
 BayesianMCMCreg (xcont, scont, xhist=NA, infhist=NA, suphist=NA, shist=NA, 
                  nbans=NA, seuil=NA, nbpas=1000, nbchaines=3,
                  confint=c(0.05, 0.95), dist="GEV",
                  apriori=function(...){1},
                  parameters0=NA, varparameters0=NA) 
 BayesianMCMCregcont (x, nbpas=NA)
 plotBayesianMCMCreg_surf (x, surf, ask=FALSE, ...)
 ## S3 method for class 'BayesianMCMC'
 plot(x, which=1, ask=FALSE, ...)
 ## S3 method for class 'BayesianMCMC'
 print(x, ...)
 ## S3 method for class 'BayesianMCMCreg'
 plot(x, which=1, ask=FALSE, ...)
 ## S3 method for class 'BayesianMCMCreg'
 print(x, ...)

Arguments

x	object of class BayesianMCMC, output of function BayesianMCMC
xcont	vector of systematic data
scont	vector of upstream catchment surfaces of systematic data
xhist	vector of historical data and/or extreme discharges at ungauged sites
infhist	vector of inferior limit for historical data and/or extreme discharges at ungauged sites
suphist	vector of superior limit for historical data and/or extreme discharges at ungauged sites
shist	vector of upstream catchment surfaces of extreme discharges at ungauged sites
nbans	period (in years) over which every threshold has not been exceeded except for the historical data and/or extreme discharges at ungauged sites. If several values of xhist for a same threshold, put the number of years associated to the threshold on the first row, then put 0 (see examples)
seuil	threshold not exceeded in the historical period except for the historical data and/or extreme discharges at ungauged sites (several thresholds allowed).
nbpas	number of iterations for the MCMC algorithm
nbchaines	number of chains for the MCMC algorithm
confint	confidence limits for the flood quantiles
dist	distribution: normal "NORM", log-normal with 2 parameters "LN", Exponential "EXP", Gumbel "GUMBEL", Generalized Extreme Value "GEV", Generalized Logistic "GENLOGIS", Generalized Pareto "GENPAR", log-normal with 3 parameters "LN3", Pearson type III "P3", (log-Pearson type III "LP3", not implemented yet)
apriori	function of the parameters of the model ‘proportional to’ their a-priori guessed distribution. The default fuction returns always 1, i.e. there is no a-priori information
parameters0	initial values of the parameters for the MCMC algorithm
varparameters0	initial values of the parameter variances for the MCMC algorithm
which	a number of a vector of numbers that defines the graph to plot (see details)
ask	if TRUE, the interactive mode is run
surf	a particular surface (number or vector), not necessarily being a surface included in the scont or shist vectors
...	other arguments

Details

Supported cases

These functions are taking 4 cases into account, depending on the type of data provided: - Using only the systematic data: xcont provided, xhist=NA, infhist=NA and suphist=NA - Using censored information, historical flood known: xcont and xhist provided, infhist=NA and suphist=NA - Using censored information, historical flood unknown precisely but its lower limit known: xcont and infhist provided, xhist=NA and suphist=NA - Taking into account flood estimation intervals: infhist and suphist (respectively lower and upper limits) provided, xcont provided, xhist=NA - Please note that every other case is NOT supported. For example, you can't have some historical flood values perfectly known as well as some other for which you only know a lower limit or an interval.

Regarding the perception thresholds: - By definition, the number of exceedances of each perception threshold within its application period has to be known precisely, and all the floods exceeding the threshold have to be included in xhist (or infhist or suphist). - Several thresholds are allowed. - It is possible to include in xhist (or infhist or suphist) historical values that do not exceed the associated perception threshold. - If for one or several thresholds you only know that this or these threshold have never been exceeded and no more information is available on floods that did not exceed the threshold(s), this case is also supported. In this case, put for the historical flood corresponding to the threshold xhist=-1 (or infhist=-1 or infhist=suphist=-1).

Bayesian inference

Bayesian inference uses a numerical estimate of the degree of belief in a hypothesis before evidence has been observed and calculates a numerical estimate of the degree of belief in the hypothesis after evidence has been observed. The name ‘Bayesian’ comes from the frequent use of Bayes' theorem in the inference process. In our case the problem is: which is the probability that a frequency function $P$ (of type defined in dist) has parameters $\theta$, given that we have observed the realizations $D$ (defined in xcont, xhist, infhist, suphist, nbans, seuil). The Bayes' theorem writes

$P(\theta|D) = \frac{P(D|\theta) \cdot P(\theta)}{P(D)}$

where $P(\theta|D)$ is the conditional probability of $\theta$, given $D$ (it is also called the posterior probability because it is derived from or depends upon the specified value of $D$) and is the result we are interested in; $P(\theta)$ is the prior probability or marginal probability of $\theta$ (‘prior’ in the sense that it does not take into account any information about $D$), and can be given using the input apriori (it can be used to account for causal information); $P(D|\theta)$ is the conditional probability of $D$ given $\theta$ and it is defined choosing dist and depending on the availability of historical data; $P(D)$ is the prior or marginal probability of $D$, and acts as a normalizing constant. Intuitively, Bayes' theorem in this form describes the way in which one's beliefs about observing $\theta$ are updated by having observed $D$.

Since complex models cannot be processed in closed form by a Bayesian analysis, namely because of the extreme difficulty in computing the normalization factor $P(D)$, simulation-based Monte Carlo techniques as the MCMC approaches are used.

MCMC Metropolis algorithm

Markov chain Monte Carlo (MCMC) methods (which include random walk Monte Carlo methods), are a class of algorithms for sampling from probability distributions based on constructing a Markov chain that has the desired distribution as its equilibrium distribution. The state of the chain after a large number of steps is then used as a sample from the desired distribution. The quality of the sample improves as a function of the number of steps.

The MCMC is performed here through a simple Metropolis algorithm, i.e. a Metropolis-Hastings algorithm with symmetric proposal density. The Metropolis-Hastings algorithm can draw samples from any probability distribution $P(x)$, requiring only that a function proportional to the density can be calculated at $x$. In Bayesian applications, the normalization factor is often extremely difficult to compute, so the ability to generate a sample without knowing this constant of proportionality is a major virtue of the algorithm. The algorithm generates a Markov chain in which each state $x_t + 1$ depends only on the previous state $x_t$. The algorithm uses a Gaussian proposal density $N(x_t, \sigma_x)$, which depends on the current state $x_t$, to generate a new proposed sample $x'$. This proposal is accepted as the next value $x_t + 1 = x'$ if $\alpha$ drawn from $U(0,1)$ satisfies

$\alpha < \frac{P(x')}{P(x_t)}$

If the proposal is not accepted, then the current value of $x$ is retained ($x_t + 1 = x_t$).

The Markov chain is started from a random initial value $x_0$ and the algorithm is run for many iterations until this initial state is forgotten. These samples, which are discarded, are known as burn-in. The remaining set of accepted values of $x$ represent a sample from the distribution $P(x)$. As a Gaussian proposal density (or a lognormal one for definite-positive parameters) is used, the variance parameter $\sigma_x^2$ has to be tuned during the burn-in period. This is done by calculating the acceptance rate, which is the fraction of proposed samples that is accepted in a window of the last $N$ samples. The desired acceptance rate depends on the target distribution, however it has been shown theoretically that the ideal acceptance rate for a one dimensional Gaussian distribution is approx 50%, decreasing to approx 23% for an N-dimensional Gaussian target distribution. If $\sigma_x^2$ is too small the chain will mix slowly (i.e., the acceptance rate will be too high, so the sampling will move around the space slowly and converge slowly to $P(x)$). If $\sigma_x^2$ is too large the acceptance rate will be very low because the proposals are likely to land in regions of much lower probability density. The desired acceptance rate is fixed here to 34%.

The MCMC algorithm is based on a code developed by Eric Gaume on Scilab. It is still unstable and not all the distributions have been tested.

Value

BayesianMCMC and BayesianMCMCcont (which just continues the simulations of BayesianMCMC for local analyses and BayesianMCMCreg and BayesianMCMCregcont for regional analyses return the following values:

BayesianMCMCreg and BayesianMCMCregcont (which just continues the simulations of BayesianMCMC starting from its output) return the following values:

parameters matrix (nbpas)x(nbchaines) with the simulated sets of parameters with the MCMC algorithm;

parametersML set of parameters correspondent to the maximum likelihood;

returnperiods return periods for which quantilesML and intervals are calculated;

quantilesML quantiles correspondent to returnperiods for the distribution whose parameters are parametersML;

logML maximum log-likelihood;

intervals confidence intervals for the quantiles quantilesML for limits confint;

varparameters matrix (nbpas)x(nbchaines)x(number of parameters) with the simulated variances for the MCMC algorithm;

vraisdist likelihoods for the sets parameters;

propsaut vector showing the evolution of the acceptance rate during the Bayesian MCMC fitting;

plot.BayesianMCMC and plot.BayesianMCMCreg (for a normalized surface of 1 km2) plot the following figures:

1 data as plotting position (the Cunanne plotting position $a = 0.4$ is used), fitted distribution (maximum likelihood) and confidence intervals;

2 diagnostic plot of the MCMC simulation (parameters);

3 diagnostic plot of the MCMC simulation (likelihood and MCMC acceptance rate);

4 posterior distribution of parameters obtained with the MCMC simulation (cloud plots);

5 a-priori distribution of parameters (contour plots);

plotBayesianMCMCreg_surf plots the same plot as the first one given by plot.BayesianMCMCreg but for each surface in argument, as well as its mean as a function of the surfaces;

Note

For information on the package and the Author, and for all the references, see nsRFA.

Author(s)

Eric Gaume, Alberto Viglione, Jose Luis Salinas, Olivier Payrastre, Chi Cong N'guyen, Karine Halbert

See Also

.

Examples

set.seed(2988)
serie <- rand.GEV(120, xi=40, alfa=20, k=-0.4)
serie100 <- serie[1:100]
serie100[serie100 < 250] <- NA
serie20 <- serie[101:120]
serie <- c(serie100, serie20)


plot(serie, type="h", ylim=c(0, 600), xlab="", 
     ylab="Annual flood peaks [m3/s]", lwd=3)
abline(h=0)
points(serie100, col=2)


## Not run: 
# Using only sistematic data
only_sist <- BayesianMCMC (xcont=serie20, nbpas=5000, nbchaines=3, varparameters0=c(70, 20, 0.5), 
                           confint=c(0.05, 0.95), dist="GEV")
plot(only_sist, which=c(1:3), ask=TRUE, ylim=c(1,600))
only_sist <- BayesianMCMCcont(only_sist)
plot(only_sist, which=c(1:3), ask=TRUE, ylim=c(1,600))
only_sist <- BayesianMCMCcont(only_sist)
plot(only_sist, which=c(1:3), ask=TRUE, ylim=c(1,600))



# Adding the information that the threshold 250 m3/s was exceeded 
#   3 times in the past 100 years
with_hist_thresh <- BayesianMCMC (xcont=serie20, infhist=rep(250,3), 
                                  nbans=100, seuil=250,
                                  nbpas=5000, nbchaines=3, 
                                  confint=c(0.05, 0.95), dist="GEV")
plot(with_hist_thresh, which=c(1:3), ask=TRUE, ylim=c(1,600))



# Assuming that the 3 historical events are known with high uncertainty
with_hist_limits <- BayesianMCMC (xcont=serie20,  
                                  infhist=c(320,320,250), 
                                  suphist=c(360,400,270), 
                                  nbans=100, seuil=250,
                                  nbpas=5000, nbchaines=3, 
                                  confint=c(0.05, 0.95), dist="GEV")
plot(with_hist_limits, which=c(1:3), ask=TRUE, ylim=c(1,600))



# Assuming that the 3 historical events are perfectly known
with_hist_known <- BayesianMCMC (xcont=serie20, xhist=serie100[!is.na(serie100)], 
                                 nbans=100, seuil=250,
                                 nbpas=5000, nbchaines=3, 
                                 confint=c(0.05, 0.95), dist="GEV")
plot(with_hist_known, which=c(1:3), ask=TRUE, ylim=c(1,600))




# Perception threshold without available information on floods
without_info <- BayesianMCMC (xcont=serie20, xhist=-1, 
                                 nbans=100, seuil=2400,
                                 nbpas=5000, nbchaines=3, 
                                 confint=c(0.05, 0.95), dist="GEV")
plot(without_info, which=c(1:3), ask=TRUE, ylim=c(1,600))




# Using one reasonable a-priori distribution
fNORM3 <- function (x) {
 # x = vector of values
 # mu = vector of means
 mu = c(44, 26, -0.40)
 # CM = covariance matrix
 CM = matrix(c(13, 7.8, -0.055,
               7.8, 15, -0.42,
               -0.055, -0.42, 0.056), nrow=3, ncol=3)
 CMm1 <- solve(CM)
 term2 <- exp(-((x - mu) %*% CMm1 %*% (x - mu))/2)
 term1 <- 1/(2*pi)^(3/2)/sqrt(det(CM))
 term1*term2
}

with_hist_known2 <- BayesianMCMC (xcont=serie20, xhist=serie100[!is.na(serie100)], 
                                  nbans=100, seuil=250,
                                  nbpas=5000, nbchaines=3, apriori=fNORM3,
                                  confint=c(0.05, 0.95), dist="GEV")
plot(with_hist_known2, 5)
plot(with_hist_known2, 4)
plot(with_hist_known, 4)
plot(with_hist_known)
plot(with_hist_known2)



# Using one non-reasonable a-priori distribution
fNORM3 <- function (x) {
 # x = vector of values
 # mu = vector of means
 mu = c(30, 50, -0.10)
 # CM = covariance matrix
 CM = matrix(c(13, 7.8, -0.055,
               7.8, 15, -0.42,
               -0.055, -0.42, 0.056), nrow=3, ncol=3)
 CMm1 <- solve(CM)
 term2 <- exp(-((x - mu) %*% CMm1 %*% (x - mu))/2)
 term2
}

with_hist_known3 <- BayesianMCMC (xcont=serie20, xhist=serie100[!is.na(serie100)], 
                                  nbans=100, seuil=250,
                                  nbpas=5000, nbchaines=3, apriori=fNORM3,
                                  confint=c(0.05, 0.95), dist="GEV")
plot(with_hist_known3, 5)
plot(with_hist_known3, 4)
plot(with_hist_known, 4)
plot(with_hist_known)
plot(with_hist_known3)

## End(Not run)

## Not run: 
# Assuming that the historical events are perfectly known and there are 4 different thresholds 
# The data file is presenting this way: 

# xhist nbans seuil 
#  6000    55  6000 
#  7400    28  7250 
#  6350     8  3050 
#  4000     0  3050 
#  4550     0  3050 
#  3950     0  3050 
#  7550    58  2400 
#  4650     0  2400 
#  3950     0  2400 

## Warning: nbans and seuil should have the same length as xhist. 

# So when a threshold is exceeded several times, replicate it as many times it is exceeded 
# and part the number of years of exceedance into the number of times of exceedance 
# (the way you part the nbans vector is not important, what is important is that you have 
# length(nbans)=length(xhist) and the total of years for one same threshold equals the number 
# of years covered by the perception threshold) 
xhist_thres <- c(6000, 7400, 6350, 4000, 4550, 3950, 7550, 4650, 3950) 
seuil_thres <- c(6000, 7250, rep(3050, 4), rep(2400, 3)) 
nbans_thres <- c(55, 28, 8, 0, 0, 0, 58, 0, 0) 

# The threshold at 6000 has been exceeded for 55 years, the one at 7250 for 28 years, 
# the one at 3050 for 8 years and the one at 2400 for 58 years 
with_hist_known_several_thresholds <- BayesianMCMC (xcont=serie20, 
                                                    xhist=xhist_thres, 
                                                    nbans=nbans_thres, seuil=seuil_thres, 
                                                    nbpas=5000, nbchaines=3, 
                                                    confint=c(0.05, 0.95), dist="GEV", 
                                                    varparameters0=c(NA, NA, 0.5)) 
plot(with_hist_known_several_thresholds, which=c(1:3), ask=TRUE)


## REGIONAL:
# Regional analysis, assuming that the 3 historical events are perfectly known and 
# there are 2 perception thresholds
regional_with_hist_known <- BayesianMCMCreg (xcont=serie20, 
                                             scont=c(rep(507,9),rep(2240,11)), 
                                             xhist=serie100[!is.na(serie100)],
				             shist=c(495, 495, 87), 
                                             nbans=c(100, 0, 50), seuil=c(312, 312, 221),
                                             nbpas=5000, nbchaines=3, 
                                             confint=c(0.05, 0.95), dist="GEV", 
                                             varparameters0=c(NA, NA, NA, 0.5))
plot(regional_with_hist_known, which=1:3, ask=TRUE, ylim=c(1,600))

surf=c(571, 2240)
plotBayesianMCMCreg_surf(regional_with_hist_known, surf)

## End(Not run)

---

Diagnostics of models

Description

Diagnostics of model results, it compares estimated values y with observed values x.

Usage

R2 (x, y, na.rm=FALSE)
 RMSE (x, y, na.rm=FALSE) 
 MAE (x, y, na.rm=FALSE)
 RMSEP (x, y, na.rm=FALSE)
 MAEP (x, y, na.rm=FALSE)

Arguments

x	observed values
y	estimated values
na.rm	logical. Should missing values be removed?

Details

If $x_i$ are the observed values, $y_i$ the estimated values, with $i=1,...,n$, and $\bar{x}$ the sample mean of $x_i$, then:

$R^2 = 1 - \frac{\sum_1^n (x_i-y_i)^2}{\sum_1^n x_i^2 - n \bar{x}^2}$

$RMSE = \sqrt{\frac{1}{n} \sum_1^n (x_i-y_i)^2}$

$MAE = \frac{1}{n} \sum_1^n |x_i-y_i|$

$RMSEP = \sqrt{\frac{1}{n} \sum_1^n ((x_i-y_i)/x_i)^2}$

$MAEP = \frac{1}{n} \sum_1^n |(x_i-y_i)/x_i|$

See https://en.wikipedia.org/wiki/Coefficient_of_determination, https://en.wikipedia.org/wiki/Mean_squared_error and https://en.wikipedia.org/wiki/Mean_absolute_error for other details.

Value

R2 returns the coefficient of determination $R^2$ of a model.

RMSE returns the root mean squared error of a model.

MAE returns the mean absolute error of a model.

RMSE returns the percentual root mean squared error of a model.

MAE returns the percentual mean absolute error of a model.

Note

For information on the package and the Author, and for all the references, see nsRFA.

See Also

lm, summary.lm, predict.lm, REGRDIAGNOSTICS

Examples

data(hydroSIMN)

datregr <- parameters
regr0 <- lm(Dm ~ .,datregr); summary(regr0)
regr1 <- lm(Dm ~ Am + Hm + Ybar,datregr); summary(regr1)

obs <- parameters[,"Dm"]
est0 <- regr0$fitted.values
est1 <- regr1$fitted.values

R2(obs, est0)
R2(obs, est1)

RMSE(obs, est0)
RMSE(obs, est1)

MAE(obs, est0)
MAE(obs, est1)

RMSEP(obs, est0)
RMSEP(obs, est1)

MAEP(obs, est0)
MAEP(obs, est1)

---

Empirical distribution plots

Description

Sample values are plotted against their empirical distribution in graphs where points belonging to a particular distribution should lie on a straight line.

Usage

plotpos (x, a=0, orient="xF", ...)
 plotposRP (x, a=0, orient="xF", ...)
 loglogplot (x, a=0, orient="xF", ...)
 unifplot (x, a=0, orient="xF", line=FALSE, ...)
 normplot (x, a=0, orient="xF", line=FALSE, ...)
 lognormplot (x, a=0, orient="xF", line=FALSE, ...)
 studentplot (x, df, a=0, orient="xF", line=FALSE,...)
 logisplot (x, a=0, orient="xF", line=FALSE,...)
 gammaplot (x, shape, a=0, orient="xF", line=FALSE,...)
 expplot (x, a=0, orient="xF", line=FALSE,...)
 paretoplot (x, a=0, orient="xF", line=FALSE,...)
 gumbelplot (x, a=0, orient="xF", line=FALSE, ...)
 frechetplot (x, a=0, orient="xF", line=FALSE,...)
 weibullplot (x, a=0, orient="xF", line=FALSE,...)
 plotposRPhist (xcont, xhist=NA, infhist=NA, suphist=NA, nbans=NA, seuil=NA, 
                col12=c(1,1), a=0, orient="xF", ...)
 pointspos (x, a=0, orient="xF", ...)
 pointsposRP (x, a=0, orient="xF", ...)
 loglogpoints (x, a=0, orient="xF", ...)
 unifpoints (x, a=0, orient="xF", ...)
 normpoints (x, a=0, orient="xF", ...)
 studentpoints (x, df, a=0, orient="xF", ...)
 logispoints (x, a=0, orient="xF", ...)
 gammapoints (x, shape, a=0, orient="xF", ...)
 exppoints (x, a=0, orient="xF", ...)
 gumbelpoints (x, a=0, orient="xF", ...)
 weibullpoints (x, a=0, orient="xF", ...)
 regionalplotpos (x, cod, a=0, orient="xF", ...)
 regionalnormplot (x, cod, a=0, orient="xF", ...)
 regionallognormplot (x, cod, a=0, orient="xF", ...)
 regionalexpplot (x, cod, a=0, orient="xF", ...)
 regionalparetoplot (x, cod, a=0, orient="xF", ...)
 regionalgumbelplot (x, cod, a=0, orient="xF", ...)
 regionalfrechetplot (x, cod, a=0, orient="xF", ...)
 pointsposRPhist (xcont, xhist=NA, infhist=NA, suphist=NA, nbans=NA, seuil=NA, 
                  col12=c(1,1), a=0, orient="xF", ...)

Arguments

x	vector representing a data-sample
xcont	vector of systematic data (see BayesianMCMC)
xhist	vector of historical data (see BayesianMCMC)
infhist	inferior limit for historical data (see BayesianMCMC)
suphist	superior limit for historical data (see BayesianMCMC)
nbans	period (in years) over which the threshold has not been exceeded except for the historical data (see BayesianMCMC)
seuil	threshold non exceeded in the historical period except for the historical data (see BayesianMCMC)
df	degrees of freedom (> 0, maybe non-integer) of the Student t distribution. 'df = Inf' is allowed.
shape	shape parameter of the distribution
a	plotting position parameter, normally between 0 and 0.5 (the default value here, corresponding to the Hazen plotting position, see details)
orient	if orient="xF" the abscissa will be x and the ordinate F
line	if TRUE (default) a straight line indicating the normal, lognormal, ..., distribution with parameters estimated from x is plotted
cod	array that defines the data subdivision among sites
col12	vector of 2 elements containing the colors for the systematic and historical data respectively
...	graphical parameters as xlab, ylab, main, ...

Details

A brief introduction on Probability Plots (or Quantile-Quantile plots) is available on https://en.wikipedia.org/wiki/Q-Q_plot. For plotting positions see https://en.wikipedia.org/wiki/Plotting_position.

For the quantiles of the comparison distribution typically the Weibull formula $k/(n + 1)$ is used (default here). Several different formulas have been used or proposed as symmetrical plotting positions. Such formulas have the form

$(k - a)/(n + 1 - 2a)$

for some value of $a$ in the range from 0 to 1/2. The above expression $k/(n+1)$ is one example of these, for $a=0$. The Filliben plotting position has $a = 0.3175$ and the Cunanne plotting position has $a = 0.4$ should be nearly quantile-unbiased for a range of distributions. The Hazen plotting position, widely used by engineers, has $a = 0.5$. The Blom's plotting position, $a = 3/8$, gives nearly unbiased quantiles for the normal distribution, while the Gringeton plotting position, $a = 0.44$, is optimized for the largest observations from a Gumbel distribution. For the generalized Pareto, the GEV and related distributions of the Type I (Gumbel) and Weibull, $a = 0.35$ is suggested.

For large sample size, $n$, there is little difference between these various expressions.

Value

Representation of the values of x vs their empirical probability function $F$ in a cartesian, uniform, normal, lognormal or Gumbel plot. plotpos and unifplot are analogous except for the axis notation, unifplot has the same notation as normplot, lognormplot, ... plotposRP is analogous to plotpos but the frequencies $F$ are expressed as Return Periods $T=1/(1-F)$. With the default settings, $F$ is defined with the Weibull plotting position $F=k/(n+1)$. The straight line (if line=TRUE) indicate the uniform, normal, lognormal or Gumbel distribution with parameters estimated from x. The regional plots draw samples of a region on the same plot.

pointspos, normpoints, ... are the analogous of points, they can be used to add points or lines to plotpos, normplot, ... normpoints can be used either on normplot or lognormplot. exppoints can be used either on expplot or paretoplot (since the log-transformed Pareto random variable is exponentially distributed). gumbelpoints can be used either on gumbelplot or frechetplot (since the log-transformed Frechet random variable is distributed as a Gumbel).

loglogplot plots the logarithm of sample vs the logarithm of the empirical exceedance probability. For the log-log plot, the tail probability is represented by a straight line for power-law distributions (e.g. log-pearson, log-logistic, Frechet, ..., HEAVY TAIL), but not for the other subexponential or exponential distributions (e.g. gumbel, gamma, Pearson type III, ..., MODERATE TAIL); see El Adlouni et al. (2008).

plotposRPhist is based on the method in Stedinger et al. (1993, pp. 18.41-42).

Note

For information on the package and the Author, and for all the references, see nsRFA.

See Also

These functons are analogous to qqnorm; for the distributions, see Normal, Lognormal, LOGNORM, GUMBEL.

Examples

x <- rnorm(30,10,2)
plotpos(x)
normplot(x)
normplot(x,xlab=expression(D[m]),ylab=expression(hat(F)),
         main="Normal plot",cex.main=1,font.main=1)
normplot(x,line=FALSE)

x <- rlnorm(30,log(100),log(10))
normplot(x)
lognormplot(x)

x <- rand.gumb(30,1000,100)
normplot(x)
gumbelplot(x)

x <- rnorm(30,10,2)
y <- rnorm(50,10,3)
z <- c(x,y)
codz <- c(rep(1,30),rep(2,50))
regionalplotpos(z,codz)
regionalnormplot(z,codz,xlab="z")
regionallognormplot(z,codz)
regionalgumbelplot(z,codz)

plotpos(x)
pointspos(y,pch=2,col=2)

x <- rnorm(50,10,2)
F <- seq(0.01,0.99,by=0.01)
qq <- qnorm(F,10,2)
plotpos(x)
pointspos(qq,type="l")

normplot(x,line=FALSE)
normpoints(x,type="l",lty=2,col=3)

lognormplot(x)
normpoints(x,type="l",lty=2,col=3)

gumbelplot(x)
gumbelpoints(x,type="l",lty=2,col=3)

# distributions comparison in probabilistic graphs
x <- rnorm(50,10,2)
F <- seq(0.001,0.999,by=0.001)
loglikelhood <- function(param) {-sum(dgamma(x, shape=param[1], 
                scale=param[2], log=TRUE))}
parameters <- optim(c(1,1),loglikelhood)$par
qq <- qgamma(F,shape=parameters[1],scale=parameters[2])
plotpos(x)
pointspos(qq,type="l")

normplot(x,line=FALSE)
normpoints(qq,type="l")

lognormplot(x,line=FALSE)
normpoints(qq,type="l")

---

Two parameter exponential distribution and L-moments

Description

EXP provides the link between L-moments of a sample and the two parameter exponential distribution.

Usage

f.exp (x, xi, alfa)
F.exp (x, xi, alfa)
invF.exp (F, xi, alfa)
Lmom.exp (xi, alfa)
par.exp (lambda1, lambda2)
rand.exp (numerosita, xi, alfa)

Arguments

x	vector of quantiles
xi	vector of exp location parameters
alfa	vector of exp scale parameters
F	vector of probabilities
lambda1	vector of sample means
lambda2	vector of L-variances
numerosita	numeric value indicating the length of the vector to be generated

Details

See https://en.wikipedia.org/wiki/Exponential_distribution for a brief introduction on the Exponential distribution.

Definition

Parameters (2): $\xi$ (lower endpoint of the distribution), $\alpha$ (scale).

Range of $x$: $\xi \le x < \infty$.

Probability density function:

$f(x) = \alpha^{-1} \exp\{-(x-\xi)/\alpha\}$

Cumulative distribution function:

$F(x) = 1 - \exp\{-(x-\xi)/\alpha\}$

Quantile function:

$x(F) = \xi - \alpha \log(1-F)$

L-moments

$\lambda_1 = \xi + \alpha$

$\lambda_2 = 1/2 \cdot \alpha$

$\tau_3 = 1/3$

$\tau_4 = 1/6$

Parameters

If $\xi$ is known, $\alpha$ is given by $\alpha = \lambda_1 - \xi$ and the L-moment, moment, and maximum-likelihood estimators are identical. If $\xi$ is unknown, the parameters are given by

$\alpha = 2 \lambda_2$

$\xi = \lambda_1 - \alpha$

For estimation based on a single sample these estimates are inefficient, but in regional frequency analysis they can give reasonable estimates of upper-tail quantiles.

Lmom.exp and par.exp accept input as vectors of equal length. In f.exp, F.exp, invF.exp and rand.exp parameters (xi, alfa) must be atomic.

Value

f.exp gives the density $f$, F.exp gives the distribution function $F$, invFexp gives the quantile function $x$, Lmom.exp gives the L-moments ($\lambda_1$, $\lambda_2$, $\tau_3$, $\tau_4$), par.exp gives the parameters (xi, alfa), and rand.exp generates random deviates.

Note

For information on the package and the Author, and for all the references, see nsRFA.

See Also

rnorm, runif, GENLOGIS, GENPAR, GEV, GUMBEL, KAPPA, LOGNORM, P3; DISTPLOTS, GOFmontecarlo, Lmoments.

Examples

data(hydroSIMN)
annualflows
summary(annualflows)
x <- annualflows["dato"][,]
fac <- factor(annualflows["cod"][,])
split(x,fac)

camp <- split(x,fac)$"45"
ll <- Lmoments(camp)
parameters <- par.exp(ll[1],ll[2])
f.exp(1800,parameters$xi,parameters$alfa)
F.exp(1800,parameters$xi,parameters$alfa)
invF.exp(0.7870856,parameters$xi,parameters$alfa)
Lmom.exp(parameters$xi,parameters$alfa)
rand.exp(100,parameters$xi,parameters$alfa)

Rll <- regionalLmoments(x,fac); Rll
parameters <- par.exp(Rll[1],Rll[2])
Lmom.exp(parameters$xi,parameters$alfa)

---

Data-sample

Description

Flood Estimation Handbook flood peak data CD-ROM.

Usage

data(FEH1000)

Format

Data.frames:

am is the data.frame of the annual maximum flows with 3 columns: number, the code of the station; date, date of the annual maximum; year, year of the annual maximum (we consider hydrologic year: 1 october - 30 september); am, the value of the annual maximum flow [m3/s].

cd is the data.frame of parameters of 1000 catchements with 24 columns: number, the code of the station; nominal_area, catchment drainage area [km2]; nominal_ngr_x, basin outflow coordinates [m]; nominal_ngr_y, basin outflow coordinates [m]; ihdtm_ngr_x, basin outflow coordinates by Institute of Hydrology digital terrain model [m]; ihdtm_ngr_y, basin outflow coordinates by Institute of Hydrology digital terrain model [m]; dtm_area, catchment drainage area [km2] derived by CEH using their DTM (IHDTM); saar4170, standard average annual rainfall 1941-1970 [mm]; bfihost, baseflow index derived from HOST soils data; sprhost, standard percentage runoff derived from HOST soils data; farl, index of flood attenuation due to reservoirs and lakes; saar, standard average annual rainfall 1961-1990 [mm]; rmed_1d, median annual maximum 1-day rainfall [mm]; rmed_2d, median annual maximum 2-days rainfall [mm]; rmed_1h, median annual maximum 1-hour rainfall [mm]; smdbar, mean SMD (soil moisture deficit) for the period 1961-1990 calculated from MORECS month-end values [mm]; propwet, proportion of time when soil moisture deficit <=6 mm during 1961-90, defined using MORECS; ldp, longest drainage path [km], defined by recording the greatest distance from a catchment node to the defined outlet; dplbar, mean drainage path length [km]; altbar, mean catchment altitude [m]; dpsbar, mean catchement slope [m/km]; aspbar, index representing the dominant aspect of catchment slopes (its values increase clockwise from zero to 360, starting from the north). Mean direction of all inter-nodal slopes with north being zero; aspvar, index describing the invariability in aspect of catchment slopes. Values close to one when all slopes face a similar direction; urbext1990, extent of urban and suburban land cover in 1990 [fraction].

Note

For information on the package and the Author, and for all the references, see nsRFA.

Source

http://www.environment-agency.gov.uk/hiflowsuk/

Examples

data(FEH1000)
names(cd)
am[1:20,]

---

Data-sample

Description

Functions for inversion calculation.

Usage

data(functionsLaio)

Format

Data.frames:

df.kgev is a data.frame with the skewness coefficient (first column) and the corresponding shape parameter of the GEV (second column)

df.polig represents the poligamma function.

Note

For information on the package and the Author, and for all the references, see nsRFA.

---

Three parameter generalized logistic distribution and L-moments

Description

GENLOGIS provides the link between L-moments of a sample and the three parameter generalized logistic distribution.

Usage

f.genlogis (x, xi, alfa, k)
F.genlogis (x, xi, alfa, k)
invF.genlogis (F, xi, alfa, k)
Lmom.genlogis (xi, alfa, k)
par.genlogis (lambda1, lambda2, tau3)
rand.genlogis (numerosita, xi, alfa, k)

Arguments

x	vector of quantiles
xi	vector of genlogis location parameters
alfa	vector of genlogis scale parameters
k	vector of genlogis shape parameters
F	vector of probabilities
lambda1	vector of sample means
lambda2	vector of L-variances
tau3	vector of L-CA (or L-skewness)
numerosita	numeric value indicating the length of the vector to be generated

Details

See https://en.wikipedia.org/wiki/Logistic_distribution for an introduction to the Logistic Distribution.

Definition

Parameters (3): $\xi$ (location), $\alpha$ (scale), $k$ (shape).

Range of $x$: $-\infty < x \le \xi + \alpha / k$ if $k>0$; $-\infty < x < \infty$ if $k=0$; $\xi + \alpha / k \le x < \infty$ if $k<0$.

Probability density function:

$f(x) = \frac{\alpha^{-1} e^{-(1-k)y}}{(1+e^{-y})^2}$

where $y = -k^{-1}\log\{1 - k(x - \xi)/\alpha\}$ if $k \ne 0$, $y = (x-\xi)/\alpha$ if $k=0$.

Cumulative distribution function:

$F(x) = 1/(1+e^{-y})$

Quantile function: $x(F) = \xi + \alpha[1-\{(1-F)/F\}^k]/k$ if $k \ne 0$, $x(F) = \xi - \alpha \log\{(1-F)/F\}$ if $k=0$.

$k=0$ is the logistic distribution.

L-moments

L-moments are defined for $-1<k<1$.

$\lambda_1 = \xi + \alpha[1/k - \pi / \sin (k \pi)]$

$\lambda_2 = \alpha k \pi / \sin (k \pi)$

$\tau_3 = -k$

$\tau_4 = (1+5 k^2)/6$

Parameters

$k=-\tau_3$, $\alpha = \frac{\lambda_2 \sin (k \pi)}{k \pi}$, $\xi = \lambda_1 - \alpha (\frac{1}{k} - \frac{\pi}{\sin (k \pi)})$.

Lmom.genlogis and par.genlogis accept input as vectors of equal length. In f.genlogis, F.genlogis, invF.genlogis and rand.genlogis parameters (xi, alfa, k) must be atomic.

Value

f.genlogis gives the density $f$, F.genlogis gives the distribution function $F$, invF.genlogis gives the quantile function $x$, Lmom.genlogis gives the L-moments ($\lambda_1$, $\lambda_2$, $\tau_3$, $\tau_4$), par.genlogis gives the parameters (xi, alfa, k), and rand.genlogis generates random deviates.

Note

For information on the package and the Author, and for all the references, see nsRFA.

See Also

rnorm, runif, EXP, GENPAR, GEV, GUMBEL, KAPPA, LOGNORM, P3; DISTPLOTS, GOFmontecarlo, Lmoments.

Examples

data(hydroSIMN)
annualflows
summary(annualflows)
x <- annualflows["dato"][,]
fac <- factor(annualflows["cod"][,])
split(x,fac)

camp <- split(x,fac)$"45"
ll <- Lmoments(camp)
parameters <- par.genlogis(ll[1],ll[2],ll[4])
f.genlogis(1800,parameters$xi,parameters$alfa,parameters$k)
F.genlogis(1800,parameters$xi,parameters$alfa,parameters$k)
invF.genlogis(0.7697433,parameters$xi,parameters$alfa,parameters$k)
Lmom.genlogis(parameters$xi,parameters$alfa,parameters$k)
rand.genlogis(100,parameters$xi,parameters$alfa,parameters$k)

Rll <- regionalLmoments(x,fac); Rll
parameters <- par.genlogis(Rll[1],Rll[2],Rll[4])
Lmom.genlogis(parameters$xi,parameters$alfa,parameters$k)

---

Three parameter generalized Pareto distribution and L-moments

Description

GENPAR provides the link between L-moments of a sample and the three parameter generalized Pareto distribution.

Usage

f.genpar (x, xi, alfa, k)
F.genpar (x, xi, alfa, k)
invF.genpar (F, xi, alfa, k)
Lmom.genpar (xi, alfa, k)
par.genpar (lambda1, lambda2, tau3)
rand.genpar (numerosita, xi, alfa, k)

Arguments

x	vector of quantiles
xi	vector of genpar location parameters
alfa	vector of genpar scale parameters
k	vector of genpar shape parameters
F	vector of probabilities
lambda1	vector of sample means
lambda2	vector of L-variances
tau3	vector of L-CA (or L-skewness)
numerosita	numeric value indicating the length of the vector to be generated

Details

See https://en.wikipedia.org/wiki/Pareto_distribution for an introduction to the Pareto distribution.

Definition

Parameters (3): $\xi$ (location), $\alpha$ (scale), $k$ (shape).

Range of $x$: $\xi < x \le \xi + \alpha / k$ if $k>0$; $\xi \le x < \infty$ if $k \le 0$.

Probability density function:

$f(x) = \alpha^{-1} e^{-(1-k)y}$

where $y = -k^{-1}\log\{1 - k(x - \xi)/\alpha\}$ if $k \ne 0$, $y = (x-\xi)/\alpha$ if $k=0$.

Cumulative distribution function:

$F(x) = 1-e^{-y}$

Quantile function: $x(F) = \xi + \alpha[1-(1-F)^k]/k$ if $k \ne 0$, $x(F) = \xi - \alpha \log(1-F)$ if $k=0$.

$k=0$ is the exponential distribution; $k=1$ is the uniform distribution on the interval $\xi < x \le \xi + \alpha$.

L-moments

L-moments are defined for $k>-1$.

$\lambda_1 = \xi + \alpha/(1+k)]$

$\lambda_2 = \alpha/[(1+k)(2+k)]$

$\tau_3 = (1-k)/(3+k)$

$\tau_4 = (1-k)(2-k)/[(3+k)(4+k)]$

The relation between $\tau_3$ and $\tau_4$ is given by

$\tau_4 = \frac{\tau_3 (1 + 5 \tau_3)}{5+\tau_3}$

Parameters

If $\xi$ is known, $k=(\lambda_1 - \xi)/\lambda_2 - 2$ and $\alpha=(1+k)(\lambda_1 - \xi)$; if $\xi$ is unknown, $k=(1 - 3 \tau_3)/(1 + \tau_3)$, $\alpha=(1+k)(2+k)\lambda_2$ and $\xi=\lambda_1 - (2+k)\lambda_2$.

Lmom.genpar and par.genpar accept input as vectors of equal length. In f.genpar, F.genpar, invF.genpar and rand.genpar parameters (xi, alfa, k) must be atomic.

Value

f.genpar gives the density $f$, F.genpar gives the distribution function $F$, invF.genpar gives the quantile function $x$, Lmom.genpar gives the L-moments ($\lambda_1$, $\lambda_2$, $\tau_3$, $\tau_4$), par.genpar gives the parameters (xi, alfa, k), and rand.genpar generates random deviates.

Note

For information on the package and the Author, and for all the references, see nsRFA.

See Also

rnorm, runif, EXP, GENLOGIS, GEV, GUMBEL, KAPPA, LOGNORM, P3; DISTPLOTS, GOFmontecarlo, Lmoments.

Examples

data(hydroSIMN)
annualflows
summary(annualflows)
x <- annualflows["dato"][,]
fac <- factor(annualflows["cod"][,])
split(x,fac)

camp <- split(x,fac)$"45"
ll <- Lmoments(camp)
parameters <- par.genpar(ll[1],ll[2],ll[4])
f.genpar(1800,parameters$xi,parameters$alfa,parameters$k)
F.genpar(1800,parameters$xi,parameters$alfa,parameters$k)
invF.genpar(0.7161775,parameters$xi,parameters$alfa,parameters$k)
Lmom.genpar(parameters$xi,parameters$alfa,parameters$k)
rand.genpar(100,parameters$xi,parameters$alfa,parameters$k)

Rll <- regionalLmoments(x,fac); Rll
parameters <- par.genpar(Rll[1],Rll[2],Rll[4])
Lmom.genpar(parameters$xi,parameters$alfa,parameters$k)

---

Three parameter generalized extreme value distribution and L-moments

Description

GEV provides the link between L-moments of a sample and the three parameter generalized extreme value distribution.

Usage

f.GEV (x, xi, alfa, k)
F.GEV (x, xi, alfa, k)
invF.GEV (F, xi, alfa, k)
Lmom.GEV (xi, alfa, k)
par.GEV (lambda1, lambda2, tau3)
rand.GEV (numerosita, xi, alfa, k)

Arguments

x	vector of quantiles
xi	vector of GEV location parameters
alfa	vector of GEV scale parameters
k	vector of GEV shape parameters
F	vector of probabilities
lambda1	vector of sample means
lambda2	vector of L-variances
tau3	vector of L-CA (or L-skewness)
numerosita	numeric value indicating the length of the vector to be generated

Details

See https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution for an introduction to the GEV distribution.

Definition

Parameters (3): $\xi$ (location), $\alpha$ (scale), $k$ (shape).

Range of $x$: $-\infty < x \le \xi + \alpha / k$ if $k>0$; $-\infty < x < \infty$ if $k=0$; $\xi + \alpha / k \le x < \infty$ if $k<0$.

Probability density function:

$f(x) = \alpha^{-1} e^{-(1-k)y - e^{-y}}$

where $y = -k^{-1}\log\{1 - k(x - \xi)/\alpha\}$ if $k \ne 0$, $y = (x-\xi)/\alpha$ if $k=0$.

Cumulative distribution function:

$F(x) = e^{-e^{-y}}$

Quantile function: $x(F) = \xi + \alpha[1-(-\log F)^k]/k$ if $k \ne 0$, $x(F) = \xi - \alpha \log(-\log F)$ if $k=0$.

$k=0$ is the Gumbel distribution; $k=1$ is the reverse exponential distribution.

L-moments

L-moments are defined for $k>-1$.

$\lambda_1 = \xi + \alpha[1 - \Gamma (1+k)]/k$

$\lambda_2 = \alpha (1-2^{-k}) \Gamma (1+k)]/k$

$\tau_3 = 2(1-3^{-k})/(1-2^{-k})-3$

$\tau_4 = [5(1-4^{-k})-10(1-3^{-k})+6(1-2^{-k})]/(1-2^{-k})$

Here $\Gamma$ denote the gamma function

$\Gamma (x) = \int_0^{\infty} t^{x-1} e^{-t} dt$

Parameters

To estimate $k$, no explicit solution is possible, but the following approximation has accurancy better than $9 \times 10^{-4}$ for $-0.5 \le \tau_3 \le 0.5$:

$k \approx 7.8590 c + 2.9554 c^2$

where

$c = \frac{2}{3+\tau_3} - \frac{\log 2}{\log 3}$

The other parameters are then given by

$\alpha = \frac{\lambda_2 k}{(1-2^{-k})\Gamma(1+k)}$

$\xi = \lambda_1 - \alpha[1 - \Gamma(1+k)]/k$

Lmom.GEV and par.GEV accept input as vectors of equal length. In f.GEV, F.GEV, invF.GEV and rand.GEV parameters (xi, alfa, k) must be atomic.

Value

f.GEV gives the density $f$, F.GEV gives the distribution function $F$, invF.GEV gives the quantile function $x$, Lmom.GEV gives the L-moments ($\lambda_1$, $\lambda_2$, $\tau_3$, $\tau_4$), par.GEV gives the parameters (xi, alfa, k), and rand.GEV generates random deviates.

Note

For information on the package and the Author, and for all the references, see nsRFA.

See Also

rnorm, runif, EXP, GENLOGIS, GENPAR, GUMBEL, KAPPA, LOGNORM, P3; DISTPLOTS, GOFmontecarlo, Lmoments.

Examples

data(hydroSIMN)
annualflows
summary(annualflows)
x <- annualflows["dato"][,]
fac <- factor(annualflows["cod"][,])
split(x,fac)

camp <- split(x,fac)$"45"
ll <- Lmoments(camp)
parameters <- par.GEV(ll[1],ll[2],ll[4])
f.GEV(1800,parameters$xi,parameters$alfa,parameters$k)
F.GEV(1800,parameters$xi,parameters$alfa,parameters$k)
invF.GEV(0.7518357,parameters$xi,parameters$alfa,parameters$k)
Lmom.GEV(parameters$xi,parameters$alfa,parameters$k)
rand.GEV(100,parameters$xi,parameters$alfa,parameters$k)

Rll <- regionalLmoments(x,fac); Rll
parameters <- par.GEV(Rll[1],Rll[2],Rll[4])
Lmom.GEV(parameters$xi,parameters$alfa,parameters$k)

---

Goodness of fit tests

Description

Anderson-Darling goodness of fit tests for extreme-value distributions, from Laio (2004).

Usage

A2_GOFlaio (x, dist="NORM")
 A2 (F)
 W2 (F)
 fw2 (w)

Arguments

x	data sample
dist	distribution: normal "NORM", log-normal "LN", Gumbel "GUMBEL", Frechet "EV2", Generalized Extreme Value "GEV", Pearson type III "P3", log-Pearson type III "LP3"
F	cumulative distribution function (that has to be sorted increasingly)
w	Transformed test statistic (Laio, 2004)

Details

An introduction on the Anderson-Darling test is available on https://en.wikipedia.org/wiki/Anderson-Darling_test and in the GOFmontecarlo help page. The original paper of Laio (2004) is available on his web site.

Value

A2_GOFlaio tests the goodness of fit of a distribution with the sample x; it return the value $A_2$ of the Anderson-Darling statistics and its non-exceedence probability $P(A_2)$. Note that $P$ is the probability of obtaining the test statistic $A_2$ lower than the one that was actually observed, assuming that the null hypothesis is true, i.e., $P$ is one minus the p-value usually employed in statistical testing (see https://en.wikipedia.org/wiki/P-value). If $P(A_2)$ is, for example, greater than 0.90, the null hypothesis at significance level $\alpha=10\%$ is rejected.

A2 is the Anderson-Darling test statistic; it is used by A2_GOFlaio.

W2 is the Cramer-von Mises test statistic.

fw2 is the approximation of the probability distribution of w (first 2 terms) when $H_0$ is true (Anderson-Darling, 1952); it is used by A2_GOFlaio.

Note

For information on the package and the Author, and for all the references, see nsRFA.

See Also

GOFmontecarlo, MLlaio2004.

Examples

sm <- rand.gumb(100, 0, 1)
ml <- ML_estimation (sm, dist="GEV"); ml
F.GEV(sm, ml[1], ml[2], ml[3])
A2(sort(F.GEV(sm, ml[1], ml[2], ml[3])))
A2_GOFlaio(sm, dist="GEV")

ml <- ML_estimation (sm, dist="P3"); ml
A2(sort(sort(F.gamma(sm, ml[1], ml[2], ml[3]))))
A2_GOFlaio(sm, dist="P3")

---

Goodness of fit tests

Description

Anderson-Darling goodness of fit tests for Regional Frequency Analysis: Monte-Carlo method.

Usage

gofNORMtest (x)
 gofEXPtest (x, Nsim=1000)
 gofGUMBELtest (x, Nsim=1000)
 gofGENLOGIStest (x, Nsim=1000)
 gofGENPARtest (x, Nsim=1000)
 gofGEVtest (x, Nsim=1000)
 gofLOGNORMtest (x, Nsim=1000)
 gofP3test (x, Nsim=1000)

Arguments

x	data sample
Nsim	number of simulated samples from the hypothetical parent distribution

Details

An introduction, analogous to the following one, on the Anderson-Darling test is available on https://en.wikipedia.org/wiki/Anderson-Darling_test.

Given a sample $x_i \ (i=1,\ldots,m)$ of data extracted from a distribution $F_R(x)$, the test is used to check the null hypothesis $H_0 : F_R(x) = F(x,\theta)$, where $F(x,\theta)$ is the hypothetical distribution and $\theta$ is an array of parameters estimated from the sample $x_i$.

The Anderson-Darling goodness of fit test measures the departure between the hypothetical distribution $F(x,\theta)$ and the cumulative frequency function $F_m(x)$ defined as:

$F_m(x) = 0 \ , \ x < x_{(1)}$

$F_m(x) = i/m \ , \ x_{(i)} \leq x < x_{(i+1)}$

$F_m(x) = 1 \ , \ x_{(m)} \leq x$

where $x_{(i)}$ is the $i$-th element of the ordered sample (in increasing order).

The test statistic is:

$Q^2 = m \! \int_x \left[ F_m(x) - F(x,\theta) \right]^2 \Psi(x) \,dF(x)$

where $\Psi(x)$, in the case of the Anderson-Darling test (Laio, 2004), is $\Psi(x) = [F(x,\theta) (1 - F(x,\theta))]^{-1}$. In practice, the statistic is calculated as:

$A^2 = -m -\frac{1}{m} \sum_{i=1}^m \left\{ (2i-1)\ln[F(x_{(i)},\theta)] + (2m+1-2i)\ln[1 - F(x_{(i)},\theta)] \right\}$

The statistic $A^2$, obtained in this way, may be confronted with the population of the $A^2$'s that one obtain if samples effectively belongs to the $F(x,\theta)$ hypothetical distribution. In the case of the test of normality, this distribution is defined (see Laio, 2004). In other cases, e.g. the Pearson Type III case, can be derived with a Monte-Carlo procedure.

Value

gofNORMtest tests the goodness of fit of a normal (Gauss) distribution with the sample x.

gofEXPtest tests the goodness of fit of a exponential distribution with the sample x.

gofGUMBELtest tests the goodness of fit of a Gumbel (EV1) distribution with the sample x.

gofGENLOGIStest tests the goodness of fit of a Generalized Logistic distribution with the sample x.

gofGENPARtest tests the goodness of fit of a Generalized Pareto distribution with the sample x.

gofGEVtest tests the goodness of fit of a Generalized Extreme Value distribution with the sample x.

gofLOGNORMtest tests the goodness of fit of a 3 parameters Lognormal distribution with the sample x.

gofP3test tests the goodness of fit of a Pearson type III (gamma) distribution with the sample x.

They return the value $A_2$ of the Anderson-Darling statistics and its non exceedence probability $P$. Note that $P$ is the probability of obtaining the test statistic $A_2$ lower than the one that was actually observed, assuming that the null hypothesis is true, i.e., $P$ is one minus the p-value usually employed in statistical testing (see https://en.wikipedia.org/wiki/P-value). If $P(A_2)$ is, for example, greater than 0.90, the null hypothesis at significance level $\alpha=10\%$ is rejected.

Note

For information on the package and the Author, and for all the references, see nsRFA.

See Also

traceWminim, roi, HOMTESTS.

Examples

x <- rnorm(30,10,1)
gofNORMtest(x)

x <- rand.gamma(50, 100, 15, 7)
gofP3test(x, Nsim=200)

x <- rand.GEV(50, 0.907, 0.169, 0.0304)
gofGEVtest(x, Nsim=200)

x <- rand.genlogis(50, 0.907, 0.169, 0.0304)
gofGENLOGIStest(x, Nsim=200)

x <- rand.genpar(50, 0.716, 0.418, 0.476)
gofGENPARtest(x, Nsim=200)

x <- rand.lognorm(50, 0.716, 0.418, 0.476)
gofLOGNORMtest(x, Nsim=200)

---