Skip to content
Snippets Groups Projects
Commit e006e776 authored by Manuela's avatar Manuela
Browse files

update output format

parent ad487d4b
No related branches found
No related tags found
No related merge requests found
......@@ -3,7 +3,7 @@
prsim.wave <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
marginal=c("kappa","empirical"), n_par=4, n_wave=100, marginalpar=TRUE,
GoFtest=NULL, out_dir, verbose=TRUE, suppWarn=FALSE, ...){
GoFtest=NULL, verbose=TRUE, suppWarn=FALSE, ...){
### function for backtransformation of continuous wavelet transform
### inverse wavelet transform
......@@ -423,28 +423,28 @@ prsim.wave <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
if (GoFtest=="NULL") {
p_vals <- NULL
}
if(!is.na(out_dir)){
### set output directory
setwd(out_dir)
if(marginal != "empirical"){
if (marginalpar) { # also return intermediate results
# return(list(simulation=data_stoch, pars=par_day, p_val=p_vals))
out_list <- list(simulation=data_stoch, pars=par_day, p_val=p_vals)
save(file=paste(l,'_stoch_sim.rda',sep=''),out_list)
} else {
# return(list(simulation=data_stoch, pars=NULL, p_val=p_vals))
out_list <- list(simulation=data_stoch, pars=NULL, p_val=p_vals)
save(file=paste(l,'_stoch_sim.rda',sep=''),out_list)
}
}else{
# return(list(simulation=data_stoch))
out_list <- list(simulation=data_stoch, pars=NULL, p_val=NULL)
save(file=paste(l,'_stoch_sim.rda',sep=''),simulation=out_list)
}
}
if(is.na(out_dir)){
# if(!is.na(out_dir)){
# ### set output directory
# setwd(out_dir)
# if(marginal != "empirical"){
# if (marginalpar) { # also return intermediate results
# # return(list(simulation=data_stoch, pars=par_day, p_val=p_vals))
# out_list <- list(simulation=data_stoch, pars=par_day, p_val=p_vals)
# save(file=paste(l,'_stoch_sim.rda',sep=''),out_list)
# } else {
# # return(list(simulation=data_stoch, pars=NULL, p_val=p_vals))
# out_list <- list(simulation=data_stoch, pars=NULL, p_val=p_vals)
# save(file=paste(l,'_stoch_sim.rda',sep=''),out_list)
# }
# }else{
# # return(list(simulation=data_stoch))
# out_list <- list(simulation=data_stoch, pars=NULL, p_val=NULL)
# save(file=paste(l,'_stoch_sim.rda',sep=''),simulation=out_list)
# }
# }
#
# if(is.na(out_dir)){
### store values in list
if(marginal != "empirical"){
if (marginalpar) { # also return intermediate results
......@@ -458,10 +458,10 @@ prsim.wave <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
# return(list(simulation=data_stoch))
out_list[[l]] <- list(simulation=data_stoch, pars=NULL, p_val=NULL)
}
}
# }
### to to next station
}
if(is.na(out_dir)){
# if(is.na(out_dir)){
return(out_list)
}
# }
}
dontrun{
# short demo
data(runoff_multi_sites)
### for empirical distribution
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="empirical",out_dir=NA)
### for kappa distribution
prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="kappa", GoFtest = "KS",out_dir=NA)
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="kappa", GoFtest = "KS",out_dir=NA)
### for GEV distribution
require("evd")
......@@ -16,7 +14,7 @@ GEV_fit <- function( xdat, ...) gev.fit( xdat, show=FALSE, ...)$mle
### GEV
prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="GEV", GoFtest = "KS", n_par=3, out_dir=NA)
out <- prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="GEV", GoFtest = "KS", n_par=3, out_dir=NA)
### load stochastically simulated time series
......@@ -30,5 +28,5 @@ plot(sim$timestamp[1:1000], sim$Qobs[1:1000], type="l",
### add simulations
matlines(sim$timestamp[1:1000], sim[1:1000, grep("r", names(sim))],
lty=1, col="gray")
}
......@@ -22,7 +22,6 @@ prsim.wave(data, station_id="Qobs", number_sim=1, win_h_length=15,
\item{verbose}{logical. Should progress be reported?}
\item{marginalpar}{logical. Should the estimated parameters of the distribution used be returned?}
\item{n_wave}{number of scales to be considered in the continuous wavelet transform.}
\item{out_dir}{path of output directory where the stochastic simulations for the stations considered are saved. If NA, list of data frames for the different sites is returned}
\item{suppWarn}{logical. See \sQuote{Details}.}
\item{...}{any other argument passed to the sub-function specifying the cdf for fitting. See \sQuote{Details} and \sQuote{Examples}. }
}
......@@ -58,8 +57,7 @@ Manuela Brunner
}
\examples{
data(runoff_multi_sites)
prsim.wave(runoff_multi_sites, "Qobs", 1, suppWarn=TRUE,
out_dir='~/PRSim-devel/data/simulations_multi_sites')
prsim.wave(runoff_multi_sites, "Qobs", 1, suppWarn=TRUE)
# warnings() # as a follow-up to `suppWarn=TRUE`
## Specifying particular CDFs:
......@@ -86,7 +84,7 @@ GB2_fit <- function( xdat, ...) ml.gb2( xdat, ...)$opt1$par
\dontrun{ # The following call requires half minute or so to execute.
Some warnings are issued
prsim.wave(runoff_multi_sites, "Qobs", 1, suppWarn=TRUE,
marginal="GB2", out_dir='dir_PRsim_wave_data/simulations_multi_sites')
marginal="GB2")
}
......
......@@ -43,15 +43,23 @@ 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')
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])
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')
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])
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment