Commit 5bed910e authored by Reinhard Furrer's avatar Reinhard Furrer

Quite some changes!

parent 5b9ae716
Package: PRSim
Type: Package
Title: Stochastic Simulation of Streamflow Time Series using Phase Randomization
Version: 1.0
Date: 2019-03-10
Version: 1.1
Date: 2019-06-21
Authors@R: c(person("Manuela", "Brunner", role = c("aut", "cre"),
email = "manuela.brunner@wsl.ch",
comment = c(ORCID = "0000-0001-8824-877X")),
......@@ -19,12 +19,12 @@ Description: Provides a simulation framework to simulate streamflow time series
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 and an application example can be found
in https://www.hydrol-earth-syst-sci-discuss.net/hess-2019-142/.
in <https://www.hydrol-earth-syst-sci-discuss.net/hess-2019-142/>.
URL: https://git.math.uzh.ch/reinhard.furrer/PRSim-devel
BugReports: https://git.math.uzh.ch/reinhard.furrer/PRSim-devel
License: GPL-3
Encoding: UTF-8
LazyData: true
Depends: homtest, goftest, graphics, ismev, evd
Suggests: lattice, GB2
Depends: homtest, goftest
Suggests: lattice, ismev, evd, GB2
Imports: stats
......@@ -2,8 +2,8 @@ importFrom("homtest", "Lmoments", "rand.kappa", "par.kappa", "F.kappa")
importFrom("stats", "rnorm", "runif", "fft", "ks.test")
importFrom("graphics", "hist")
importFrom("goftest", "ad.test")
importFrom("ismev","gev.fit")
importFrom("evd","rgev","pgev")
importFrom("GB2","mlfit.gb2","rgb2","pgb2")
# importFrom("ismev","gev.fit")
# importFrom("evd","rgev","pgev")
# importFrom("GB2","mlfit.gb2","rgb2","pgb2")
export("prsim")
prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
marginal=c("kappa","empirical","gev"),n_par=4,
verbose=TRUE, marginalpar=TRUE, suppWarn=FALSE, GoFtest=c(NULL, "AD", "KS"), ...){
### Reinhard: Do we somehow need to include "others" under marginal. Or is this fine?
marginal=c("kappa","empirical"), n_par=4, marginalpar=TRUE,
GoFtest=NULL, verbose=TRUE, suppWarn=FALSE, ...) {
ifft <- function (x) fft(x, inverse = TRUE)/length(x)
if (verbose) cat(paste0("Detrending with (half-)length ",win_h_length,"...\n"))
marginal <- tolower(marginal)[1]
if (!(marginal %in% c("kappa","empirical"))) { # check if distributions exist
if (!is.character(marginal)) stop("'marginal' should be a character string.")
r_CDF <- get(paste0("r_",marginal), mode = "function")
p_CDF <- get(paste0("p_",marginal), mode = "function")
CDF_fit <- get(paste0(marginal,"_fit"), mode = "function")
### I had to change the function names to not yet existing function names. Otherwise, the functions get stuck in recursive loops.
### functions themselves are defines outside the prsim function because they have diff
}
ifft <- function (x) fft(x, inverse = TRUE)/length(x)
## start preparing arguments.
if (!is.null(GoFtest)) {
GoFtest <- toupper(GoFtest)[1]
if (!(GoFtest %in% c("AD","KS"))) stop("'GoFtest' should be either 'NULL', 'AD' or 'KS'.")
} else GoFtest <- "NULL"
marginal <- marginal[1] # take only the first element
if (!(marginal %in% c("kappa","empirical"))) { # check if distributions exist
if (!is.character(marginal)) stop("'marginal' should be a character string.")
rCDF <- get(paste0("r",marginal), mode = "function")
CDF_fit <- get(paste0(marginal,"_fit"), mode = "function")
if (GoFtest=="AD") pCDF <- get(paste0("p",marginal), mode = "function")
}
op <- options("warn")$warn
### input data needs to be of the format year (four digits), month (two digits), day (one digit), input discharge time series
......@@ -51,7 +48,10 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
data <- data[format(data$timestamp, "%m %d") != "02 29",]
if ((nrow( data) %% 365)>0) stop("No missing values allowed. Some days are missing.")
if (verbose) cat(paste0("Detrending with (half-)length ",win_h_length,"...\n"))
# ### a) Detrend
# plot(data$Qobs[1:1000])
# # data$smooth <- loess(c(1:length(data$Qobs))~data$Qobs,span=0.0001)$y
......@@ -113,9 +113,9 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
### test whether Kappa distribution can be fit
if (suppWarn) {
suppressWarnings( test <- try(par.kappa(ll[1],ll[2],ll[4],ll[5])) )
suppressWarnings( test <- try(par.kappa(ll[1],ll[2],ll[4],ll[5]), silent = TRUE) )
} else {
test <- try(par.kappa(ll[1],ll[2],ll[4],ll[5]))
test <- try(par.kappa(ll[1],ll[2],ll[4],ll[5]), silent = TRUE)
}
if(length(test)>1){
......@@ -138,12 +138,15 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
if (tolower(GoFtest)=="ks")
p_vals[d] <- ks.test(data_window, data_kap)$p.value ### kappa distribution not rejected at alpha=0.05
if (tolower(GoFtest)=="ad")
if(length(try(ad.test(data_window,F.kappa,xi=kap_par$xi,alfa=kap_par$alfa,k=kap_par$k,h=kap_par$h)))==1){
if (tolower(GoFtest)=="ad") {
try_ad_test <- try(ad.test(data_window,F.kappa,xi=kap_par$xi,alfa=kap_par$alfa,k=kap_par$k,h=kap_par$h), silent=TRUE)
if(length(try_ad_test)==1){
p_vals[d] <- NA
}else{
p_vals[d] <- ad.test(data_window, F.kappa, xi=kap_par$xi, alfa=kap_par$alfa, k=kap_par$k, h=kap_par$h)$p.value
p_vals[d] <- try_ad_test$p.value
}
}
} else{
if(d==1){
p_vals[d] <- NA
......@@ -181,17 +184,10 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
ids <- c(before,after)
### determine values in window around day i
data_window <- data$Qobs[which(data$index%in%ids)]
### fit generalized gamma distribution
# fitdistr(x=data_window,dgengamma.stacy,start=list("scale"=1,"k"=6,"d"=2)) ### package VGAM: fitting does not work
### use function from package felxsurv
# theta <- CDF.fit(xdat=data_window)
theta <- CDF_fit(xdat=data_window)
# do.call(paste(CDF,".fit",sep=""),list(xdat=data_window))
theta <- CDF_fit(xdat=data_window, ...)
### goodness of fit test
# data_random <- rCDF(n=length(data_window),theta=theta)
data_random <- r_CDF(n=length(data_window),theta=theta)
data_random <- rCDF(n=length(data_window), theta)
# density_gengam[[d]] <- density(data_gengam)
hist(data_window)
......@@ -200,7 +196,7 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
p_vals[d] <- ks.test(data_window,data_random)$p.value
}
if (tolower(GoFtest)=="ad"){
p_vals[d] <- ad.test(data_window,p_CDF,theta)$p.value
p_vals[d] <- ad.test(data_window,pCDF,theta)$p.value
}
### store parameters
par_day[d,] <- theta
......@@ -227,6 +223,7 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
### repeat stochastic simulation several times
### list for storing results
data_sim <- list()
if(verbose) cat(paste0("Starting ",number_sim," simulations:\n"))
for (r in 1:number_sim){
### random generation of new phase sequence.
......@@ -357,7 +354,7 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
if(marginal!="kappa" & marginal!="empirical"){
### use monthly distribution for backtransformation
### simulate random sample of size n from disribution
data_day$cdf <- r_CDF(n=length(data_day$Qobs),par_day[d,])
data_day$cdf <- rCDF(n=length(data_day$Qobs), par_day[d,])
data_day$rank <- rank(data_day$cdf)
data_new$rank <- rank(data_new$seasonal)
......@@ -386,17 +383,10 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
deseaonalized=data$des,
data_sim)
if (is.null(GoFtest)) {
if (GoFtest=="NULL") {
p_vals <- NULL
}
if(marginal=="kappa"){
if (marginalpar) { # also return intermediate results
return(list( simulation=data_stoch, pars=par_day, p_val=p_vals))
p_vals <- NULL
}
}
if(marginal != "empirical"){
if (marginalpar) { # also return intermediate results
return(list( simulation=data_stoch, pars=par_day, p_val=p_vals))
......
......@@ -2,35 +2,18 @@
data(runoff)
out <- prsim(data=runoff, number_sim=1, marginal="empirical")
out <- prsim(data=runoff, number_sim=1, marginal="kappa",GoFtest = "KS")
out <- prsim(data=runoff, number_sim=1, marginal="kappa", GoFtest = "KS")
### GEV distribution
r_gev <- function(n, theta){
rgev(n, theta[1], theta[2], theta[3])
}
p_gev <- function(x, theta){
pgev(x, theta[1], theta[2], theta[3])
}
gev_fit <- function( xdat, ...) {
gev.fit( xdat, ...)$mle
}
### generalized Pareto distribution of the second kind
r_gb2 <- function(n, theta){
rgb2(n, theta[1], theta[2], theta[3], theta[4])
}
p_gb2 <- function(x, theta){
pgb2(x, theta[1], theta[2], theta[3], theta[4])
}
gb2_fit <- function( xdat, ...) {
mlfit.gb2( xdat, ...)[[2]]$par
}
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
out <- prsim(data=runoff, number_sim=1, marginal="gev",GoFtest = "KS",n_par=3)
### Generalized Beta of the second kind
out <- prsim(data=runoff, number_sim=1, marginal="gb2",GoFtest = "KS",n_par=4)
out <- prsim(data=runoff, number_sim=1, marginal="GEV", GoFtest = "KS", n_par=3)
sim <- out$simulation
# p_val <- out$p_val
......
......@@ -6,21 +6,22 @@
Applies the algorithm to a single station
}
\usage{
prsim(data, station_id="Qobs", number_sim=1, win_h_length=15,marginal="kappa",n_par=4,
verbose=TRUE, marginalpar=TRUE, suppWarn=FALSE, GoFtest=FALSE, ...)
prsim(data, station_id="Qobs", number_sim=1, win_h_length=15,
marginal=c("kappa","empirical"), n_par=4, marginalpar=TRUE,
GoFtest=NULL, verbose=TRUE, suppWarn=FALSE, ...)
}
\arguments{
\item{data}{data frame containing the time indications and runoff of at least one station. See \sQuote{Details}.}
\item{station_id}{identifies the station in case several runoffs are present in \code{data}. See \sQuote{Details}.}
\item{number_sim}{number of simulations to be carried out.}
\item{marginal}{marginal distribution to be used for the backtransformation. Can be either "kappa", "empirical", or any type of CDF (see \sQuote{Details}). "kappa" uses the four-parameter kappa distribution for backtransformation, "empirical" uses the empirical distribution. CDF allows for specifying any distribution \sQuote{Examples}.}
\item{n_par}{number of parameters of the marginal distribution used}
\item{win_h_length}{(half-)length of moving window size.}
\item{marginal}{marginal distribution to be used for the backtransformation. Can be either \code{"kappa"}, \code{"empirical"}, or any type of CDF (see \sQuote{Details}). \code{"kappa"} uses the four-parameter kappa distribution for backtransformation, \code{"empirical"} uses the empirical distribution. CDF allows for specifying any distribution \sQuote{Examples}.}
\item{n_par}{number of parameters of the marginal distribution used}
\item{GoFtest}{If (non-null) a GoF test for daily data should be performed: \code{"KS"} performs a Kolmogorof-Smirnov test, and \code{"AD"} performs an Anderson-Darling test. see \sQuote{Details})}
\item{verbose}{logical. Should progress be reported?}
\item{marginalpar}{logical. Should the estimated parameters of the distribution used be returned?}
\item{suppWarn}{logical. See \sQuote{Details}.}
\item{GoFtest}{NULL, KS or AD. Should a goodness-of-fit test for daily data be performed and the p-values be passed back? NULL performs no test, KS performs a Kolmogorof-Smirnov test, and AD performs an Anderson-Darling test. If (non-null) a GoF test for daily data should be performed, which one (\code{"AD"}, \code{KS}, see \sQuote{Details})}
\item{...}{any other argument passed to the sub-function specifying the cdf for fitting and drawing. See \sQuote{Details} and \sQuote{Examples}. }
\item{...}{any other argument passed to the sub-function specifying the cdf for fitting. See \sQuote{Details} and \sQuote{Examples}. }
}
\details{
Time can be given with three columns named \code{"YYYY"}, \code{"MM"}, \code{"DD"}, or as in POSIXct format \code{YYYY-MM-DD}.
......@@ -31,7 +32,7 @@ Stations are identified by column name (default \code{"Qobs"}), or by column ind
The function \code{homtest::par.kappa} might issue quite a few warnings of type \code{In fn(par, ...) : value out of range in 'gammafn'}. The argument \code{suppWarn} allows to silence warnings for the specific function call via \code{suppressWarnings()}. Of course, a subsequent check via \code{warnings()} is recommended.
Alternative distributions can be specified by providing three functions: (1) a function fitting the parameters of a distributions and providing a vector of these parameters as output (CDF_fit), (2) a function simulating random numbers from this distribution (r_CDF), and (3) a function specifying the distribution (p_CDF). \sQuote{Examples} for the generalized beta for the second kind and for the Generalized Extreme Values (GEV) distribution.
Alternative distributions can be specified by providing three functions: (1) a function fitting the parameters of a distributions and providing a vector of these parameters as output (CDF_fit), (2) a function simulating random numbers from this distribution (rCDF), and (3) a function specifying the distribution (pCDF). See \sQuote{Examples} for the generalized beta for the second kind and for the Generalized Extreme Values (GEV) distribution.
When using the kappa distribution, the AD test can for certain values of the parameter h not be performed.
......@@ -43,56 +44,45 @@ 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, Bardossy, Furrer (2019) Technical note: Stochastic simulation of streamflow time series using phase randomization. under review. https://www.hydrol-earth-syst-sci-discuss.net/hess-2019-142/.
Brunner, Bardossy, Furrer (2019) Technical note: Stochastic simulation of streamflow time series using phase randomization. Under review. \url{https://www.hydrol-earth-syst-sci-discuss.net/hess-2019-142/}.
}
\author{
Manuela Brunner
}
%\note{
%% ~~further notes~~
%}
%% ~Make other sections like Warning with \section{Warning }{....} ~
%\seealso{
%% ~~objects to See Also as \code{\link{help}}, ~~~
%}
\seealso{
\code{ks.test}
}
\examples{
data( runoff)
out <- prsim( runoff[ runoff$YYYY<1980, ], "Qobs", 1, suppWarn=TRUE)
# warnings() # as a follow-up to `suppWarn=TRUE`
## Specifying particular CDF:
## example with generalized Beta distribution of the second kind
require( "GB2")
CDF <- "gb2"
r_gb2 <- function(n, theta){
rgb2(n, theta[1], theta[2], theta[3], theta[4])
}
p_gb2 <- function(x, theta){
pgb2(x, theta[1], theta[2], theta[3], theta[4])
}
gb2_fit <- function( xdat, ...) {
mlfit.gb2( xdat, ...)[[2]]$par
## Specifying particular CDFs:
## (1) example with the Generalized Extreme Value (GEV) distribution
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, ...)$mle
\dontrun{ # The following call requires 5 seconds to execute
out <- prsim( runoff[ runoff$YYYY<1978, ], "Qobs", 1,
marginal="GEV", n_par=3, verbose=FALSE, marginalpar=FALSE,
show=FALSE) # Supress 'gev.fit' output.
}
out <- prsim( runoff[ runoff$YYYY<1980, ], "Qobs", 1, suppWarn=TRUE,
marginal="gb2")
## example with the Generalized Extreme Value (GEV) distribution
### test with GEV distribution
CDF <- "gev"
r_gev <- function(n, theta){
rgev(n, theta[1], theta[2], theta[3])
}
p_gev <- function(x, theta){
pgev(x, theta[1], theta[2], theta[3])
}
gev_fit <- function( xdat, ...) {
gev.fit( xdat, ...)$mle
}
out <- prsim( runoff[ runoff$YYYY<1980, ], "Qobs", 1, suppWarn=TRUE,
marginal="gev",n_par=3)
## (2) example with generalized Beta distribution of the second kind
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
\dontrun{ # The following call requires half minute or so to execute. Some warnings are issued
out <- prsim( runoff[ runoff$YYYY<1987, ], "Qobs", 1, suppWarn=TRUE,
marginal="GB2")
}
}
\keyword{ts}
R Under development (unstable) (2019-03-24 r76263) -- "Unsuffered Consequences"
R version 3.6.0 (2019-04-26) -- "Planting of a Tree"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
......@@ -7,6 +7,8 @@ R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
......@@ -20,6 +22,7 @@ Type 'q()' to quit R.
> options(warn = 1)
> library('PRSim')
Loading required package: homtest
Loading required package: goftest
>
> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
> base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv')
......@@ -65,10 +68,48 @@ Starting 1 simulations:
Finished.
> # warnings() # as a follow-up to `suppWarn=TRUE`
>
> ## Specifying particular CDFs:
> ## (1) example with the Generalized Extreme Value (GEV) distribution
> require("evd")
Loading required package: evd
> require("ismev")
Loading required package: ismev
Loading required package: mgcv
Loading required package: nlme
This is mgcv 1.8-28. For overview type 'help("mgcv-package")'.
> 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, ...)$mle
>
> ## Not run:
> ##D # The following call requires 5 seconds to execute
> ##D out <- prsim( runoff[ runoff$YYYY<1978, ], "Qobs", 1,
> ##D marginal="GEV", n_par=3, verbose=FALSE, marginalpar=FALSE,
> ##D show=FALSE) # Supress 'gev.fit' output.
> ## End(Not run)
>
> ## (2) example with generalized Beta distribution of the second kind
> require( "GB2")
Loading required package: 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
>
> ## Not run:
> ##D # The following call requires half minute or so to execute. Some warnings are issued
> ##D out <- prsim( runoff[ runoff$YYYY<1987, ], "Qobs", 1, suppWarn=TRUE,
> ##D marginal="GB2")
> ##D
> ## End(Not run)
>
>
>
>
> cleanEx()
detaching ‘package:GB2’, ‘package:ismev’, ‘package:mgcv’,
‘package:nlme’, ‘package:evd’
> nameEx("runoff")
> ### * runoff
>
......@@ -133,7 +174,7 @@ Finished.
> cleanEx()
> options(digits = 7L)
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed: 1.34 0.032 1.377 0 0
Time elapsed: 2.156 0.072 2.227 0 0
> grDevices::dev.off()
null device
1
......
......@@ -15,7 +15,7 @@ unique(runoff$YYYY)
try( prsim( runoff[1:130,] )) # At least one year of data required.
try( prsim( runoff[1:730,] )) # At least one year of data required.
try( prsim( runoff[1:730,] )) # No missing values allowed. Some days are missing.
try( prsim( runoff[1:1445,] )) # No missing values allowed. Some days are missing.
try( prsim( runoff[runoff$YYYY<1976,] )) # At least one year of data required.
......@@ -26,26 +26,52 @@ suppressWarnings( out <- prsim( runoff[runoff$YYYY<1977,] ) )
runof <- runoff[runoff$YYYY<1980,]
set.seed(1)
out1 <- prsim( runof, kappapar=FALSE, suppWarn=TRUE)
str(out1 <- prsim( runof, marginalpar=FALSE, suppWarn=TRUE))
runo <- runof
names( runo) <- tolower( names(runof))
try( prsim( runo, kappapar=FALSE, suppWarn=TRUE)) # Wrong column for observations selected.
try( prsim( runo, marginalpar=FALSE, suppWarn=TRUE)) # Wrong column for observations selected.
runo <- runof[,4:1]
set.seed(1)
out3 <- prsim( runo, kappapar=FALSE, suppWarn=TRUE) # ok
out3 <- prsim( runo, marginalpar=FALSE, suppWarn=TRUE) # ok
identical(out1,out3)
runo <- runof[,4:1]
set.seed(1)
out4 <- prsim( runo, station_id=1, kappapar=FALSE, suppWarn=TRUE) # ok
out4 <- prsim( runo, station_id=1, marginalpar=FALSE, suppWarn=TRUE) # ok
identical(out1,out4)
tmp <- paste(runof$YYYY, runof$MM, runof$DD,sep=" ")
runo <- data.frame(time=as.POSIXct(strptime(tmp, format="%Y %m %d", tz="GMT")), Qobs=runof$Qobs)
set.seed(1)
out5 <- prsim( runo, kappapar=FALSE, suppWarn=TRUE) # ok
out5 <- prsim( runo, marginalpar=FALSE, suppWarn=TRUE) # ok
identical(out1,out5)
#
\ No newline at end of file
#
######################
# Test 'kappa' distribution with manual construction:
rKappa <- function(n, theta) rand.kappa(n, theta[1], theta[2], theta[3], theta[4])
Kappa_fit <- function(xdat, ...) {
ll <- Lmoments(xdat)
unlist(par.kappa(ll[1],ll[2],ll[4],ll[5]))
}
set.seed(1)
out6a <- prsim( runo, marginalpar=TRUE)
set.seed(1)
out6b <- prsim( runo, marginal="Kappa", marginalpar=TRUE)
identical(out6a$pars, out6b$pars) # columns are differently named...
colSums( (as.matrix(out6a$pars)-out6b$pars)^2)
summary(out6a$simulation-out6b$simulation)
plot(out6a$simulation$r1, type='l')
lines.default(out6b$simulation$r1, col=3)
rug( which(out6b$simulation$r1 != out6a$simulation$r1))
days_diff <- matrix(out6a$simulation$r1 != out6b$simulation$r1,nrow=365)
kap_par <- data.frame(out6a$pars)
thresh <- kap_par$xi + kap_par$alfa*(1 - kap_par$h^(-kap_par$k))/kap_par$k
image( cbind( is.na(thresh), days_diff))
R Under development (unstable) (2019-03-24 r76263) -- "Unsuffered Consequences"
R version 3.6.0 (2019-04-26) -- "Planting of a Tree"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
......@@ -21,6 +21,7 @@ Type 'q()' to quit R.
> require(PRSim)
Loading required package: PRSim
Loading required package: homtest
Loading required package: goftest
>
> # raw demos
> demo( "PRSim", ask=FALSE)
......@@ -32,14 +33,44 @@ Loading required package: homtest
> # short demo
> data(runoff)
> out <- prsim(runoff, 1, 10)
> out <- prsim(data=runoff, number_sim=1, marginal="empirical")
Detrending with (half-)length 15...
Starting 10 simulations:
..........
Starting 1 simulations:
.
Finished.
> out <- prsim(data=runoff, number_sim=1, marginal="kappa", GoFtest = "KS")
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
> ### GEV distribution
> require("evd")
Loading required package: evd
> require("ismev")
Loading required package: ismev
Loading required package: mgcv
Loading required package: nlme
This is mgcv 1.8-28. For overview type 'help("mgcv-package")'.
> 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
> out <- prsim(data=runoff, number_sim=1, marginal="GEV", GoFtest = "KS", n_par=3)
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
> sim <- out$simulation
> # p_val <- out$p_val
> par(mai=c(.9,.9,.1,.1))
> plot(sim$timestamp[1:1000], sim$Qobs[1:1000], type="l",
......@@ -180,18 +211,14 @@ Finished.
>
>
> try( prsim( runoff[1:130,] )) # At least one year of data required.
Detrending with (half-)length 15...
Error in prsim(runoff[1:130, ]) : At least one year of data required.
> try( prsim( runoff[1:730,] )) # At least one year of data required.
Detrending with (half-)length 15...
> try( prsim( runoff[1:730,] )) # No missing values allowed. Some days are missing.
Error in prsim(runoff[1:730, ]) :
No missing values allowed. Some days are missing.
> try( prsim( runoff[1:1445,] )) # No missing values allowed. Some days are missing.
Detrending with (half-)length 15...
Error in prsim(runoff[1:1445, ]) :
No missing values allowed. Some days are missing.
> try( prsim( runoff[runoff$YYYY<1976,] )) # At least one year of data required.
Detrending with (half-)length 15...
Error in prsim(runoff[runoff$YYYY < 1976, ]) :
At least one year of data required.
>
......@@ -206,22 +233,32 @@ Finished.
> runof <- runoff[runoff$YYYY<1980,]
>
> set.seed(1)
> out1 <- prsim( runof, kappapar=FALSE, suppWarn=TRUE)
> str(out1 <- prsim( runof, marginalpar=FALSE, suppWarn=TRUE))
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
List of 3
$ simulation:'data.frame': 1825 obs. of 7 variables:
..$ YYYY : int [1:1825] 1975 1975 1975 1975 1975 1975 1975 1975 1975 1975 ...
..$ MM : int [1:1825] 1 1 1 1 1 1 1 1 1 1 ...
..$ DD : int [1:1825] 1 2 3 4 5 6 7 8 9 10 ...
..$ timestamp : POSIXct[1:1825], format: "1975-01-01" "1975-01-02" ...
..$ Qobs : num [1:1825] 2.05 1.75 1.62 1.58 1.47 ...
..$ deseaonalized: num [1:1825] 1.595 0.738 1.512 0.944 0.919 ...
..$ r1 : num [1:1825] 0.812 0.665 0.449 0.842 0.928 ...
$ pars : NULL
$ p_val : NULL
>
> runo <- runof
> names( runo) <- tolower( names(runof))
> try( prsim( runo, kappapar=FALSE, suppWarn=TRUE)) # Wrong column for observations selected.
Detrending with (half-)length 15...
Error in prsim(runo, kappapar = FALSE, suppWarn = TRUE) :
> try( prsim( runo, marginalpar=FALSE, suppWarn=TRUE)) # Wrong column for observations selected.
Error in prsim(runo, marginalpar = FALSE, suppWarn = TRUE) :
Wrong column (name) for observations selected.
>
> runo <- runof[,4:1]
> set.seed(1)
> out3 <- prsim( runo, kappapar=FALSE, suppWarn=TRUE) # ok
> out3 <- prsim( runo, marginalpar=FALSE, suppWarn=TRUE) # ok
Detrending with (half-)length 15...
Starting 1 simulations:
.
......@@ -231,7 +268,7 @@ Finished.
>
> runo <- runof[,4:1]
> set.seed(1)
> out4 <- prsim( runo, station_id=1, kappapar=FALSE, suppWarn=TRUE) # ok
> out4 <- prsim( runo, station_id=1, marginalpar=FALSE, suppWarn=TRUE) # ok
Detrending with (half-)length 15...
Starting 1 simulations:
.
......@@ -242,7 +279,7 @@ Finished.
> tmp <- paste(runof$YYYY, runof$MM, runof$DD,sep=" ")
> runo <- data.frame(time=as.POSIXct(strptime(tmp, format="%Y %m %d", tz="GMT")), Qobs=runof$Qobs)
> set.seed(1)
> out5 <- prsim( runo, kappapar=FALSE, suppWarn=TRUE) # ok
> out5 <- prsim( runo, marginalpar=FALSE, suppWarn=TRUE) # ok
Detrending with (half-)length 15...
Starting 1 simulations:
.
......@@ -252,6 +289,57 @@ Finished.
>
> #
>
> ######################
> # Test 'kappa' distribution with manual construction:
> rKappa <- function(n, theta) rand.kappa(n, theta[1], theta[2], theta[3], theta[4])
> Kappa_fit <- function(xdat, ...) {
+ ll <- Lmoments(xdat)
+ unlist(par.kappa(ll[1],ll[2],ll[4],ll[5]))
+ }
> set.seed(1)
> out6a <- prsim( runo, marginalpar=TRUE)
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
> set.seed(1)
> out6b <- prsim( runo, marginal="Kappa", marginalpar=TRUE)
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
> identical(out6a$pars, out6b$pars) # columns are differently named...
[1] FALSE
> colSums( (as.matrix(out6a$pars)-out6b$pars)^2)
xi alfa k h
0 0 0 0
> summary(out6a$simulation-out6b$simulation)
YYYY MM DD timestamp Qobs
Min. :0 Min. :0 Min. :0 Length:1825 Min. :0
1st Qu.:0 1st Qu.:0 1st Qu.:0 Class :difftime 1st Qu.:0
Median :0 Median :0 Median :0 Mode :numeric Median :0
Mean :0 Mean :0 Mean :0 Mean :0
3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0
Max. :0 Max. :0 Max. :0 Max. :0
deseaonalized r1
Min. :0 Min. :-8.52690
1st Qu.:0 1st Qu.: 0.00000
Median :0 Median : 0.00000
Mean :0 Mean : 0.00043
3rd Qu.:0 3rd Qu.: 0.00000
Max. :0 Max. : 9.88258
>
> plot(out6a$simulation$r1, type='l')
> lines.default(out6b$simulation$r1, col=3)
> rug( which(out6b$simulation$r1 != out6a$simulation$r1))
>
> days_diff <- matrix(out6a$simulation$r1 != out6b$simulation$r1,nrow=365)
> kap_par <- data.frame(out6a$pars)
> thresh <- kap_par$xi + kap_par$alfa*(1 - kap_par$h^(-kap_par$k))/kap_par$k
> image( cbind( is.na(thresh), days_diff))
>
>
>
> proc.time()
user system elapsed
21.328 0.164 21.508
24.912 0.094 24.997
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