Commit 6df011ef authored by Manuela Brunner's avatar Manuela Brunner

Merge branch 'master' of

# Conflicts:
#	R/fun_stoch_sim.R
#	R/todo.R
#	man/fun_stoch_sim.Rd
parents 160dc3fe f0241bb5
......@@ -25,6 +25,6 @@ BugReports:
License: GPL-3
Encoding: UTF-8
LazyData: true
Depends: homtest
Suggests: lattice
Depends: homtest, goftest
Suggests: lattice, GB2
Imports: stats
importFrom("homtest", "Lmoments", "rand.kappa", "par.kappa")
importFrom("homtest", "Lmoments", "rand.kappa", "par.kappa", "F.kappa")
importFrom("stats", "rnorm", "runif", "fft", "ks.test")
<<<<<<< HEAD
importFrom("graphics", "hist")
importFrom("goftest", "ad.test")
>>>>>>> f0241bb562c45952cb5f9f249dbd9b2585c469b0
......@@ -6,6 +6,20 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
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.")
rCDF <- get(paste0("r",marginal), mode = "function")
pCDF <- get(paste0("p",marginal), mode = "function") <- get(paste0(marginal,".fit"), mode = "function")
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"
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
......@@ -122,6 +136,7 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
# density_kap[[d]] <- density(data_kap)
# hist(data_window)
# 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
if (tolower(GoFtest)=="ad")
......@@ -373,16 +388,23 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
if (is.null(GoFtest)) {
p_vals <- NULL
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=kap_par_day, p_val=p_vals))
} else {
return(list( simulation=data_stoch, pars=NULL, p_val=p_vals))
return(list( simulation=data_stoch) )
return(list(simulation=data_stoch) )
<<<<<<< HEAD
rCDF <- function(n, theta) {
rnorm(n, theta[1], sqrt(theta[2]))
} <- function(xdat,...) {,..., show=FALSE)$mle}
>>>>>>> f0241bb562c45952cb5f9f249dbd9b2585c469b0
......@@ -6,20 +6,34 @@
Applies the algorithm to a single station
<<<<<<< HEAD
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"), marginalpar=TRUE, GoFtest=NULL,
verbose=TRUE, suppWarn=FALSE, ...)
>>>>>>> f0241bb562c45952cb5f9f249dbd9b2585c469b0
\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{win_h_length}{(half-)length of moving window size.}
\item{number_sim}{number of simulations to be carried out.}
<<<<<<< HEAD
\item{marginal}{marginal distribution to be used for the backtransformation. Can be either "kappa", "empirical", or CDF. "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, see \sQuote{Details}.}
>>>>>>> f0241bb562c45952cb5f9f249dbd9b2585c469b0
\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}.}
<<<<<<< HEAD
\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.}
\item{GoFtest}{If (non-null) a GoF test for daily data should be performed, which one (\code{"AD"}, \code{KS}, see \sQuote{Details})}
>>>>>>> f0241bb562c45952cb5f9f249dbd9b2585c469b0
\item{...}{any other argument passed to the sub-function specifying the cdf for fitting and drawing. See \sQuote{Details} and \sQuote{Examples}. }
......@@ -34,6 +48,9 @@ The function \code{homtest::par.kappa} might issue quite a few warnings of type
To specify alternatives distr ... we need a function iwth input giving output and a funciton ... see \sQuote{Examples} for generalized beta for the second kind.
Add remark about AD test!
\value{A list with elements
\item{simulation }{A data frame with time information, observations, deseaonalized observations and
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