Commit bd14e72a authored by Reinhard Furrer's avatar Reinhard Furrer

some work...

parent 3b204ce0
......@@ -6,3 +6,6 @@
^README\.md$
^.*\.Rproj$
^\.Rproj\.user$
R/todo.R
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, KStest=FALSE, ...){
verbose=TRUE, marginalpar=TRUE, suppWarn=FALSE, GoFtest=c(NULL, "AD", "KS"), ...){
ifft <- function (x) fft(x, inverse = TRUE)/length(x)
......@@ -120,8 +120,11 @@ 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 (KStest)
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")
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{
if(d==1){
p_vals[d] <- NA
......@@ -148,8 +151,7 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
### Deseasonalization was not found to be necessary. Use normalized data directly.
data$des <- data$norm
# }
### 3) compute fast Fourier transform
ft <- fft(data$des)
......@@ -244,7 +246,7 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
# ### identify value corresponding to rank in the original time series
# data_ordered <- data_month[order(data_month$rank),]
# data_new$simulated_seasonal[which(data_new$MM%in%c(m))] <- data_ordered$Qobs[data_new$rank[which(data$MM%in%c(m))]]
# }
### use kappa distribution for backtransformation
if(marginal=="kappa"){
colnames(kap_par_day) <- names(kap_par)
......@@ -254,12 +256,12 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
data_day$kappa <- rand.kappa(length(data_day$Qobs),
xi=kap_par_day[d,"xi"],alfa=kap_par_day[d,"alfa"],
k=kap_par_day[d,"k"],h=kap_par_day[d,"h"])
# }
# if(marginal=="wakeby"){
# ### use Wakeby distribution for backtrasformation
# ### simulate random sample of size n from Wakeby distribution
# data_day$kappa <- rlmomco(length(data_day$Qobs),wak_par_day[[d]])
# }
#}
data_day$rank <- rank(data_day$kappa)
......@@ -272,7 +274,7 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
data_new$simulated_seasonal[which(data_new$index%in%c(d))] <- data_ordered$kappa[data_new$rank[which(data$index%in%c(d))]]
### if error was applied, replace negative values by 0 values
### in any case, replace negative values by 0. Corresponds to a bounded Kappa distribution
# if(error==TRUE){
if(length(which(data_new$simulated_seasonal<0))>0){
### do not use 0 as a replacement value directly
# data_new$simulated_seasonal[which(data_new$simulated_seasonal<0)] <- 0
......@@ -281,8 +283,7 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
rep_value <- runif(n=1,min=0,max=min(data_day$Qobs))
data_new$simulated_seasonal[which(data_new$simulated_seasonal<0)]<-rep_value
}
# }
}
}
### use empirical distribution for backtransformation
if(marginal=="empirical"){
data_day$rank <- rank(data_day$Qobs)
......@@ -294,7 +295,7 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
data_new$simulated_seasonal[which(data_new$index%in%c(d))] <- data_ordered$Qobs[data_new$rank[which(data$index%in%c(d))]]
# }
}
}
} # end for loop
data_sim[[r]] <- data_new$simulated_seasonal
if(verbose) cat(".")
......
new argument
marginal=c("kappa", "empirical", "gev.fit")
drawmarginal="rgev"
rCDF <- function(n, theta) {
rnorm(n, theta[1], sqrt(theta[2]))
}
CDF.fit <- function(xdat,...) { gev.fit(xdat,..., show=FALSE)$mle}
> data(portpirie)
> gev.fttt(portpirie[,2])
\ No newline at end of file
......@@ -7,7 +7,7 @@ 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)
verbose=TRUE, kappapar=TRUE, suppWarn=FALSE, KStest=FALSE, ...)
}
\arguments{
\item{data}{data frame containing the time indications and runoff of at least one station. See \sQuote{Details}.}
......@@ -16,9 +16,10 @@ prsim(data, station_id="Qobs", number_sim=1, win_h_length=15,marginal="kappa",
\item{number_sim}{number of simulations to be carried out.}
\item{marginal}{marginal distribution to be used for the backtransformation.}
\item{verbose}{logical. Should progress be reported?}
\item{kappapar}{logical. Should the estimated parameters of the kappa distribution be returned?}
\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{...}{any other argument passed to the sub-function specifying the cdf for fitting and drawing. 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}.
......@@ -30,11 +31,13 @@ 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.
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.
}
\value{A list with elements
\item{simulation }{A data frame with time information, observations, deseaonalized observations and
\code{number_sim} columns containing the simulated runoff.}
\item{kappa_pars}{A matrix containing the estimated parameters of the kappa distribution (if \code{kappapar}).}
\item{pars}{A matrix containing the estimated parameters of the marginal distribution (if \code{marginalpar}).}
\item{p_val}{A vector containing the p-values of \code{ks.test} applied to the daily detrended data (if \code{KStest}).}
}
\references{
......@@ -57,5 +60,21 @@ data( runoff)
out <- prsim( runoff[ runoff$YYYY<1980, ], "Qobs", 1, suppWarn=TRUE)
# warnings() # as a follow-up to `suppWarn=TRUE`
## Specifying particular CDF:
require( "GB2")
rCDF <- function(n, theta){
rgb2(n, theta[1], theta[2], theta[3], theta[4])
}
pCDF <- function(x, theta){
pgb2(x, theta[1], theta[2], theta[3], theta[4])
}
CDF.fit <- function( xdat, ...) {
mlfit.gb2( xdat, ...)[[2]]$par
}
out <- prsim( runoff[ runoff$YYYY<1980, ], "Qobs", 1, suppWarn=TRUE,
marginal="CDF.fit", draw="rCDF")
}
\keyword{ts}
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