Commit 0950b6f8 authored by Reinhard Furrer's avatar Reinhard Furrer

addression MLK

parent fee84d63
......@@ -6,6 +6,7 @@
^README\.md$
^.*\.Rproj$
^\.Rproj\.user$
^.Renviron
R/todo.R
R_LIBS=../lib
## Environmental variables to simulate R-devel CRAN errors
R_KEEP_PKG_SOURCE=yes
_R_S3_METHOD_LOOKUP_USE_TOPENV_AS_DEFENV_=true
_R_CHECK_INSTALL_DEPENDS_=true
_R_CHECK_SUGGESTS_ONLY_=false
_R_CHECK_NO_RECOMMENDED_=true
_R_CHECK_DOC_SIZES2_=true
_R_CHECK_DEPRECATED_DEFUNCT_=true
_R_CHECK_SCREEN_DEVICE_=warn
_R_CHECK_REPLACING_IMPORTS_=true
_R_CHECK_TOPLEVEL_FILES_=true
_R_CHECK_DOT_FIRSTLIB_=true
_R_CHECK_RD_LINE_WIDTHS_=true
_R_CHECK_S3_METHODS_NOT_REGISTERED_=true
_R_CHECK_OVERWRITE_REGISTERED_S3_METHODS_=true
_R_CHECK_CODE_USAGE_WITH_ONLY_BASE_ATTACHED_=TRUE
_R_CHECK_NATIVE_ROUTINE_REGISTRATION_=true
_R_CHECK_FF_CALLS_=registration
_R_CHECK_PRAGMAS_=true
_R_CHECK_COMPILATION_FLAGS_=true
_R_CHECK_R_DEPENDS_=true
_R_CHECK_PACKAGES_USED_IN_TESTS_USE_SUBDIRS_=true
_R_CHECK_SHLIB_OPENMP_FLAGS_=true
_R_CHECK_LIMIT_CORES_=true
_R_CHECK_LENGTH_1_CONDITION_=package:_R_CHECK_PACKAGE_NAME_,abort,verbose
_R_S3_METHOD_LOOKUP_BASEENV_AFTER_GLOBALENV_=true
#_R_CHECK_COMPILATION_FLAGS_KNOWN_="-Wno-deprecated-declarations -Wno-ignored-attributes -Wno-parentheses-Werror=format-security -Wp,-D_FORTIFY_SOURCE=2"
_R_CHECK_COMPILATION_FLAGS_KNOWN_="-Wp,-D_FORTIFY_SOURCE=2"
_R_CHECK_AUTOCONF_=true
Package: PRSim
Type: Package
Title: Stochastic Simulation of Streamflow Time Series using Phase Randomization
Version: 1.2
Date: 2019-11-20
Version: 1.2-1
Date: 2020-01-10
Authors@R: c(person("Manuela", "Brunner", role = c("aut", "cre"),
email = "manuela.i.brunner@gmail.com",
comment = c(ORCID = "0000-0001-8824-877X")),
......@@ -29,4 +29,3 @@ LazyData: true
Depends: R (>= 3.5.0), homtest, goftest, wmtsa
Suggests: lattice, ismev, evd, GB2
Imports: stats
RoxygenNote: 7.0.2
......@@ -16,4 +16,3 @@ BuildType: Package
PackageUseDevtools: Yes
PackageInstallArgs: --no-multiarch --with-keep.source
PackageCheckArgs: --as-cran
PackageRoxygenize: rd,collate,namespace
......@@ -137,7 +137,8 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
# hist(data_kap,add=T,col="red")
if (tolower(GoFtest)=="ks")
p_vals[d] <- ks.test(data_window, data_kap)$p.value ### kappa distribution not rejected at alpha=0.05
p_vals[d] <- ks_test(data_window, data_kap) ### kappa distribution not rejected at alpha=0.05
# p_vals[d] <- ks.test(data_window, data_kap)$p.value ### kappa distribution not rejected at alpha=0.05
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)
......@@ -193,7 +194,8 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
# hist(data_window)
# hist(data_random,add=T,col="red")
if (tolower(GoFtest)=="ks"){
p_vals[d] <- ks.test(data_window,data_random)$p.value
p_vals[d] <- ks_test(data_window,data_random)
# p_vals[d] <- ks.test(data_window,data_random)$p.value
}
if (tolower(GoFtest)=="ad"){
p_vals[d] <- ad.test(data_window,pCDF,theta)$p.value
......
# this is a stripped down version of the original stats::ks.test function, R Under development (unstable) (2019-12-10 r77548)
# we only use x,y (two.sided and exact=FALSE are default), pvalue is returned only.
# ks.test(x,y, exact=FALSE)$p.value = ks_test(x,y)
ks_test <- function (x, y)
{
x <- x[!is.na(x)]
n <- length(x)
y <- y[!is.na(y)]
n.x <- as.double(n)
n.y <- length(y)
n <- n.x * n.y/(n.x + n.y)
w <- c(x, y)
z <- cumsum(ifelse(order(w) <= n.x, 1/n.x, -1/n.y))
STATISTIC <- max(abs(z))
pkstwo <- function(x, tol = 1e-06) {
if (is.numeric(x))
x <- as.double(x)
else stop("argument 'x' must be numeric")
p <- rep(0, length(x))
p[is.na(x)] <- NA
IND <- which(!is.na(x) & (x > 0))
if (length(IND))
p[IND] <- .Call(stats:::C_pKS2, p = x[IND], tol)
p
}
PVAL <- 1 - pkstwo(sqrt(n) * STATISTIC)
return( min(1, max(0, PVAL)) )
}
R version 3.6.0 (2019-04-26) -- "Planting of a Tree"
R Under development (unstable) (2019-12-10 r77548) -- "Unsuffered Consequences"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
......@@ -23,6 +23,7 @@ Type 'q()' to quit R.
> library('PRSim')
Loading required package: homtest
Loading required package: goftest
Loading required package: wmtsa
>
> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
> base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv')
......@@ -43,6 +44,8 @@ Loading required package: goftest
> ## Not run:
> ##D demo("PRSim")
> ##D demo("PRSim-validate")
> ##D demo("PRSim_wave")
> ##D demo("PRSim_wave-validate")
> ## End(Not run)
>
>
......@@ -61,12 +64,10 @@ Loading required package: goftest
> ### ** Examples
>
> data( runoff)
> out <- prsim( runoff[ runoff$YYYY<1980, ], "Qobs", 1, suppWarn=TRUE)
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
> # warnings() # as a follow-up to `suppWarn=TRUE`
> ## Not run:
> ##D out <- prsim( runoff[ runoff$YYYY<1980, ], "Qobs", 1, suppWarn=TRUE)
> ##D # warnings() # as a follow-up to `suppWarn=TRUE`
> ## End(Not run)
>
> ## Specifying particular CDFs:
> ## (1) example with the Generalized Extreme Value (GEV) distribution
......@@ -76,7 +77,7 @@ Loading required package: evd
Loading required package: ismev
Loading required package: mgcv
Loading required package: nlme
This is mgcv 1.8-28. For overview type 'help("mgcv-package")'.
This is mgcv 1.8-31. 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
......@@ -107,6 +108,69 @@ Loading required package: GB2
>
> cleanEx()
detaching ‘package:GB2’, ‘package:ismev’, ‘package:mgcv’,
‘package:nlme’, ‘package:evd’
> nameEx("fun_stoch_sim_wave")
> ### * fun_stoch_sim_wave
>
> flush(stderr()); flush(stdout())
>
> ### Name: pRsim.wave
> ### Title: Simulate for multiple stations
> ### Aliases: PRsim.wave prsim.wave prsim_wave
> ### Keywords: ts
>
> ### ** Examples
>
> data(runoff_multi_sites)
> ## Not run:
> ##D # The following call requires half minute or so to execute.
> ##D prsim.wave(runoff_multi_sites, "Qobs", 1, suppWarn=TRUE)
> ## End(Not run)
> # 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-31. 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 prsim.wave(runoff_multi_sites, "Qobs", 1,
> ##D marginal="GEV", n_par=3, verbose=FALSE, marginalpar=FALSE,
> ##D show=FALSE)
> ##D # 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.
> ##D Some warnings are issued
> ##D prsim.wave(runoff_multi_sites, "Qobs", 1, suppWarn=TRUE,
> ##D marginal="GB2")
> ##D
> ## End(Not run)
>
>
>
>
> cleanEx()
detaching ‘package:GB2’, ‘package:ismev’, ‘package:mgcv’,
‘package:nlme’, ‘package:evd’
......@@ -138,6 +202,51 @@ detaching ‘package:GB2’, ‘package:ismev’, ‘package:mgcv’,
>
>
> cleanEx()
> nameEx("runoff_multi_sites")
> ### * runoff_multi_sites
>
> flush(stderr()); flush(stdout())
>
> ### Name: runoff_multi_sites
> ### Title: Sample runoff of four catchments with a similar discharge regime
> ### Aliases: runoff_multi_sites 'runoff multi sites'
> ### Keywords: datasets
>
> ### ** Examples
>
> data(runoff_multi_sites)
> str(runoff_multi_sites)
List of 4
$ :'data.frame': 14235 obs. of 4 variables:
..$ YYYY: Factor w/ 36 levels "1980","1984",..: 1 1 1 1 1 1 1 1 1 1 ...
..$ MM : Factor w/ 12 levels "01","02","03",..: 1 1 1 1 1 1 1 1 1 1 ...
..$ DD : Factor w/ 31 levels "01","02","03",..: 1 2 3 4 5 6 7 8 9 10 ...
..$ Qobs: num [1:14235] 2520 2460 2130 1650 1350 1120 976 892 814 748 ...
$ :'data.frame': 14235 obs. of 4 variables:
..$ YYYY: Factor w/ 39 levels "1980","1981",..: 1 1 1 1 1 1 1 1 1 1 ...
..$ MM : Factor w/ 12 levels "01","02","03",..: 1 1 1 1 1 1 1 1 1 1 ...
..$ DD : Factor w/ 31 levels "01","02","03",..: 1 2 3 4 5 6 7 8 9 10 ...
..$ Qobs: int [1:14235] 2970 2770 2250 1750 1520 1320 1210 1140 1080 1050 ...
$ :'data.frame': 14235 obs. of 4 variables:
..$ YYYY: Factor w/ 39 levels "1980","1981",..: 1 1 1 1 1 1 1 1 1 1 ...
..$ MM : Factor w/ 12 levels "01","02","03",..: 1 1 1 1 1 1 1 1 1 1 ...
..$ DD : Factor w/ 31 levels "01","02","03",..: 1 2 3 4 5 6 7 8 9 10 ...
..$ Qobs: num [1:14235] 5180 5690 5610 5150 5240 5190 4600 4700 7950 8550 ...
$ :'data.frame': 14235 obs. of 4 variables:
..$ YYYY: Factor w/ 39 levels "1980","1981",..: 1 1 1 1 1 1 1 1 1 1 ...
..$ MM : Factor w/ 12 levels "01","02","03",..: 1 1 1 1 1 1 1 1 1 1 ...
..$ DD : Factor w/ 31 levels "01","02","03",..: 1 2 3 4 5 6 7 8 9 10 ...
..$ Qobs: num [1:14235] 836 848 735 708 2210 1990 1220 944 854 944 ...
> runoff_multi_sites[[1]]$timestamp <- paste(runoff_multi_sites[[1]]$YYYY,
+ runoff_multi_sites[[1]]$MM, runoff_multi_sites[[1]]$DD, sep=" ")
> runoff_multi_sites[[1]]$timestamp <-
+ as.POSIXct(strptime(runoff_multi_sites[[1]]$timestamp,format="%Y %m %d", tz="GMT"))
> plot(runoff_multi_sites[[1]]$timestamp[1:1000], runoff_multi_sites[[1]]$Qobs[1:1000], type="l",
+ xlab="Time [d]", ylab=expression(paste("Discharge [m"^3,"/s]")))
>
>
>
> cleanEx()
> nameEx("simulations")
> ### * simulations
>
......@@ -169,12 +278,62 @@ detaching ‘package:GB2’, ‘package:ismev’, ‘package:mgcv’,
>
>
>
> cleanEx()
> nameEx("simulations_multi_sites")
> ### * simulations_multi_sites
>
> flush(stderr()); flush(stdout())
>
> ### Name: simulations_multi_sites
> ### Title: Simulated runoff for four catchments
> ### Aliases: simulations.multi.sites simulations_multi_sites
> ### Keywords: datasets
>
> ### ** Examples
>
> ### greys
> col_vect_obs <- c('#cccccc','#969696','#636363','#252525')
> ### oranges
> col_vect_sim <- c('#fdbe85','#fd8d3c','#e6550d','#a63603')
> data(simulations_multi_sites)
> sim <- simulations_multi_sites
> dim(sim[[1]])
[1] 14235 10
> ### plot time series for multiple sites
> par(mfrow=c(2,1),mar=c(3,3,2,1))
> ### determine ylim
> ylim_max <- max(sim[[1]]$Qobs)*1.5
> ### observed
> plot(sim[[1]]$Qobs[1:1000],
+ ylab=expression(bold(
+ paste("Specific discharge [mm/d]"))),
+ xlab="Time [d]",type="l",col=col_vect_obs[1],
+ ylim=c(0,ylim_max),main='Observations')
> for(l in 2:4){
+ lines(sim[[l]]$Qobs[1:1000],col=col_vect_obs[l])
+ }
> legend('topleft',legend=c('Station 1','Station 2',
+ 'Station 3','Station 4'),
+ lty=1,col=col_vect_obs[1:4])
> ### simulated (one run)
> plot(sim[[1]]$r1[1:1000],
+ ylab=expression(bold(paste("Specific discharge [mm/d]"))),
+ xlab="Time [d]",type="l",col=col_vect_sim[1],
+ ylim=c(0,ylim_max),
+ main='Stochastic simulations')
> for(l in 2:4){
+ lines(sim[[l]]$r1[1:1000],col=col_vect_sim[l])
+ }
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> ### * <FOOTER>
> ###
> cleanEx()
> options(digits = 7L)
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed: 2.156 0.072 2.227 0 0
Time elapsed: 1.963 0.093 2.062 0 0
> grDevices::dev.off()
null device
1
......
R version 3.6.0 (2019-04-26) -- "Planting of a Tree"
R Under development (unstable) (2019-12-10 r77548) -- "Unsuffered Consequences"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
......@@ -22,6 +22,7 @@ Type 'q()' to quit R.
Loading required package: PRSim
Loading required package: homtest
Loading required package: goftest
Loading required package: wmtsa
>
> # raw demos
> demo( "PRSim", ask=FALSE)
......@@ -53,7 +54,7 @@ Loading required package: evd
Loading required package: ismev
Loading required package: mgcv
Loading required package: nlme
This is mgcv 1.8-28. For overview type 'help("mgcv-package")'.
This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
> rGEV <- function(n, theta) rgev(n, theta[1], theta[2], theta[3])
......@@ -342,4 +343,4 @@ Finished.
>
> proc.time()
user system elapsed
24.912 0.094 24.997
25.623 0.072 25.725
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