Commit f0241bb5 authored by Reinhard Furrer's avatar Reinhard Furrer

extensions

parent 6d468b45
......@@ -25,6 +25,6 @@ BugReports: https://git.math.uzh.ch/reinhard.furrer/PRSim-devel
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")
importFrom("goftest", "ad.test")
export("prsim")
prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
marginal=c("kappa","empirical","CDF.fit"), draw="rCDF",
verbose=TRUE, marginalpar=TRUE, suppWarn=FALSE, GoFtest=c(NULL, "AD", "KS"), ...){
marginal=c("kappa","empirical"), 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.")
rCDF <- get(paste0("r",marginal), mode = "function")
pCDF <- get(paste0("p",marginal), mode = "function")
CDF.fit <- 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
......@@ -120,9 +134,9 @@ 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")
if (GoFtest=="KS")
p_vals[d] <- ks.test(data_window, data_kap)$p.value ### kappa distribution not rejected at alpha=0.05
if (GoFtest=="AD")
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 ### NOTE: TRY
} else{
......@@ -310,12 +324,12 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
deseaonalized=data$des,
data_sim)
if (!KStest) {
p_vals <- NULL
if (is.null(GoFtest)) {
p_vals <- NULL
}
if(marginal=="kappa"){
if (marginalpar) { # also return intermediate results
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))
......
new argument
marginal=c("kappa", "empirical", "gev.fit")
drawmarginal="rgev"
rCDF <- function(n, theta) {
rnorm(n, theta[1], sqrt(theta[2]))
......@@ -12,5 +9,5 @@ rCDF <- function(n, theta) {
CDF.fit <- function(xdat,...) { gev.fit(xdat,..., show=FALSE)$mle}
> data(portpirie)
> gev.fttt(portpirie[,2])
\ No newline at end of file
#data(portpirie)
#gev.fttt(portpirie[,2])
\ No newline at end of file
......@@ -6,19 +6,20 @@
Applies the algorithm to a single station
}
\usage{
prsim(data, station_id="Qobs", number_sim=1, win_h_length=15,marginal="kappa",
verbose=TRUE, kappapar=TRUE, suppWarn=FALSE, KStest=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, ...)
}
\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{win_h_length}{(half-)length of moving window size.}
\item{number_sim}{number of simulations to be carried out.}
\item{marginal}{marginal distribution to be used for the backtransformation.}
\item{win_h_length}{(half-)length of moving window size.}
\item{marginal}{marginal distribution to be used for the backtransformation, see \sQuote{Details}.}
\item{verbose}{logical. Should progress be reported?}
\item{marginalpar}{logical. Should the estimated parameters of the kappa distribution be returned?}
\item{suppWarn}{logical. See \sQuote{Details}.}
\item{KStest}{logical. Should a \code{ks.test} for daily data be performed and the p-values be passed back?}
\item{GoFtest}{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}. }
}
\details{
......@@ -33,6 +34,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