Commit a22a1941 authored by Manuela's avatar Manuela

small adjustments

parent 2917e402
......@@ -7,6 +7,7 @@
^.*\.Rproj$
^\.Rproj\.user$
^.Renviron
^.fun_stoch_sim_wave_diff_margins$
R/todo.R
Package: PRSim
Type: Package
Title: Stochastic Simulation of Streamflow Time Series using Phase Randomization
Version: 1.2-1
Date: 2020-01-10
Version: 1.2-2
Date: 2020-05-27
Authors@R: c(person("Manuela", "Brunner", role = c("aut", "cre"),
email = "manuela.i.brunner@gmail.com",
comment = c(ORCID = "0000-0001-8824-877X")),
......@@ -16,12 +16,17 @@ Description: Provides a simulation framework to simulate streamflow time series
main characteristics as observed data. These characteristics include the distribution of daily
streamflow values and their temporal correlation as expressed by short- and long-range
dependence. The approach is based on the randomization of the phases of the Fourier
transform or the phases of the wavelet transform. The function prsim() is applicable to single site simulation and uses the Fourier transform.
The function prsim.wave() extends the approach to multiple sites and is based on the complex wavelet transform. We further use the flexible four-parameter Kappa distribution, which allows
for the extrapolation to yet unobserved low and high flows. Alternatively, the empirical or any other distribution can be used. A detailed description of
the simulation approach for single sites and an application example can be found
in <https://www.hydrol-earth-syst-sci.net/23/3175/2019/>. A detailed description and evaluation of the wavelet-based multi-site approach can be found in
<https://www.hydrol-earth-syst-sci-discuss.net/hess-2019-658/>.
transform or the phases of the wavelet transform. The function prsim() is applicable
to single site simulation and uses the Fourier transform.
The function prsim.wave() extends the approach to multiple sites
and is based on the complex wavelet transform.
We further use the flexible four-parameter Kappa distribution,
which allows for the extrapolation to yet unobserved low and high flows.
Alternatively, the empirical or any other distribution can be used.
A detailed description of the simulation approach for single sites
and an application example can be found in <https://www.hydrol-earth-syst-sci.net/23/3175/2019/>.
A detailed description and evaluation of the wavelet-based multi-site approach
can be found in <https://www.hydrol-earth-syst-sci-discuss.net/hess-2019-658/>.
URL: https://git.math.uzh.ch/reinhard.furrer/PRSim-devel
BugReports: https://git.math.uzh.ch/reinhard.furrer/PRSim-devel
License: GPL-3
......@@ -30,3 +35,4 @@ LazyData: true
Depends: R (>= 3.5.0), homtest, goftest, wavScalogram, splus2R
Suggests: lattice, ismev, evd, GB2
Imports: stats, methods
RoxygenNote: 7.0.2
This diff is collapsed.
# short demo
###===============================###===============================###
### Application of PRSim.wave for different (extreme value) distributions
###===============================###===============================###
### for empirical distribution
### (1) Empirical distribution
### does not allow for extrapolation to yet unobserved values
###===============================###===============================###
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="empirical")
### for kappa distribution
### (2) Kappa distribution (default)
### 4 parameters and flexible. Found to be performing best for
### a set of 671 nearly natural catchments in the United States (Brunner and Gilleland 2020)
###===============================###===============================###
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="kappa", GoFtest = "KS")
# out <- prsim.wave(data=runoff_multi_sites, number_sim=5, marginal="kappa",
# GoFtest = NULL,pars=NULL, p_val=NULL)
......@@ -10,22 +18,25 @@ out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="kappa", GoFte
# setwd("C:/Users/mbrunner/Documents/PRSim-devel/data")
# save(simulations_multi_sites,file='simulations_multi_sites.rda')
### for GEV distribution
### (3) GEV distribution
### 3 parameters, classical extreme value distribution (Coles 2001)
###===============================###===============================###
require("evd")
require("ismev")
rGEV <- function(n, theta) rgev(n, theta[1], theta[2], theta[3])
pGEV <- function(x, theta) pgev(x, theta[1], theta[2], theta[3])
GEV_fit <- function( xdat, ...) gev.fit( xdat, show=FALSE, ...)$mle
GEV_fit <- function( xdat, ...) gev.fit(xdat, show=FALSE, ...)$mle
### GEV
### simulate using GEV
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="GEV", GoFtest = "KS", n_par=3)
### load stochastically simulated time series
### visualize simulations
### store stochastically simulated time series
### station 2
sim <- out[[2]]$simulation
### plot example of simulated time series
par(mai=c(.9,.9,.1,.1))
### observed time series
plot(sim$timestamp[1:1000], sim$Qobs[1:1000], type="l",
......@@ -34,4 +45,129 @@ plot(sim$timestamp[1:1000], sim$Qobs[1:1000], type="l",
matlines(sim$timestamp[1:1000], sim[1:1000, grep("r", names(sim))],
lty=1, col="gray")
### compare distributions without outliers
### without outliers
boxplot(sim$Qobs,sim$r1,outline=F,col=c('black','grey'),names=c('Obs','Sim'))
### with outliers
boxplot(sim$Qobs,sim$r1,col=c('black','grey'),names=c('Obs','Sim'))
### (4) generalized Beta distribution of the second kind (Mielke 1974)
### 4 parameters: can lead to very extreme outliers, which are potentially implausible
###===============================###===============================###
require( "GB2")
rGB2 <- function(n, theta) rgb2(n, theta[1], theta[2], theta[3], theta[4])
pGB2 <- function(x, theta) pgb2(x, theta[1], theta[2], theta[3], theta[4])
GB2_fit <- function( xdat, ...) ml.gb2( xdat, ...)$opt1$par
### simulate using GB2
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="GB2", GoFtest = "KS", n_par=4)
### (5) Normal distribution (just as a reference, will not perform well because it is not extreme)
### 2 parameters, not flexible enough for daily flow
###===============================###===============================###
library(fitdistrplus)
rNORM <- function(n, theta) rnorm(n, theta[1], theta[2])
pNORM <- function(x, theta) pnorm(x, theta[1], theta[2])
NORM_fit <- function( xdat, ...) fitdistr( xdat, 'normal', show=FALSE, ...)$estimate
### run simulations using normal distribution
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="NORM", GoFtest = "KS", n_par=2)
### (6) generalized Gamma distribution (Prentice 1974)
### 3 parameters
### non-stable in optimization
###===============================###===============================###
library(flexsurv)
rGENGAM <- function(n, theta) rgengamma(n, theta[1], theta[2], theta [3])
pGENGAM <- function(x, theta) pgengamma(x, theta[1], theta[2], theta [3])
GENGAM_fit <- function( xdat, ...) fitdistr(xdat, dgengamma, start=list("mu"=0,"sigma"=0.9,"Q"=0.5))$estimate
### depending on your data, it may be necessary to experiment with different start values
### as optimization may produce some 'random' errors.
### As optimization is not very stable, this distribution is maybe not
### very useful for stochastic simulation.
### run simulations with generalized Gamma
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="GENGAM", GoFtest = "KS", n_par=3)
### (7) Burr type XII distribution (Burr 1942)
### 3 parameters
### according to my test simulations performing very poorly.
###===============================###===============================###
library(actuar)
rBURR <- function(n, theta) rburr(n, theta[1], theta[2], theta [3])
pBURR <- function(x, theta) pburr(x, theta[1], theta[2], theta [3])
BURR_fit <- function( xdat, ...) fitdist(xdat, 'burr', start = list(shape1 = 0.3, shape2 = 0.9, rate = 0.9),method='mge')$estimate
### using maximum gooness-of-fit estimation instead of maximum likelihood estimation
### because the latter is numerically unstable.
### run simulations with Burr type XII distribution
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="BURR", GoFtest = "KS", n_par=3)
### (8) Wakeby distribution (Hosking 1986)
### 5 parameters, very flexible
###===============================###===============================###
library(lmomco)
rWAK <- function(n, theta) rlmomco(n, list(type='wak',para=theta,source='pwarwak',ifail=0,ifailtext="Successful parameter estimation."))
pWAK <- function(x, theta) pdfwak(x, list(type='wak',para=theta,source='pwarwak',ifail=0,ifailtext="Successful parameter estimation."))
WAK_fit <- function( xdat, ...){
lmom_X <- lmoms(xdat)
### fit Wakeby distribution
wakeby_X <- parwak(lmom_X)$para
return(wakeby_X)
}
### run simulations with Wakeby distribution
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="WAK", GoFtest = "KS", n_par=5)
### (9) Generalized Pareto
###===============================###===============================###
### produces NA values
rGPD <- function(n, theta) rgpd(n, theta[1], theta[2], theta[3])
pGPD <- function(x, theta) pgpd(x, theta[1], theta[2], theta[3])
GPD_fit <- function( xdat, ...){
param <- gpd.fit( xdat, threshold=min(xdat))
theta <- c(param$threshold,param$mle)
return(theta)
}
### run simulations with Wakeby distribution
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="GPD", GoFtest = "KS", n_par=3)
### (10) Lognormal distribution
###===============================###===============================###
rLNORM <- function(n, theta) rlnorm(n, theta[1], theta[2])
pLNORM <- function(x, theta) plnorm(x, theta[1], theta[2])
LNORM_fit <- function( xdat, ...) fitdistr( xdat, 'lognormal', show=FALSE, ...)$estimate
### run simulations with Wakeby distribution
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="LNORM", GoFtest = "KS", n_par=2)
### (11) Generalized logistic
### 3 parameters
### not worth testing as it produced negative values
###===============================###===============================###
# library(glogis)
# rGLOGIS <- function(n, theta) rglogis(n, theta[1], theta[2], theta[3])
# pGLOGIS <- function(x, theta) plnorm(x, theta[1], theta[2], theta[3])
# GLOGIS_fit <- function(xdat, ...) glogisfit(xdat)$coefficients
#
# ### run simulations with Wakeby distribution
# out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="GLOGIS", GoFtest = "KS", n_par=3)
### (12) Pearson type III
###===============================###===============================###
library(FAdist)
rPEARSON <- function(n, theta) rgamma3(n, theta[1], theta[2], theta[3])
pPEARSON <- function(x, theta) pgamma3(x, theta[1], theta[2], theta[3])
PEARSON_fit <- function( xdat, ...) fitdist(xdat, 'gamma3',start=list(shape=0.9,scale=0.9,thres=mean(xdat)),method='mge')$estimate
### using maximum gooness-of-fit estimation instead of maximum likelihood estimation
### because the latter is numerically unstable.
### run simulations with Wakeby distribution
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="PEARSON", GoFtest = "KS", n_par=3)
......@@ -27,6 +27,8 @@ Maintainer: \packageMaintainer{PRSim}
}
\references{
Brunner, M. I., A. Bárdossy, and R. Furrer (2019). Technical note: Stochastic simulation of streamflow time series using phase randomization. Hydrology and Earth System Sciences, 23, 3175-3187, doi:10.5194/hess-23-3175-2019
Brunner, M. I., and E. Gilleland (2020), Stochastic simulation of streamflow and spatial extremes: a continuous, wavelet-based approach, Hydrology and Earth System Sciences Discussion, https://doi.org/10.5194/hess-2019-658, in review, 2020.
}
\keyword{ package }
\examples{
......
......@@ -46,8 +46,7 @@ When using the kappa distribution, the AD test can for certain values of the par
\item{p_val}{A vector containing the p-values of \code{ks.test} or \code{ad.test} applied to the daily detrended data (if \code{GoFtest} is not NULL)}
}
\references{
Brunner and Gilleland to be submitted
}
Brunner, M. I., and E. Gilleland (2020), Stochastic simulation of streamflow and spatial extremes: a continuous, wavelet-based approach, Hydrology and Earth System Sciences Discussion, https://doi.org/10.5194/hess-2019-658, in review, 2020.}
\author{
Manuela Brunner
}
......
......@@ -26,7 +26,7 @@ The actual discharge data can be ordered from
\url{http://www.bafu.admin.ch/wasser/13462/13494/15076/index}.
}
\references{
Brunner, Bardossy, Furrer (2019) Technical note: Stochastic simulation of streamflow time series using phase randomization. Submitted.
Brunner, M. I., A. Bárdossy, and R. Furrer (2019). Technical note: Stochastic simulation of streamflow time series using phase randomization. Hydrology and Earth System Sciences, 23, 3175-3187, doi:10.5194/hess-23-3175-2019
}
\examples{
data(runoff)
......
......@@ -26,7 +26,7 @@ The actual discharge data were downloaded from
\url{https://waterdata.usgs.gov/nwis}.
}
\references{
Brunner and Gilleland to be submitted.
Brunner, M. I., and E. Gilleland (2020), Stochastic simulation of streamflow and spatial extremes: a continuous, wavelet-based approach, Hydrology and Earth System Sciences Discussion, https://doi.org/10.5194/hess-2019-658, in review, 2020.
}
\examples{
data(runoff_multi_sites)
......
......@@ -35,7 +35,7 @@ The data has been generated with
(default values for all other arguments).
}
\references{
Brunner, Bardossy, Furrer (2019) Technical note: Stochastic simulation of streamflow time series using phase randomization. Submitted.
Brunner, M. I., A. Bárdossy, and R. Furrer (2019). Technical note: Stochastic simulation of streamflow time series using phase randomization. Hydrology and Earth System Sciences, 23, 3175-3187, doi:10.5194/hess-23-3175-2019
}
\examples{
data(simulations)
......
......@@ -33,7 +33,7 @@ The data has been generated with
(default values for all other arguments).
}
\references{
Brunner and Gilleland to be submitted
Brunner, M. I., and E. Gilleland (2020), Stochastic simulation of streamflow and spatial extremes: a continuous, wavelet-based approach, Hydrology and Earth System Sciences Discussion, https://doi.org/10.5194/hess-2019-658, in review, 2020.
}
\examples{
### greys
......
File added
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment