| Title: | Bayesian Modeling: Estimate a Computer Model and Make Uncertain Predictions |
|---|---|
| Description: | An interface to the 'BaM' (Bayesian Modeling) engine, a 'Fortran'-based executable aimed at estimating a model with a Bayesian approach and using it for prediction, with a particular focus on uncertainty quantification. Classes are defined for the various building blocks of 'BaM' inference (model, data, error models, Markov Chain Monte Carlo (MCMC) samplers, predictions). The typical usage is as follows: (1) specify the model to be estimated; (2) specify the inference setting (dataset, parameters, error models...); (3) perform Bayesian-MCMC inference; (4) read, analyse and use MCMC samples; (5) perform prediction experiments. Technical details are available (in French) in Renard (2017) <https://hal.science/hal-02606929v1>. Examples of applications include Mansanarez et al. (2019) <doi:10.1029/2018WR023389>, Le Coz et al. (2021) <doi:10.1002/hyp.14169>, Perret et al. (2021) <doi:10.1029/2020WR027745>, Darienzo et al. (2021) <doi:10.1029/2020WR028607> and Perret et al. (2023) <doi:10.1061/JHEND8.HYENG-13101>. |
| Authors: | Benjamin Renard [aut, cre, cph] (ORCID: <https://orcid.org/0000-0001-8447-5430>), INRAE [fnd], Ministère de la Transition Ecologique - SCHAPI [fnd] |
| Maintainer: | Benjamin Renard <[email protected]> |
| License: | GPL-3 |
| Version: | 1.1.2 |
| Built: | 2026-06-08 07:32:21 UTC |
| Source: | https://github.com/bam-tools/rbam |
Run BaM.exe
BaM( workspace, mod, data, remnant = rep(list(remnantErrorModel()), mod$nY), mcmc = mcmcOptions(), cook = mcmcCooking(), summary = mcmcSummary(), residuals = residualOptions(), pred = NULL, doCalib = TRUE, doPred = FALSE, na.value = -9999, run = TRUE, preClean = FALSE, dir.exe = .BAM_PATH, name.exe = "BaM", predMaster_fname = "Config_Pred_Master.txt", stout = "" )BaM( workspace, mod, data, remnant = rep(list(remnantErrorModel()), mod$nY), mcmc = mcmcOptions(), cook = mcmcCooking(), summary = mcmcSummary(), residuals = residualOptions(), pred = NULL, doCalib = TRUE, doPred = FALSE, na.value = -9999, run = TRUE, preClean = FALSE, dir.exe = .BAM_PATH, name.exe = "BaM", predMaster_fname = "Config_Pred_Master.txt", stout = "" )
workspace |
Character, directory where config and result files are stored. |
mod |
model object, the model to be calibrated |
data |
dataset object, calibration data |
remnant |
list of remnantErrorModel objects. WARNING: make sure you use a list of length mod$nY (even if mod$nY=1!) |
mcmc |
mcmcOptions object, MCMC simulation number and options |
cook |
mcmcCooking object, properties of MCMC cooking (burn and slice) |
summary |
mcmcSummary object, properties of MCMC summary |
residuals |
residualOptions object, properties of residual analysis |
pred |
list of prediction objects, properties of prediction experiments |
doCalib |
Logical, do Calibration? (mcmc+cooking+summary+residuals) |
doPred |
Logical, do Prediction? |
na.value |
numeric, value used for NAs when writing the dataset in BaM format |
run |
Logical, run BaM? if FALSE, just write config files. |
preClean |
Logical, start by cleaning up workspace? Be careful, this will delete all files in the workspace, including old results! |
dir.exe |
Character, directory where BaM executable stands. |
name.exe |
Character, name of the executable without extension ('BaM' by default). |
predMaster_fname |
Character, name of configuration file pointing to all prediction experiments. |
stout |
Character string, standard output (see ?system2). In particular, stout="" (default) shows BaM messages in the console, stout=NULL discards BaM messages, stout='log.txt' saves BaM messages in file "log.txt". |
Nothing: just write config files and runs the executable.
# Fitting a rating curve - see https://github.com/BaM-tools/RBaM workspace=tempdir() D=dataset(X=SauzeGaugings['H'],Y=SauzeGaugings['Q'],Yu=SauzeGaugings['uQ'],data.dir=workspace) # Parameters of the low flow section control: activation stage k, coefficient a and exponent c k1=parameter(name='k1',init=-0.5,prior.dist='Uniform',prior.par=c(-1.5,0)) a1=parameter(name='a1',init=50,prior.dist='LogNormal',prior.par=c(log(50),1)) c1=parameter(name='c1',init=1.5,prior.dist='Gaussian',prior.par=c(1.5,0.05)) # Parameters of the high flow channel control: activation stage k, coefficient a and exponent c k2=parameter(name='k2',init=1,prior.dist='Gaussian',prior.par=c(1,1)) a2=parameter(name='a2',init=100,prior.dist='LogNormal',prior.par=c(log(100),1)) c2=parameter(name='c2',init=1.67,prior.dist='Gaussian',prior.par=c(1.67,0.05)) # Define control matrix: columns are controls, rows are stage ranges. controlMatrix=rbind(c(1,0),c(0,1)) # Stitch it all together into a model object M=model(ID='BaRatin', nX=1,nY=1, # number of input/output variables par=list(k1,a1,c1,k2,a2,c2), # list of model parameters xtra=xtraModelInfo(object=controlMatrix)) # use xtraModelInfo() to pass the control matrix # Call BaM to write configuration files. To actually run BaM, use run=TRUE, # but BaM executable needs to be downloaded first (use downloadBaM()) BaM(workspace=workspace,mod=M,data=D,run=FALSE)# Fitting a rating curve - see https://github.com/BaM-tools/RBaM workspace=tempdir() D=dataset(X=SauzeGaugings['H'],Y=SauzeGaugings['Q'],Yu=SauzeGaugings['uQ'],data.dir=workspace) # Parameters of the low flow section control: activation stage k, coefficient a and exponent c k1=parameter(name='k1',init=-0.5,prior.dist='Uniform',prior.par=c(-1.5,0)) a1=parameter(name='a1',init=50,prior.dist='LogNormal',prior.par=c(log(50),1)) c1=parameter(name='c1',init=1.5,prior.dist='Gaussian',prior.par=c(1.5,0.05)) # Parameters of the high flow channel control: activation stage k, coefficient a and exponent c k2=parameter(name='k2',init=1,prior.dist='Gaussian',prior.par=c(1,1)) a2=parameter(name='a2',init=100,prior.dist='LogNormal',prior.par=c(log(100),1)) c2=parameter(name='c2',init=1.67,prior.dist='Gaussian',prior.par=c(1.67,0.05)) # Define control matrix: columns are controls, rows are stage ranges. controlMatrix=rbind(c(1,0),c(0,1)) # Stitch it all together into a model object M=model(ID='BaRatin', nX=1,nY=1, # number of input/output variables par=list(k1,a1,c1,k2,a2,c2), # list of model parameters xtra=xtraModelInfo(object=controlMatrix)) # use xtraModelInfo() to pass the control matrix # Call BaM to write configuration files. To actually run BaM, use run=TRUE, # but BaM executable needs to be downloaded first (use downloadBaM()) BaM(workspace=workspace,mod=M,data=D,run=FALSE)
This function creates a square bloc-diagonal matrix from a list of square blocs.
blocDiag(blocs)blocDiag(blocs)
blocs |
list, each element is a square matrix. |
A square matrix.
blocDiag(list(1,cbind(c(1,2),c(3,4)),25))blocDiag(list(1,cbind(c(1,2),c(3,4)),25))
Creates a new instance of a 'dataset' object
dataset( X, Y, data.dir = getwd(), data.fname = "CalibrationData.txt", fname = "Config_Data.txt", Xu = NULL, Xb = NULL, Xb.indx = NULL, Yu = NULL, Yb = NULL, Yb.indx = NULL, VAR.indx = NULL )dataset( X, Y, data.dir = getwd(), data.fname = "CalibrationData.txt", fname = "Config_Data.txt", Xu = NULL, Xb = NULL, Xb.indx = NULL, Yu = NULL, Yb = NULL, Yb.indx = NULL, VAR.indx = NULL )
X |
data frame, observed input variables. |
Y |
data frame, observed output variables (same number of rows as X). |
data.dir |
Character, directory where a copy of the dataset will be written if required. Default is the current working directory, but you may prefer to use the BaM workspace. |
data.fname |
Character, data file name. |
fname |
Character, configuration file name. |
Xu |
data frame, random uncertainty in X, expressed as a standard deviation. Same dimension as X. |
Xb |
data frame, systematic uncertainty in X, expressed as a standard deviation. Same dimension as X. |
Xb.indx |
data frame, index of systematic errors in X. Same dimension as X. |
Yu |
data frame, random uncertainty in Y, expressed as a standard deviation. Same dimension as Y. |
Yb |
data frame, systematic uncertainty in Y, expressed as a standard deviation. Same dimension as Y. |
Yb.indx |
data frame, index of systematic errors in Y. Same dimension as Y. |
VAR.indx |
data frame, indices used for defining how VAR parameters vary. |
An object of class 'dataset'.
X=data.frame(input1=rnorm(100),input2=rnorm(100)) Y=data.frame(output=X$input1+0.8*X$input2+0.1*rnorm(100)) workspace=tempdir() d <- dataset(X=X,Y=Y,data.dir=workspace)X=data.frame(input1=rnorm(100),input2=rnorm(100)) Y=data.frame(output=X$input1+0.8*X$input2+0.1*rnorm(100)) workspace=tempdir() d <- dataset(X=X,Y=Y,data.dir=workspace)
returns a histogram+density ggplot (or a list thereof if several columns in sim)
densityPlot(sim, xlab = "values", col = "black")densityPlot(sim, xlab = "values", col = "black")
sim |
vector or matrix or data frame, MCMC simulations |
xlab |
Character, label of x-axis to be used if sim has no names |
col |
Color |
A ggplot (or a list thereof if several columns in sim)
# Create Monte Carlo samples n=1000 sim=data.frame(p1=rnorm(n),p2=rlnorm(n),p3=runif(n)) # create density plot for each component figures=densityPlot(sim)# Create Monte Carlo samples n=1000 sim=data.frame(p1=rnorm(n),p2=rlnorm(n),p3=runif(n)) # create density plot for each component figures=densityPlot(sim)
Download BaM executable
downloadBaM( destFolder, url = NULL, os = Sys.info()["sysname"], quiet = FALSE, ... )downloadBaM( destFolder, url = NULL, os = Sys.info()["sysname"], quiet = FALSE, ... )
destFolder |
character string, folder where BaM executable will be downloaded. |
url |
character string, the url from which BaM should be downloaded. When NULL, the url is determined automatically by using GitHub API to determine the latest release and the file corresponding to the OS. |
os |
character string, operating system, e.g. 'Linux', 'Windows' or 'Darwin'. |
quiet |
logical, if TRUE, suppress status messages. |
... |
arguments passed to function 'download.file' |
nothing - just download the file.
try(downloadBaM(destFolder=tempdir()))try(downloadBaM(destFolder=tempdir()))
Compute Area A(h), width w(h) and wet perimeter P(h) from a bathymetry profile (a,z).
getAwPfromBathy( bathy, hgrid = seq(min(bathy[, 2]), max(bathy[, 2]), diff(range(bathy[, 2]))/1000), segmentLength = sum(sqrt(apply(apply(bathy, 2, diff)^2, 1, sum)))/1000 )getAwPfromBathy( bathy, hgrid = seq(min(bathy[, 2]), max(bathy[, 2]), diff(range(bathy[, 2]))/1000), segmentLength = sum(sqrt(apply(apply(bathy, 2, diff)^2, 1, sum)))/1000 )
bathy |
data frame, 2 columns containing abscissa (increasing values) and stage. |
hgrid |
numeric vector, grid of h values where A, w and P are computed. By default 1000 values in the range of bathymetry's z. |
segmentLength |
numeric, segment length for bathymetry subsampling. By default 1/1000 of the total bathymetry's perimeter. |
A 4-column dataframe containing h, A(h), w(h) and P(h)
bathy=data.frame(a=c(0,0,0,1,2,2,4,6,8),h=c(3,2,0,-0.5,0,2,2.0001,2.3,3)) plot(bathy,type='l') df=getAwPfromBathy(bathy) plot(df$h,df$A,type='l') plot(df$h,df$w,type='l') plot(df$h,df$P,type='l')bathy=data.frame(a=c(0,0,0,1,2,2,4,6,8),h=c(3,2,0,-0.5,0,2,2.0001,2.3,3)) plot(bathy,type='l') df=getAwPfromBathy(bathy) plot(df$h,df$A,type='l') plot(df$h,df$w,type='l') plot(df$h,df$P,type='l')
Distributions and models available in BaM
getCatalogue(printOnly = FALSE)getCatalogue(printOnly = FALSE)
printOnly |
Logical, should the catalogue be returned or only printed? |
If printOnly==FALSE, a list with the following fields:
available univariate distributions.
available models.
catalogue <- getCatalogue() getCatalogue(printOnly=TRUE)catalogue <- getCatalogue() getCatalogue(printOnly=TRUE)
Get the initial values of a list of parameters and return them as a vector.
getInitPar(parameters)getInitPar(parameters)
parameters |
List of parameter objects, parameters whose initial values are sought. |
A numeric vector containing the initial values
ps <- list(parameter(name='par1',init=0),parameter(name='par2',init=1)) getInitPar(ps)ps <- list(parameter(name='par1',init=0),parameter(name='par2',init=1)) getInitPar(ps)
getNames from an object or a list of objects having a $name field (e.g. parameters)
getNames(loo, name = "name")getNames(loo, name = "name")
loo |
List Of Objects |
name |
character, string denoting the name field |
A character vector containing names
pars <- list(parameter(name='par1',init=0), parameter(name='par2',init=0), parameter(name='Third parameter',init=0)) getNames(pars)pars <- list(parameter(name='par1',init=0), parameter(name='par2',init=0), parameter(name='Third parameter',init=0)) getNames(pars)
Get parameter names for a distribution d
getParNames(d)getParNames(d)
d |
Character (possibly vector), distribution (possibly distributions) |
A character vector with parameter names.
parnames <- getParNames('GEV') npar <- length(getParNames('Gumbel'))parnames <- getParNames('GEV') npar <- length(getParNames('Gumbel'))
Get path to BaM executable
getPathToBaM()getPathToBaM()
A string containing the path of the directory containing BaM executable, or NULL if this directory has not been set yet.
getPathToBaM()getPathToBaM()
Computes the log-likelihood from model-simulated values, based on a Gaussian iid error model:
Yobs = Ysim + delta + epsilon
Measurement errors: delta~ N(0,sdev=Yu)
Structural errors: epsilon~ N(0,sdev=gamma)
If Yobs/Ysim are multi-variate, this error model is applied independently to each component.
llfunk_iid_Gaussian(Ysim, Yobs, Yu, gamma)llfunk_iid_Gaussian(Ysim, Yobs, Yu, gamma)
Ysim |
data frame, model-simulated values. |
Yobs |
data frame, corresponding observed values, same dimensions as Ysim. NAs are skipped. |
Yu |
data frame, measurement uncertainties (standard deviations), same dimensions as Ysim and Yobs. |
gamma |
numeric vector, structural error parameters. length(gamma) = number of columns in Ysim. |
A numeric value equal to the log-likelihood.
Yobs=SauzeGaugings['Q'] Yu=SauzeGaugings['uQ'] Ysim=100*(SauzeGaugings['H']+0.5)^1.6 llfunk_iid_Gaussian(Ysim,Yobs,Yu,gamma=100)Yobs=SauzeGaugings['Q'] Yu=SauzeGaugings['uQ'] Ysim=100*(SauzeGaugings['H']+0.5)^1.6 llfunk_iid_Gaussian(Ysim,Yobs,Yu,gamma=100)
Computes the log-likelihood from model-simulated values based on a Gaussian independent error model with linearly-varying standard deviation:
Yobs = Ysim + delta + epsilon
Measurement errors: delta~ N(0,sdev=Yu)
Structural errors: epsilon~ N(0,sdev=g1+g2*|Ysim|)
If Yobs/Ysim are multi-variate, this error model is applied independently to each component.
llfunk_iLinear_Gaussian(Ysim, Yobs, Yu, gamma)llfunk_iLinear_Gaussian(Ysim, Yobs, Yu, gamma)
Ysim |
data frame, model-simulated values. |
Yobs |
data frame, corresponding observed values, same dimensions as Ysim. NAs are skipped. |
Yu |
data frame, measurement uncertainties (standard deviations), same dimensions as Ysim and Yobs. |
gamma |
numeric vector, structural error parameters, organized as: gamma=c((g1,g2) for the 1st component of Ysim,(g1,g2) for the 2nd component of Ysim, etc.) => length(gamma) = 2*(number of columns in Ysim). |
A numeric value equal to the log-likelihood.
Yobs=SauzeGaugings['Q'] Yu=SauzeGaugings['uQ'] Ysim=100*(SauzeGaugings['H']+0.5)^1.6 llfunk_iLinear_Gaussian(Ysim,Yobs,Yu,gamma=c(1,0.1))Yobs=SauzeGaugings['Q'] Yu=SauzeGaugings['uQ'] Ysim=100*(SauzeGaugings['H']+0.5)^1.6 llfunk_iLinear_Gaussian(Ysim,Yobs,Yu,gamma=c(1,0.1))
Log-likelihood engine for a model available in BaM. Unlike functions llfunk_***, which compute the log-likelihood from model-simulated values Ysim (see e.g. llfunk_iid_Gaussian), this function computes the log-likelihood from (parameters + inputs X + model mod).
logLikelihood_BaM( parvector, X, Yobs, Yu, llfunk, mod, Ysim = NULL, llargs = NULL )logLikelihood_BaM( parvector, X, Yobs, Yu, llfunk, mod, Ysim = NULL, llargs = NULL )
parvector |
numeric vector, parameter vector, including thetas (model parameters) and gammas (structural errors parameters). |
X |
data frame, model inputs. |
Yobs |
data frame, corresponding observed values. |
Yu |
data frame, measurement uncertainties (standard deviations), same dimensions as Yobs. |
llfunk |
function, function computing the log-likelihood given Ysim, see e.g. llfunk_iid_Gaussian. |
mod |
model object, the model to be calibrated. |
Ysim |
data frame, model-simulated values. When NULL (default), the model is run to provide simulations. When a non-NULL data frame is provided, it is used as pre-computed simulations, and the model is hence not run. This is useful to speed-up some MCMC strategies. |
llargs |
object, any other arguments to be passed to llfunk. |
A list with the following components:
logLikelihood |
numeric, the log-likelihood. |
Ysim |
data frame, the model-simulated values. |
# Single-control rating curve model - see https://github.com/BaM-tools/RBaM # Parameters are activation stage k, coefficient a and exponent c k=parameter(name='k',init=-0.5) a=parameter(name='a',init=100) c=parameter(name='c',init=1.6) # Define control matrix: columns are controls, rows are stage ranges. controlMatrix=matrix(1,nrow=1,ncol=1) # Stitch it all together into a model object M=model(ID='BaRatin', nX=1,nY=1, # number of input/output variables par=list(k,a,c), # list of model parameters xtra=xtraModelInfo(object=controlMatrix)) # use xtraModelInfo() to pass the control matrix # Define calibration data X=SauzeGaugings['H'] Yobs=SauzeGaugings['Q'] Yu=SauzeGaugings['uQ'] # Define the parameter vector (model parameters + structural error parameters) parvector=c(RBaM::getInitPar(M$par),c(1,0.1)) # Compute log-likelihood logLikelihood_BaM(parvector=parvector,X=X,Yobs=Yobs,Yu=Yu, llfunk=llfunk_iLinear_Gaussian,mod=M)# Single-control rating curve model - see https://github.com/BaM-tools/RBaM # Parameters are activation stage k, coefficient a and exponent c k=parameter(name='k',init=-0.5) a=parameter(name='a',init=100) c=parameter(name='c',init=1.6) # Define control matrix: columns are controls, rows are stage ranges. controlMatrix=matrix(1,nrow=1,ncol=1) # Stitch it all together into a model object M=model(ID='BaRatin', nX=1,nY=1, # number of input/output variables par=list(k,a,c), # list of model parameters xtra=xtraModelInfo(object=controlMatrix)) # use xtraModelInfo() to pass the control matrix # Define calibration data X=SauzeGaugings['H'] Yobs=SauzeGaugings['Q'] Yu=SauzeGaugings['uQ'] # Define the parameter vector (model parameters + structural error parameters) parvector=c(RBaM::getInitPar(M$par),c(1,0.1)) # Compute log-likelihood logLikelihood_BaM(parvector=parvector,X=X,Yobs=Yobs,Yu=Yu, llfunk=llfunk_iLinear_Gaussian,mod=M)
Log-posterior engine for a model available in BaM.
logPosterior_BaM( parvector, X, Yobs, Yu, lpfunk, llfunk, mod, Ysim = NULL, logLikelihood_engine = logLikelihood_BaM, llargs = NULL )logPosterior_BaM( parvector, X, Yobs, Yu, lpfunk, llfunk, mod, Ysim = NULL, logLikelihood_engine = logLikelihood_BaM, llargs = NULL )
parvector |
numeric vector, parameter vector, including thetas (model parameters) and gammas (structural errors parameters). |
X |
data frame, model inputs. |
Yobs |
data frame, corresponding observed values. |
Yu |
data frame, measurement uncertainties (standard deviations), same dimensions as Yobs. |
lpfunk |
function, function computing the log-prior density of parvector. |
llfunk |
function, function computing the log-likelihood given Ysim, see e.g. llfunk_iid_Gaussian. |
mod |
model object, the model to be calibrated. |
Ysim |
data frame, model-simulated values. When NULL (default), the model is run to provide simulations. When a non-NULL data frame is provided, it is used as pre-computed simulations, and the model is hence not run. This is useful to speed-up some MCMC strategies. |
logLikelihood_engine |
function, engine function used to compute the log-likelihood, see e.g. logLikelihood_BaM. Unlike functions llfunk, which computes the log-likelihood from model-simulated values Ysim (see e.g. llfunk_iid_Gaussian, logLikelihood_engine computes the log-likelihood from (parameters + inputs X + model mod). |
llargs |
object, any other arguments to be passed to llfunk. |
A list with the following components:
logPosterior |
numeric, the unnormalized posterior log-pdf. |
logPrior |
numeric, the prior log-pdf. |
logLikelihood |
numeric, the log-likelihood. |
Ysim |
data frame, the model-simulated values. |
# Single-control rating curve model - see https://github.com/BaM-tools/RBaM # Parameters are activation stage k, coefficient a and exponent c k=parameter(name='k',init=-0.5) a=parameter(name='a',init=100) c=parameter(name='c',init=1.6) # Define control matrix: columns are controls, rows are stage ranges. controlMatrix=matrix(1,nrow=1,ncol=1) # Stitch it all together into a model object M=model(ID='BaRatin', nX=1,nY=1, # number of input/output variables par=list(k,a,c), # list of model parameters xtra=xtraModelInfo(object=controlMatrix)) # use xtraModelInfo() to pass the control matrix # Define calibration data X=SauzeGaugings['H'] Yobs=SauzeGaugings['Q'] Yu=SauzeGaugings['uQ'] # Define the parameter vector (model parameters + structural error parameters) parvector=c(RBaM::getInitPar(M$par),c(1,0.1)) # Define prior function - here for instance just an informative prior on the third parameter myPrior <-function(parvector){dnorm(parvector[3],1.6,0.1,log=TRUE)} # Compute log-likelihood logPosterior_BaM(parvector=parvector,X=X,Yobs=Yobs,Yu=Yu, lpfunk=myPrior,llfunk=llfunk_iLinear_Gaussian,mod=M)# Single-control rating curve model - see https://github.com/BaM-tools/RBaM # Parameters are activation stage k, coefficient a and exponent c k=parameter(name='k',init=-0.5) a=parameter(name='a',init=100) c=parameter(name='c',init=1.6) # Define control matrix: columns are controls, rows are stage ranges. controlMatrix=matrix(1,nrow=1,ncol=1) # Stitch it all together into a model object M=model(ID='BaRatin', nX=1,nY=1, # number of input/output variables par=list(k,a,c), # list of model parameters xtra=xtraModelInfo(object=controlMatrix)) # use xtraModelInfo() to pass the control matrix # Define calibration data X=SauzeGaugings['H'] Yobs=SauzeGaugings['Q'] Yu=SauzeGaugings['uQ'] # Define the parameter vector (model parameters + structural error parameters) parvector=c(RBaM::getInitPar(M$par),c(1,0.1)) # Define prior function - here for instance just an informative prior on the third parameter myPrior <-function(parvector){dnorm(parvector[3],1.6,0.1,log=TRUE)} # Compute log-likelihood logPosterior_BaM(parvector=parvector,X=X,Yobs=Yobs,Yu=Yu, lpfunk=myPrior,llfunk=llfunk_iLinear_Gaussian,mod=M)
Log-posterior engine for a model available in BaM. This is a wrapped version of logPosterior_BaM, returning a single numeric value (the log-posterior). This function can hence be passed to standard optimization (e.g. optim) or MCMC (e.g. metrop) tools.
logPosterior_BaM_wrapped( parvector, X, Yobs, Yu, lpfunk, llfunk, mod, Ysim = NULL, logLikelihood_engine = logLikelihood_BaM, llargs = NULL )logPosterior_BaM_wrapped( parvector, X, Yobs, Yu, lpfunk, llfunk, mod, Ysim = NULL, logLikelihood_engine = logLikelihood_BaM, llargs = NULL )
parvector |
numeric vector, parameter vector, including thetas (model parameters) and gammas (structural errors parameters). |
X |
data frame, model inputs. |
Yobs |
data frame, corresponding observed values. |
Yu |
data frame, measurement uncertainties (standard deviations), same dimensions as Yobs. |
lpfunk |
function, function computing the log-prior density of parvector. |
llfunk |
function, function computing the log-likelihood given Ysim, see e.g. llfunk_iid_Gaussian. |
mod |
model object, the model to be calibrated. |
Ysim |
data frame, model-simulated values. When NULL (default), the model is run to provide simulations. When a non-NULL data frame is provided, it is used as pre-computed simulations, and the model is hence not run. This is useful to speed-up some MCMC strategies. |
logLikelihood_engine |
function, engine function used to compute the log-likelihood, see e.g. logLikelihood_BaM. Unlike functions llfunk, which computes the log-likelihood from model-simulated values Ysim (see e.g. llfunk_iid_Gaussian, logLikelihood_engine computes the log-likelihood from (parameters + inputs X + model mod). |
llargs |
object, any other arguments to be passed to llfunk. |
the unnormalized posterior log-pdf (numeric).
# Single-control rating curve model - see https://github.com/BaM-tools/RBaM # Parameters are activation stage k, coefficient a and exponent c k=parameter(name='k',init=-0.5) a=parameter(name='a',init=100) c=parameter(name='c',init=1.6) # Define control matrix: columns are controls, rows are stage ranges. controlMatrix=matrix(1,nrow=1,ncol=1) # Stitch it all together into a model object M=model(ID='BaRatin', nX=1,nY=1, # number of input/output variables par=list(k,a,c), # list of model parameters xtra=xtraModelInfo(object=controlMatrix)) # use xtraModelInfo() to pass the control matrix # Define calibration data X=SauzeGaugings['H'] Yobs=SauzeGaugings['Q'] Yu=SauzeGaugings['uQ'] # Define the parameter vector (model parameters + structural error parameters) parvector=c(RBaM::getInitPar(M$par),c(1,0.1)) # Define prior function - here for instance just an informative prior on the third parameter myPrior <-function(parvector){dnorm(parvector[3],1.6,0.1,log=TRUE)} # Compute log-likelihood logPosterior_BaM_wrapped(parvector=parvector,X=X,Yobs=Yobs,Yu=Yu, lpfunk=myPrior,llfunk=llfunk_iLinear_Gaussian,mod=M)# Single-control rating curve model - see https://github.com/BaM-tools/RBaM # Parameters are activation stage k, coefficient a and exponent c k=parameter(name='k',init=-0.5) a=parameter(name='a',init=100) c=parameter(name='c',init=1.6) # Define control matrix: columns are controls, rows are stage ranges. controlMatrix=matrix(1,nrow=1,ncol=1) # Stitch it all together into a model object M=model(ID='BaRatin', nX=1,nY=1, # number of input/output variables par=list(k,a,c), # list of model parameters xtra=xtraModelInfo(object=controlMatrix)) # use xtraModelInfo() to pass the control matrix # Define calibration data X=SauzeGaugings['H'] Yobs=SauzeGaugings['Q'] Yu=SauzeGaugings['uQ'] # Define the parameter vector (model parameters + structural error parameters) parvector=c(RBaM::getInitPar(M$par),c(1,0.1)) # Define prior function - here for instance just an informative prior on the third parameter myPrior <-function(parvector){dnorm(parvector[3],1.6,0.1,log=TRUE)} # Compute log-likelihood logPosterior_BaM_wrapped(parvector=parvector,X=X,Yobs=Yobs,Yu=Yu, lpfunk=myPrior,llfunk=llfunk_iLinear_Gaussian,mod=M)
Computes the log-density of an improper flat prior distribution.
logPrior_Flat(parvector)logPrior_Flat(parvector)
parvector |
numeric vector, parameter vector, including thetas (model parameters) and gammas (structural errors parameters). |
A numeric value equal to the prior log-density.
logPrior_Flat(c(1,1,0.2))logPrior_Flat(c(1,1,0.2))
An adaptive Metropolis sampler largely inspired by Haario et al. (2001, https://www.jstor.org/stable/3318737). The jump covariance is adapted using the empirical covariance of previously-sampled values, and the scaling factor is adapted in order to comply with a specified move rate interval.
MCMC_AM( logPdf, x0, C0 = diag((0.01 * (abs(x0) + 0.1))^2), scaleFactor = 2.4/sqrt(length(x0)), nAdapt = 50, nCycles = 20, minMoveRate = 0.2, maxMoveRate = 0.5, downMult = 0.9, upMult = 1.1, burnCov = 0.2, dofCovMin = 10, nCovMax = 1000, ... )MCMC_AM( logPdf, x0, C0 = diag((0.01 * (abs(x0) + 0.1))^2), scaleFactor = 2.4/sqrt(length(x0)), nAdapt = 50, nCycles = 20, minMoveRate = 0.2, maxMoveRate = 0.5, downMult = 0.9, upMult = 1.1, burnCov = 0.2, dofCovMin = 10, nCovMax = 1000, ... )
logPdf |
function, evaluating the log-density of the distribution to sample from (up to a proportionality constant). logPdf can return either a single numeric value, interpreted as the target log-pdf, or a list containing components named 'logPosterior', 'logLikelihood' and 'logPrior'. |
x0 |
numeric vector, starting point |
C0 |
numeric matrix, covariance matrix of the Gaussian jump distribution (up to a scale factor, see next). |
scaleFactor |
numeric >0, used to scale the jump covariance. The covariance of the jump distribution is equal to (scaleFactor^2)*C0 |
nAdapt |
integer > 1, number of iterations before adapting covariance C and scaleFactor. |
nCycles |
integer > 1, number of adaption cycles. Total number of iterations is hence equal to nAdapt*nCycles. nCycles=1 leads to the standard non-adaptive Metropolis sampler. |
minMoveRate |
numeric in (0;1), lower bound for the desired move rate interval. |
maxMoveRate |
numeric in (0;1), upper bound for the desired move rate interval. |
downMult |
numeric in (0;1), multiplicative factor used to decrease scaleFactor when move rate is too low. |
upMult |
numeric (>1, avoid 1/downMult) multiplicative factor used to increase scaleFactor when move rate is too high. |
burnCov |
numeric in (0;1), fraction of initial values to be discarded before computing the empirical covariance of sampled vectors, which is used to adapt the jump covariance. |
dofCovMin |
integer, minimum number of degrees of freedom required to compute the empirical covariance of sampled vectors and hence to adapt the jump covariance. If D denotes the length of x0, at least dofCovMin*(D+0.5*(D-1)*(D-2)) iterations are required before adapting the jump covariance (i.e. dofCovMin times the number of unknown elements in the covariance matrix). |
nCovMax |
integer, maximum number of iterations used to compute the empirical covariance. If the number of available iterations is larger than nCovMax, iterations are 'slimmed' to reach nCovMax. |
... |
other arguments passed to function logPdf |
A list with the following components:
samples |
data frame, MCMC simulations. |
components |
data frame, corresponding values of the log-posterior, the log-prior and the log-likelihood. |
C |
matrix, the adapted jump covariance matrix. |
scaleFactor |
numeric, the adapted scaling factor. |
# Define a 2-dimensional target log-pdf logPdf <- function(x){ p1=log(0.6*dnorm(x[1],0,1)+0.4*dnorm(x[1],2,0.5)) p2=log(dlnorm(x[2],0,1)) return(p1+p2) } # Sample from it mcmc=MCMC_AM(logPdf,c(1,1)) plot(mcmc$samples)# Define a 2-dimensional target log-pdf logPdf <- function(x){ p1=log(0.6*dnorm(x[1],0,1)+0.4*dnorm(x[1],2,0.5)) p2=log(dlnorm(x[2],0,1)) return(p1+p2) } # Sample from it mcmc=MCMC_AM(logPdf,c(1,1)) plot(mcmc$samples)
An adaptive Metropolis sampler that updates the parameter vector one component at a time using a 1-dimensional jump. This allows easily adapting the jump standard deviation for each component in order to comply with a specified move rate interval.
MCMC_OAAT( logPdf, x0, s0 = 0.05 * (abs(x0) + 0.1), nTheta = length(x0), nAdapt = 50, nCycles = 20, minMoveRate = 0.2, maxMoveRate = 0.5, downMult = 0.9, upMult = 1.1, ... )MCMC_OAAT( logPdf, x0, s0 = 0.05 * (abs(x0) + 0.1), nTheta = length(x0), nAdapt = 50, nCycles = 20, minMoveRate = 0.2, maxMoveRate = 0.5, downMult = 0.9, upMult = 1.1, ... )
logPdf |
function, evaluating the log-density of the distribution to sample from (up to a proportionality constant). logPdf can return either a single numeric value, interpreted as the target log-pdf, or a list containing components named 'logPosterior', 'logLikelihood' and 'logPrior'. |
x0 |
numeric vector, starting point. |
s0 |
numeric vector, starting jump standard deviations. |
nTheta |
integer>0, size of the "theta" part of x0, i.e. components that represent the model parameters rather than structural errors parameters (gamma). This is used to speed-up the sampler by avoiding running the model for gamma components. nTheta=length(x0) (default) implies no attempt at speeding up. |
nAdapt |
integer > 1, number of iterations before adapting the jump standard deviations. |
nCycles |
integer > 1, number of adaption cycles. Total number of iterations is hence equal to nAdapt*nCycles. nCycles=1 leads to a non-adaptive one-at-a-time Metropolis sampler. |
minMoveRate |
numeric in (0;1), lower bound for the desired move rate interval. |
maxMoveRate |
numeric in (0;1), upper bound for the desired move rate interval. |
downMult |
numeric in (0;1), multiplication factor used to decrease the jump standard deviation when move rate is too low. |
upMult |
numeric (>1, avoid 1/downMult) multiplication factor used to increase the jump standard deviations when move rate is too high. |
... |
other arguments passed to function logPdf |
A list with the following components:
samples |
data frame, MCMC simulations. |
components |
data frame, corresponding values of the log-posterior, the log-prior and the log-likelihood. |
sjump |
numeric vector, the adapted jump standard deviations. |
# Define a 2-dimensional target log-pdf logPdf <- function(x){ p1=log(0.6*dnorm(x[1],0,1)+0.4*dnorm(x[1],2,0.5)) p2=log(dlnorm(x[2],0,1)) return(p1+p2) } # Sample from it mcmc=MCMC_OAAT(logPdf,c(1,1)) plot(mcmc$samples)# Define a 2-dimensional target log-pdf logPdf <- function(x){ p1=log(0.6*dnorm(x[1],0,1)+0.4*dnorm(x[1],2,0.5)) p2=log(dlnorm(x[2],0,1)) return(p1+p2) } # Sample from it mcmc=MCMC_OAAT(logPdf,c(1,1)) plot(mcmc$samples)
Creates a new instance of a 'mcmcCooking' object
mcmcCooking( fname = "Config_Cooking.txt", result.fname = "Results_Cooking.txt", burn = 0.5, nSlim = 10 )mcmcCooking( fname = "Config_Cooking.txt", result.fname = "Results_Cooking.txt", burn = 0.5, nSlim = 10 )
fname |
Character, configuration file name. |
result.fname |
Character, result file name. |
burn |
numeric, burn factor, >=0 and <1. 0.4 means the first 40 percent of MCMC samples are discarded). |
nSlim |
Integer, slimming period: 10 means only one MCMC sample every 10 is kept (after burning). |
An object of class 'mcmcCooking'.
m <- mcmcCooking()m <- mcmcCooking()
Creates a new instance of a 'mcmcOptions' object
mcmcOptions( fname = "Config_MCMC.txt", result.fname = "Results_MCMC.txt", nAdapt = 100, nCycles = 100, minMoveRate = 0.1, maxMoveRate = 0.5, downMult = 0.9, upMult = 1.1, multFactor = 0.1, manualMode = FALSE, thetaStd = 9999, gammaStd = list(9999) )mcmcOptions( fname = "Config_MCMC.txt", result.fname = "Results_MCMC.txt", nAdapt = 100, nCycles = 100, minMoveRate = 0.1, maxMoveRate = 0.5, downMult = 0.9, upMult = 1.1, multFactor = 0.1, manualMode = FALSE, thetaStd = 9999, gammaStd = list(9999) )
fname |
Character, configuration file name. |
result.fname |
Character, result file name. |
nAdapt |
Integer, adaptation period: jump sizes are increased/decreased every Nadapt iterations to comply with the desired moving rates. |
nCycles |
Integer, number of adaptation cycles (total number of iterations is hence Nadapt * Ncycles). |
minMoveRate |
Numeric in (0;1), lower bound for the desired move rate interval. |
maxMoveRate |
Numeric in (0;1), upper bound for the desired move rate interval. |
downMult |
Numeric in (0:1), multiplication factor used to decrease jump size when move rate is too low. |
upMult |
Numeric (>1, avoid 1/dowMult) multiplication factor used to increase jump size when move rate is too high. |
multFactor |
Numeric >0, multiplicative factor to set initial jump standard deviations to multFactor*|initValue| (AUTO mode). |
manualMode |
logical, should jump standard deviations be entered manually? |
thetaStd |
Numeric vector (>0), jump standard deviations for model parameters theta (MANUAL mode). |
gammaStd |
list of numeric vectors (>0), size = number of output variables of the model. Jump standard deviations for structural error parameters gamma of each output variable (MANUAL mode). |
An object of class 'mcmcOptions'.
m <- mcmcOptions()m <- mcmcOptions()
Creates a new instance of a 'mcmcSummary' object
mcmcSummary( fname = "Config_Summary.txt", result.fname = "Results_Summary.txt", DIC.fname = "Results_DIC.txt", xtendedMCMC.fname = "" )mcmcSummary( fname = "Config_Summary.txt", result.fname = "Results_Summary.txt", DIC.fname = "Results_DIC.txt", xtendedMCMC.fname = "" )
fname |
Character, configuration file name. |
result.fname |
Character, summary file name. |
DIC.fname |
Character, DIC file name. Not computed if empty string. |
xtendedMCMC.fname |
Character, xtended MCMC file name. Not written if empty string. |
An object of class 'mcmcSummary'.
m <- mcmcSummary()m <- mcmcSummary()
Stage-discharge gaugings from the hydrometric station 'the Ardèche River at Meyras'. See https://en.wikipedia.org/wiki/Ardèche_(river) for a description of the river See https://doi.org/10.1029/2018WR023389 for an article using this dataset
MeyrasGaugingsMeyrasGaugings
A data frame with 104 rows and 4 variables:
Stage (m)
Discharge (m3/s)
Discharge uncertainty (m3/s) expressed as a standard deviation
Stability period on which a single rating curve can be used
Creates a new instance of a 'model' object
model( fname = "Config_Model.txt", ID = "Linear", nX = 1, nY = 1, par = list(parameter("Xeffect", 1, prior.dist = "FlatPrior")), xtra = xtraModelInfo() )model( fname = "Config_Model.txt", ID = "Linear", nX = 1, nY = 1, par = list(parameter("Xeffect", 1, prior.dist = "FlatPrior")), xtra = xtraModelInfo() )
fname |
Character, configuration file name. |
ID |
Character, model ID. Type 'getCatalogue()' for available models. |
nX |
Integer, number of input variables. |
nY |
Integer, number of output variables. |
par |
list of parameter objects, parameters of the model. |
xtra |
xtraModelInfo object. |
An object of class 'model'.
# defaut linear regression model Y=aX+b mod <- model() # BaRatin model for a single-control rating curve Y=a(X-b)^c mod <- model(ID='BaRatin',nX=1,nY=1, par=list(parameter('a',10,prior.dist='LogNormal',prior.par=c(log(10),0.1)), parameter('b',-1,prior.dist='Gaussian',prior.par=c(-1,1)), parameter('c',5/3,prior.dist='Gaussian',prior.par=c(5/3,0.05))), xtra=xtraModelInfo(object=matrix(1,nrow=1,ncol=1)))# defaut linear regression model Y=aX+b mod <- model() # BaRatin model for a single-control rating curve Y=a(X-b)^c mod <- model(ID='BaRatin',nX=1,nY=1, par=list(parameter('a',10,prior.dist='LogNormal',prior.par=c(log(10),0.1)), parameter('b',-1,prior.dist='Gaussian',prior.par=c(-1,1)), parameter('c',5/3,prior.dist='Gaussian',prior.par=c(5/3,0.05))), xtra=xtraModelInfo(object=matrix(1,nrow=1,ncol=1)))
Creates a new instance of a 'parameter' object
parameter(name, init, prior.dist = "FlatPrior", prior.par = NULL)parameter(name, init, prior.dist = "FlatPrior", prior.par = NULL)
name |
character, parameter name. |
init |
numeric, initial guess. |
prior.dist |
character, prior distribution. |
prior.par |
numeric vector, prior parameters |
An object of class 'parameter'.
p <- parameter(name='par',init=0,prior.dist='Gaussian',prior.par=c(0,1))p <- parameter(name='par',init=0,prior.dist='Gaussian',prior.par=c(0,1))
Creates a new instance of a 'parameter_VAR' object
parameter_VAR( name, index, d, init, prior.dist = rep("FlatPrior", length(init)), prior.par = rep(list(NULL), length(init)) )parameter_VAR( name, index, d, init, prior.dist = rep("FlatPrior", length(init)), prior.par = rep(list(NULL), length(init)) )
name |
character, parameter name. |
index |
character, name of column in VAR.indx (see ?dataset) containing the index for this varying parameter |
d |
dataset object, the dataset containing (amongst other things) the index above |
init |
numeric vector, initial guesses for each instance of the VAR parameter. |
prior.dist |
character vector, prior distribution for each instance of the VAR parameter. |
prior.par |
list of numeric vectors, prior parameters for each instance of the VAR parameter |
An object of class 'parameter_VAR'.
X=data.frame(input1=rnorm(100),input2=rnorm(100)) Y=data.frame(output=X$input1+0.8*X$input2+0.1*rnorm(100)) VAR.indx=data.frame(indx=c(rep(1,50),rep(2,50))) workspace=tempdir() d <- dataset(X=X,Y=Y,data.dir=workspace,VAR.indx=VAR.indx) p <- parameter_VAR(name='par',index='indx',d=d, init=c(-1,1,2), prior.dist=c('Gaussian','FlatPrior','Triangle'), prior.par=list(c(-1,1),NULL,c(2,0,5)))X=data.frame(input1=rnorm(100),input2=rnorm(100)) Y=data.frame(output=X$input1+0.8*X$input2+0.1*rnorm(100)) VAR.indx=data.frame(indx=c(rep(1,50),rep(2,50))) workspace=tempdir() d <- dataset(X=X,Y=Y,data.dir=workspace,VAR.indx=VAR.indx) p <- parameter_VAR(name='par',index='indx',d=d, init=c(-1,1,2), prior.dist=c('Gaussian','FlatPrior','Triangle'), prior.par=list(c(-1,1),NULL,c(2,0,5)))
Creates a new instance of a 'prediction' object
prediction( X, spagFiles, data.dir = getwd(), data.fnames = paste0("X", 1:length(X), ".pred"), fname = paste0("Config_Pred_", paste0(sample(c(letters, LETTERS, 0:9), 6), collapse = ""), ".txt"), doParametric = FALSE, doStructural = rep(FALSE, length(spagFiles)), transposeSpag = TRUE, priorNsim = NULL, envFiles = paste0(tools::file_path_sans_ext(spagFiles), ".env"), consoleProgress = TRUE, spagFiles_state = NULL, transposeSpag_state = TRUE, envFiles_state = switch(is.null(spagFiles_state) + 1, paste0(tools::file_path_sans_ext(spagFiles_state), ".env"), NULL), parSamples = NULL )prediction( X, spagFiles, data.dir = getwd(), data.fnames = paste0("X", 1:length(X), ".pred"), fname = paste0("Config_Pred_", paste0(sample(c(letters, LETTERS, 0:9), 6), collapse = ""), ".txt"), doParametric = FALSE, doStructural = rep(FALSE, length(spagFiles)), transposeSpag = TRUE, priorNsim = NULL, envFiles = paste0(tools::file_path_sans_ext(spagFiles), ".env"), consoleProgress = TRUE, spagFiles_state = NULL, transposeSpag_state = TRUE, envFiles_state = switch(is.null(spagFiles_state) + 1, paste0(tools::file_path_sans_ext(spagFiles_state), ".env"), NULL), parSamples = NULL )
X |
data frame or list of dataframes / matrices, representing the values taken by the input variables.
|
spagFiles |
Character vector (size nY, the number of output variables). Name of the files containing the spaghettis for each output variable. NOTE: provide file names only, not full paths. Using a '.spag' extension is a good practice. |
data.dir |
Character, directory where a copies of the dataset X will be written if required (1 file per variable in X). Default is the current working directory, but you may prefer to use the BaM workspace. |
data.fnames |
Character, data file names. |
fname |
Character, configuration file name. |
doParametric |
Logical, propagate parametric uncertainty? If FALSE, maxpost parameters are used. |
doStructural |
Logical, propagate structural uncertainty for each output variable? (size nY) |
transposeSpag |
Logical. If FALSE, spaghettis are written horizontally (row-wise), otherwise they will be transposed so that each spaghetti is a column. |
priorNsim |
Integer, number of samples from the prior distribution for 'prior prediction' experiments. If negative or NULL (default), posterior samples are used. |
envFiles |
Character vector (size nY, the number of output variables). Name of the files containing the envelops (e.g. prediction intervals) computed from the spaghettis for each output variable. By default, same name as spaghetti files but with a '.env' extension. If NULL, envelops are not computed. |
consoleProgress |
Logical, print progress in BaM.exe console? |
spagFiles_state |
Character vector (size nState, the number of state variables), same as spagFiles but for states rather than outputs. If NULL, states are not predicted. Note that only parametric uncertainty is propagated for state variables since they are not observed. Consequently, structural uncertainty = 0 and total uncertainty = parametric uncertainty. |
transposeSpag_state |
Logical. Same as transposeSpag, but for states rather than outputs. |
envFiles_state |
Character vector (size nState, the number of state variables), same as envFiles, but for states rather than outputs. |
parSamples |
data frame, parameter samples that will replace the MCMC-generated one for this prediction. |
An object of class 'prediction'.
#-------------- # Example using the twoPopulations dataset, containing 101 values for # 3 input variables (time t, temperature at site 1 T1, temperature at site 2 T2) # and 2 output variables (population at site 1 P1, population at site 2 P2). pred=prediction(X=twoPopulations[,1:3],spagFiles=c('P1.spag','P2.spag')) #-------------- # Alternative example showing how to propagate uncertainty in some of # the input variables (here, temperatures T1 and T2) # Create 100 noisy replicates for T1, representing uncertainty T1rep=matrix(rnorm(101*100,mean=twoPopulations$T1,sd=0.1),nrow=101,ncol=100) # Same for T2 T2rep=matrix(rnorm(101*100,mean=twoPopulations$T2,sd=0.1),nrow=101,ncol=100) # Create prediction object pred=prediction(X=list(twoPopulations$t,T1rep,T2rep),spagFiles=c('P1.spag','P2.spag'))#-------------- # Example using the twoPopulations dataset, containing 101 values for # 3 input variables (time t, temperature at site 1 T1, temperature at site 2 T2) # and 2 output variables (population at site 1 P1, population at site 2 P2). pred=prediction(X=twoPopulations[,1:3],spagFiles=c('P1.spag','P2.spag')) #-------------- # Alternative example showing how to propagate uncertainty in some of # the input variables (here, temperatures T1 and T2) # Create 100 noisy replicates for T1, representing uncertainty T1rep=matrix(rnorm(101*100,mean=twoPopulations$T1,sd=0.1),nrow=101,ncol=100) # Same for T2 T2rep=matrix(rnorm(101*100,mean=twoPopulations$T2,sd=0.1),nrow=101,ncol=100) # Create prediction object pred=prediction(X=list(twoPopulations$t,T1rep,T2rep),spagFiles=c('P1.spag','P2.spag'))
Read raw MCMC samples, return cooked (burnt & sliced) ones
readMCMC( file = "Results_Cooking.txt", burnFactor = 0, slimFactor = 1, sep = "", reportFile = NULL, panelPerCol = 10, panelHeight = 3, panelWidth = 23/panelPerCol )readMCMC( file = "Results_Cooking.txt", burnFactor = 0, slimFactor = 1, sep = "", reportFile = NULL, panelPerCol = 10, panelHeight = 3, panelWidth = 23/panelPerCol )
file |
Character, full path to MCMC file. |
burnFactor |
Numeric, burn factor. 0.1 means the first 10 are discarded. |
slimFactor |
Integer, slim factor. 10 means that only one iteration every 10 is kept. |
sep |
Character, separator used in MCMC file. |
reportFile |
Character, full path to pdf report file, not created if NULL |
panelPerCol |
Integer, max number of panels per column |
panelHeight |
Numeric, height of each panel |
panelWidth |
Numeric, width of each panel |
A data frame containing the cooked mcmc samples.
# Create Monte Carlo samples and write them to file n=4000 sim=data.frame(p1=rnorm(n),p2=rlnorm(n),p3=runif(n)) workspace=tempdir() write.table(sim,file=file.path(workspace,'MCMC.txt'),row.names=FALSE) # Read file, burn the first half and keep every other row M=readMCMC(file=file.path(workspace,'MCMC.txt'),burnFactor=0.5,slimFactor=2) dim(M)# Create Monte Carlo samples and write them to file n=4000 sim=data.frame(p1=rnorm(n),p2=rlnorm(n),p3=runif(n)) workspace=tempdir() write.table(sim,file=file.path(workspace,'MCMC.txt'),row.names=FALSE) # Read file, burn the first half and keep every other row M=readMCMC(file=file.path(workspace,'MCMC.txt'),burnFactor=0.5,slimFactor=2) dim(M)
Creates a new instance of a 'remnantErrorModel' object
remnantErrorModel( fname = "Config_RemnantSigma.txt", funk = "Linear", par = list(parameter("g1", 1, prior.dist = "FlatPrior+"), parameter("g2", 0.1, prior.dist = "FlatPrior+")) )remnantErrorModel( fname = "Config_RemnantSigma.txt", funk = "Linear", par = list(parameter("g1", 1, prior.dist = "FlatPrior+"), parameter("g2", 0.1, prior.dist = "FlatPrior+")) )
fname |
Character, configuration file name. |
funk |
Character, function f used in remnant sdev = f(Ysim). Available: 'Constant', 'Proportional', 'Linear' (default), 'Exponential', 'Gaussian'. |
par |
list of parameter objects, parameters of the function above. respectively, npar= 1,1,2,3,3 |
An object of class 'remnantErrorModel'.
r <- remnantErrorModel()r <- remnantErrorModel()
Creates a new instance of a 'residualOptions' object
residualOptions( fname = "Config_Residuals.txt", result.fname = "Results_Residuals.txt" )residualOptions( fname = "Config_Residuals.txt", result.fname = "Results_Residuals.txt" )
fname |
Character, configuration file name. |
result.fname |
Character, result file name. |
An object of class 'residualOptions'.
r <- residualOptions()r <- residualOptions()
Perform a single run of a model, using the initial values specified for the model's parameters.
runModel( workspace, mod, X, na.value = -666.666, run = TRUE, preClean = FALSE, dir.exe = .BAM_PATH, name.exe = "BaM", stout = "", Inputs_fname = "Config_Inputs.txt", X_fname = "X.txt", Y_fname = "Y.txt" )runModel( workspace, mod, X, na.value = -666.666, run = TRUE, preClean = FALSE, dir.exe = .BAM_PATH, name.exe = "BaM", stout = "", Inputs_fname = "Config_Inputs.txt", X_fname = "X.txt", Y_fname = "Y.txt" )
workspace |
Character, directory where config and result files are stored. workspace = tempdir() is recommended. |
mod |
model object, the model to be run. |
X |
data frame, containing the inputs of the model |
na.value |
numeric, value used by BaM to denote impossible runs, that will be changed to NA in RBaM. |
run |
Logical, run the model? if FALSE, just write config files and returns NULL. |
preClean |
Logical, start by cleaning up workspace? Be careful, this will delete all files in the workspace, including old results! |
dir.exe |
Character, directory where BaM executable stands. |
name.exe |
Character, name of the executable without extension ('BaM' by default). |
stout |
Character string, standard output (see ?system2). In particular, stout="" (default) shows BaM messages in the console, stout=NULL discards BaM messages, stout='log.txt' saves BaM messages in file "log.txt". |
Inputs_fname |
Character, name of configuration file used to specify the format of X. |
X_fname |
Character, name of file containing a copy of X. |
Y_fname |
Character, name of file where simulations are written. |
A data frame containing the outputs simulated by the model.
# Rating curve model - see https://github.com/BaM-tools/RBaM # Parameters of the low flow section control: activation stage k, coefficient a and exponent c k1=parameter(name='k1',init=-0.5) a1=parameter(name='a1',init=50) c1=parameter(name='c1',init=1.5) # Parameters of the high flow channel control: activation stage k, coefficient a and exponent c k2=parameter(name='k2',init=1) a2=parameter(name='a2',init=100) c2=parameter(name='c2',init=1.67) # Define control matrix: columns are controls, rows are stage ranges. controlMatrix=rbind(c(1,0),c(0,1)) # Stitch it all together into a model object M=model(ID='BaRatin', nX=1,nY=1, # number of input/output variables par=list(k1,a1,c1,k2,a2,c2), # list of model parameters xtra=xtraModelInfo(object=controlMatrix)) # use xtraModelInfo() to pass the control matrix # Define the model input X=data.frame(stage=seq(-1,7,0.1)) # Write the config files used to run the model. To actually run it, # use run=TRUE, but BaM executable needs to be downloaded first (use downloadBaM()) runModel(workspace=tempdir(),mod=M,X=X,run=FALSE)# Rating curve model - see https://github.com/BaM-tools/RBaM # Parameters of the low flow section control: activation stage k, coefficient a and exponent c k1=parameter(name='k1',init=-0.5) a1=parameter(name='a1',init=50) c1=parameter(name='c1',init=1.5) # Parameters of the high flow channel control: activation stage k, coefficient a and exponent c k2=parameter(name='k2',init=1) a2=parameter(name='a2',init=100) c2=parameter(name='c2',init=1.67) # Define control matrix: columns are controls, rows are stage ranges. controlMatrix=rbind(c(1,0),c(0,1)) # Stitch it all together into a model object M=model(ID='BaRatin', nX=1,nY=1, # number of input/output variables par=list(k1,a1,c1,k2,a2,c2), # list of model parameters xtra=xtraModelInfo(object=controlMatrix)) # use xtraModelInfo() to pass the control matrix # Define the model input X=data.frame(stage=seq(-1,7,0.1)) # Write the config files used to run the model. To actually run it, # use run=TRUE, but BaM executable needs to be downloaded first (use downloadBaM()) runModel(workspace=tempdir(),mod=M,X=X,run=FALSE)
Creates a new instance of a 'runOptions' object
runOptions( fname = "Config_RunOptions.txt", doMCMC = TRUE, doSummary = TRUE, doResiduals = TRUE, doPrediction = FALSE )runOptions( fname = "Config_RunOptions.txt", doMCMC = TRUE, doSummary = TRUE, doResiduals = TRUE, doPrediction = FALSE )
fname |
Character, configuration file name. |
doMCMC |
logical, do MCMC sampling? |
doSummary |
logical, do MCMC summarizing? |
doResiduals |
logical, do residuals analysis? |
doPrediction |
logical, do prediction experiments? |
An object of class 'runOptions'.
o <- runOptions()o <- runOptions()
Stage-discharge gaugings from the hydrometric station 'the Ardèche River at Sauze-St-Martin'. See https://en.wikipedia.org/wiki/Ardèche_(river) for a description of the river See https://hal.science/hal-00934237 for an article using this dataset
SauzeGaugingsSauzeGaugings
A data frame with 38 rows and 3 variables:
Stage (m)
Discharge (m3/s)
Discharge uncertainty (m3/s) expressed as a standard deviation
Set the initial values of a list of parameters.
setInitPar(parameters, values)setInitPar(parameters, values)
parameters |
List of parameter objects, parameters whose initial values are to be set. |
values |
Numeric vector, initial values. |
a list of parameter objects equal to the input list 'parameters', except for their initial values.
ps <- list(parameter(name='par1',init=0),parameter(name='par2',init=1)) ps=setInitPar(ps,c(10,100)) print(ps)ps <- list(parameter(name='par1',init=0),parameter(name='par2',init=1)) ps=setInitPar(ps,c(10,100)) print(ps)
Set path to BaM executable
setPathToBaM(dir.exe, quiet = FALSE)setPathToBaM(dir.exe, quiet = FALSE)
dir.exe |
character string, folder where BaM executable is located. A NULL value resets BaM directory to 'unknown' by removing the config folder where it was stored on your computer (use 'tools::R_user_dir(package="RBaM",which="config")' to locate this folder). |
quiet |
logical, if TRUE, suppress status messages. |
nothing - just write a config file.
setPathToBaM(dir.exe=tempdir())setPathToBaM(dir.exe=tempdir())
Run BaM to estimate a BaRatin-SPD model. The following assumptions hold:
Only parameters b/k's and a's can be variable. Exponents c's are stable.
Incremental changes affecting parameters b/k's are additive, while incremental changes affecting parameters a's are multiplicative.
The prior distribution for parameters b/k's and associated incremental changes has to be 'Gaussian'.
The prior distribution for parameters a's and associated incremental changes has to be 'LogNormal'.
SPD_estimate( workspace, controlMatrix, pars, bVAR, aVAR, deltaPars, periods, H, Q, uQ = 0 * Q, nPeriods = lapply(periods, max), BaRatinFlavor = "BaRatinBAC", remnant = remnantErrorModel(funk = "Linear", par = list(parameter("g1", 1, "LogNormal", c(0, 10)), parameter("g2", 0.1, "LogNormal", c(log(0.1), 10)))), mcmcOpt = mcmcOptions(), mcmcCook = mcmcCooking() )SPD_estimate( workspace, controlMatrix, pars, bVAR, aVAR, deltaPars, periods, H, Q, uQ = 0 * Q, nPeriods = lapply(periods, max), BaRatinFlavor = "BaRatinBAC", remnant = remnantErrorModel(funk = "Linear", par = list(parameter("g1", 1, "LogNormal", c(0, 10)), parameter("g2", 0.1, "LogNormal", c(log(0.1), 10)))), mcmcOpt = mcmcOptions(), mcmcCook = mcmcCooking() )
workspace |
Character, directory where config and result files are stored. |
controlMatrix |
Integer matrix, control matrix, dimension nControl*nControl. |
pars |
list of parameter objects, parameters of the model. For VAR parameters, they will be interpreted as the prior for period 1. |
bVAR |
Logical vector, size nControl. bVAR[i]=TRUE means that the b/k parameter of control i is variable, otherwise it is stable. |
aVAR |
Logical vector, size nControl. aVaR[i]=TRUE means that the a parameter of control i is variable, otherwise it is stable. |
deltaPars |
list, prior parameters of incremental changes for each VAR parameter. deltaPars should be a named list, with the names corresponding to the names of the parameters that have been declared variable. Each element of deltaPars is then a numeric vector of size 2. For b/k's, the 2 values are the mean/sd of the Gaussian prior for ADDITIVE incremental changes. For a's, the 2 parameters are the meanlog/sdlog of the LogNormal prior for MULTIPLICATIVE incremental changes. |
periods |
list, period index for each VAR parameter. periods should be a named list as previously. Each element of the list is an integer vector (starting at 1) with same length as the calibration data. Periods do not need to be the same for all VAR parameters. |
H |
numeric vector, gauging stages. |
Q |
numeric vector, gauging discharges. |
uQ |
numeric vector, gauging discharge uncertainties. |
nPeriods |
list, number of periods for each VAR parameter. nPeriods should be a named list as deltaPars and periods. In general and by default, nPeriods[[i]] is just the max of periods[[i]], but this is not compulsory: there could be one or several additional periods with no gaugings. |
BaRatinFlavor |
character, either 'BaRatinBAC' (default) or 'BaRatin' (the original k-a-c parameterization). It is in general easier to specify priors on changes affecting b's than k's, hence the default choice. However, 'BaRatinBAC' requires some numerical resolution and it is hence a bit slower and more prone to failures than 'BaRatin'. |
remnant |
remnantErrorModel object, by default the structural standard deviation varies as an affine function of simulated discharges, with very wide priors on coefficients g1 and g2. |
mcmcOpt |
mcmcOptions object, MCMC options passed to BaM. |
mcmcCook |
mcmcCooking object, MCMC cooking options (burn and slice) passed to BaM. |
A data frame containing the MCMC simulations performed by BaM.
# Calibration data H=MeyrasGaugings$h Q=MeyrasGaugings$Q uQ=MeyrasGaugings$uQ # Control matrix controlMatrix=rbind(c(1,0,0),c(0,1,0),c(0,1,1)) # Declare variable parameters. bVAR=c(TRUE, TRUE, FALSE) # b's for first 2 controls (k1 and k2) are VAR aVAR=c(TRUE, FALSE, FALSE) # a for first control (a1) is VAR # Define priors. b1=parameter(name='b1',init=-0.6,prior.dist='Gaussian',prior.par=c(-0.6,0.5)) a1=parameter(name='a1',init=exp(2.65),prior.dist='LogNormal',prior.par=c(2.65,0.35)) c1=parameter(name='c1',init=1.5,prior.dist='Gaussian',prior.par=c(1.5,0.025)) b2=parameter(name='b2',init=0,prior.dist='Gaussian',prior.par=c(-0.6,0.5)) a2=parameter(name='a2',init=exp(3.28),prior.dist='LogNormal',prior.par=c(3.28,0.33)) c2=parameter(name='c2',init=1.67,prior.dist='Gaussian',prior.par=c(1.67,0.025)) b3=parameter(name='b3',init=1.2,prior.dist='Gaussian',prior.par=c(1.2,0.2)) a3=parameter(name='a3',init=exp(3.48),prior.dist='LogNormal',prior.par=c(3.46,0.38)) c3=parameter(name='c3',init=1.67,prior.dist='Gaussian',prior.par=c(1.67,0.025)) pars=list(b1,a1,c1,b2,a2,c2,b3,a3,c3) # Define properties of VAR parameters. deltaPars=list(b1=c(0,0.25),a1=c(0,0.2),b2=c(0,0.5)) periods=list(b1=MeyrasGaugings$Period,a1=c(rep(1,49),rep(2,55)),b2=MeyrasGaugings$Period) # Run BaM and estimate SPD parameters mcmcOpt=mcmcOptions(nAdapt=20,nCycles=25) # only few iterations so that the example runs fast. mcmc=SPD_estimate(workspace=tempdir(),controlMatrix=controlMatrix,pars=pars, bVAR=bVAR,aVAR=aVAR,deltaPars=deltaPars,periods=periods, H=H,Q=Q,uQ=uQ,mcmcOpt=mcmcOpt)# Calibration data H=MeyrasGaugings$h Q=MeyrasGaugings$Q uQ=MeyrasGaugings$uQ # Control matrix controlMatrix=rbind(c(1,0,0),c(0,1,0),c(0,1,1)) # Declare variable parameters. bVAR=c(TRUE, TRUE, FALSE) # b's for first 2 controls (k1 and k2) are VAR aVAR=c(TRUE, FALSE, FALSE) # a for first control (a1) is VAR # Define priors. b1=parameter(name='b1',init=-0.6,prior.dist='Gaussian',prior.par=c(-0.6,0.5)) a1=parameter(name='a1',init=exp(2.65),prior.dist='LogNormal',prior.par=c(2.65,0.35)) c1=parameter(name='c1',init=1.5,prior.dist='Gaussian',prior.par=c(1.5,0.025)) b2=parameter(name='b2',init=0,prior.dist='Gaussian',prior.par=c(-0.6,0.5)) a2=parameter(name='a2',init=exp(3.28),prior.dist='LogNormal',prior.par=c(3.28,0.33)) c2=parameter(name='c2',init=1.67,prior.dist='Gaussian',prior.par=c(1.67,0.025)) b3=parameter(name='b3',init=1.2,prior.dist='Gaussian',prior.par=c(1.2,0.2)) a3=parameter(name='a3',init=exp(3.48),prior.dist='LogNormal',prior.par=c(3.46,0.38)) c3=parameter(name='c3',init=1.67,prior.dist='Gaussian',prior.par=c(1.67,0.025)) pars=list(b1,a1,c1,b2,a2,c2,b3,a3,c3) # Define properties of VAR parameters. deltaPars=list(b1=c(0,0.25),a1=c(0,0.2),b2=c(0,0.5)) periods=list(b1=MeyrasGaugings$Period,a1=c(rep(1,49),rep(2,55)),b2=MeyrasGaugings$Period) # Run BaM and estimate SPD parameters mcmcOpt=mcmcOptions(nAdapt=20,nCycles=25) # only few iterations so that the example runs fast. mcmc=SPD_estimate(workspace=tempdir(),controlMatrix=controlMatrix,pars=pars, bVAR=bVAR,aVAR=aVAR,deltaPars=deltaPars,periods=periods, H=H,Q=Q,uQ=uQ,mcmcOpt=mcmcOpt)
This function computes the properties of VAR parameters (prior means, sds and correlation) given the prior for period 1 and the prior for incremental changes.
SPD_getVARproperties(p1Par, deltaPar, nPeriod)SPD_getVARproperties(p1Par, deltaPar, nPeriod)
p1Par |
numeric vector of length 2, parameters of the prior for period 1. |
deltaPar |
numeric vector of length 2, parameters of the prior for incremental changes. |
nPeriod |
integer, number of periods. |
A list with the following components:
mean |
numeric vector of length nPeriod, prior means. |
sd |
numeric vector of length nPeriod, prior standard deviations. |
cor |
matrix of dim nPeriod*nPeriod, prior correlation matrix. |
res=SPD_getVARproperties(p1Par=c(2,0.1),deltaPar=c(0,0.3),nPeriod=25) image(res$cor)res=SPD_getVARproperties(p1Par=c(2,0.1),deltaPar=c(0,0.3),nPeriod=25) image(res$cor)
Convert an object of class 'dataset' into a ready-to-write vector of string
## S3 method for class 'dataset' toString(x, ...)## S3 method for class 'dataset' toString(x, ...)
x |
dataset object, object to be converted. |
... |
Optional arguments. |
A string ready to be printed or written.
X=data.frame(input1=rnorm(100),input2=rnorm(100)) Y=data.frame(output=X$input1+0.8*X$input2+0.1*rnorm(100)) workspace=tempdir() d <- dataset(X=X,Y=Y,data.dir=workspace) toString(d)X=data.frame(input1=rnorm(100),input2=rnorm(100)) Y=data.frame(output=X$input1+0.8*X$input2+0.1*rnorm(100)) workspace=tempdir() d <- dataset(X=X,Y=Y,data.dir=workspace) toString(d)
Convert an object of class 'mcmcCooking' into a ready-to-write vector of string
## S3 method for class 'mcmcCooking' toString(x, ...)## S3 method for class 'mcmcCooking' toString(x, ...)
x |
mcmcCooking object, object to be converted. |
... |
Optional arguments. |
A string ready to be printed or written.
toString(mcmcCooking())toString(mcmcCooking())
Convert an object of class 'mcmcOptions' into a ready-to-write vector of string
## S3 method for class 'mcmcOptions' toString(x, ...)## S3 method for class 'mcmcOptions' toString(x, ...)
x |
mcmcOptions object, object to be converted. |
... |
Optional arguments. |
A string ready to be printed or written.
toString(mcmcOptions())toString(mcmcOptions())
Convert an object of class 'mcmcSummary' into a ready-to-write vector of string
## S3 method for class 'mcmcSummary' toString(x, ...)## S3 method for class 'mcmcSummary' toString(x, ...)
x |
mcmcSummary object, object to be converted. |
... |
Optional arguments. |
A string ready to be printed or written.
toString(mcmcSummary())toString(mcmcSummary())
Convert an object of class 'model' into a ready-to-write vector of string
## S3 method for class 'model' toString(x, ...)## S3 method for class 'model' toString(x, ...)
x |
model object, object to be converted. |
... |
Optional arguments. |
A string ready to be printed or written.
toString(model())toString(model())
Convert an object of class 'parameter' into a ready-to-write vector of string
## S3 method for class 'parameter' toString(x, ...)## S3 method for class 'parameter' toString(x, ...)
x |
parameter object, object to be converted. |
... |
Optional arguments. |
A string ready to be printed or written.
p <- parameter(name='par',init=0,prior.dist='Gaussian',prior.par=c(0,1)) toString(p)p <- parameter(name='par',init=0,prior.dist='Gaussian',prior.par=c(0,1)) toString(p)
Convert an object of class 'parameter_VAR' into a ready-to-write vector of string
## S3 method for class 'parameter_VAR' toString(x, ...)## S3 method for class 'parameter_VAR' toString(x, ...)
x |
parameter_VAR object, object to be converted. |
... |
Optional arguments. |
A string ready to be printed or written.
X=data.frame(input1=rnorm(100),input2=rnorm(100)) Y=data.frame(output=X$input1+0.8*X$input2+0.1*rnorm(100)) VAR.indx=data.frame(indx=c(rep(1,50),rep(2,50))) workspace=tempdir() d <- dataset(X=X,Y=Y,data.dir=workspace,VAR.indx=VAR.indx) p <- parameter_VAR(name='par',index='indx',d=d, init=c(-1,1,2), prior.dist=c('Gaussian','FlatPrior','Triangle'), prior.par=list(c(-1,1),NULL,c(2,0,5))) toString(p)X=data.frame(input1=rnorm(100),input2=rnorm(100)) Y=data.frame(output=X$input1+0.8*X$input2+0.1*rnorm(100)) VAR.indx=data.frame(indx=c(rep(1,50),rep(2,50))) workspace=tempdir() d <- dataset(X=X,Y=Y,data.dir=workspace,VAR.indx=VAR.indx) p <- parameter_VAR(name='par',index='indx',d=d, init=c(-1,1,2), prior.dist=c('Gaussian','FlatPrior','Triangle'), prior.par=list(c(-1,1),NULL,c(2,0,5))) toString(p)
Convert an object of class 'prediction' into a ready-to-write vector of string
## S3 method for class 'prediction' toString(x, ...)## S3 method for class 'prediction' toString(x, ...)
x |
prediction object, object to be converted. |
... |
Optional arguments. |
A string ready to be printed or written.
pred=prediction(X=twoPopulations[,1:3],spagFiles=c('P1.spag','P2.spag')) toString(pred)pred=prediction(X=twoPopulations[,1:3],spagFiles=c('P1.spag','P2.spag')) toString(pred)
Convert an object of class 'remnantErrorModel' into a ready-to-write vector of string
## S3 method for class 'remnantErrorModel' toString(x, ...)## S3 method for class 'remnantErrorModel' toString(x, ...)
x |
remnantErrorModel object, object to be converted. |
... |
Optional arguments. |
A string ready to be printed or written.
toString(remnantErrorModel())toString(remnantErrorModel())
Convert an object of class 'residualOptions' into a ready-to-write vector of string
## S3 method for class 'residualOptions' toString(x, ...)## S3 method for class 'residualOptions' toString(x, ...)
x |
residualOptions object, object to be converted. |
... |
Optional arguments. |
A string ready to be printed or written.
toString(residualOptions())toString(residualOptions())
Convert an object of class 'runOptions' into a ready-to-write vector of string
## S3 method for class 'runOptions' toString(x, ...)## S3 method for class 'runOptions' toString(x, ...)
x |
runOptions object, object to be converted. |
... |
Optional arguments. |
A string ready to be printed or written.
toString(runOptions())toString(runOptions())
2DO (adapt from STooDs): Generate pdf report files summarizing mcmc samples
tracePlot(sim, ylab = "values", keep = NULL, col = "black", psize = 0.5)tracePlot(sim, ylab = "values", keep = NULL, col = "black", psize = 0.5)
sim |
vector or matrix or data frame, MCMC simulations |
ylab |
Character, label of y-axis to be used if sim has no names |
keep |
Integer vector, indices of samples to be kept in cooked MCMC sample |
col |
Color |
psize |
Numeric, point size |
tracePlot
returns a trace plot ggplot (or a list thereof if several columns in sim)
A ggplot (or a list thereof if several columns in sim)
# Create Monte Carlo samples n=1000 sim=data.frame(p1=rnorm(n),p2=rlnorm(n),p3=runif(n)) # create trace plot for each component figures=tracePlot(sim)# Create Monte Carlo samples n=1000 sim=data.frame(p1=rnorm(n),p2=rlnorm(n),p3=runif(n)) # create trace plot for each component figures=tracePlot(sim)
Size of two populations of the same species put in two different environments, as a function of time and local temperature.
Input variables: time t, temperature at site 1 T1, temperature at site 2 T2.
Output variables: population size at site 1 P1, population size at site 2 P2.
Data are synthetically generated from a logistic model.
twoPopulationstwoPopulations
An object of class data.frame with 101 rows and 5 columns.
returns a violinplot ggplot
violinPlot(sim, ylab = "values", col = "black")violinPlot(sim, ylab = "values", col = "black")
sim |
vector or matrix or data frame, MCMC simulations |
ylab |
Character, label of y-axis |
col |
Color |
A ggplot
# Create Monte Carlo samples n=1000 sim=data.frame(p1=rnorm(n),p2=rlnorm(n),p3=runif(n)) # create violin plot comparing all components figure=violinPlot(sim)# Create Monte Carlo samples n=1000 sim=data.frame(p1=rnorm(n),p2=rlnorm(n),p3=runif(n)) # create violin plot comparing all components figure=violinPlot(sim)
Write input data of the prediction into files
writePredInputs(o)writePredInputs(o)
o |
prediction object |
nothing - just write to files.
temp=tempdir() pred=prediction(X=twoPopulations[,1:3],spagFiles=c('P1.spag','P2.spag'), data.dir=temp) writePredInputs(pred)temp=tempdir() pred=prediction(X=twoPopulations[,1:3],spagFiles=c('P1.spag','P2.spag'), data.dir=temp) writePredInputs(pred)
Creates a new instance of a 'xtraModelInfo' object containing extra model information
xtraModelInfo(fname = "Config_Xtra.txt", object = NULL)xtraModelInfo(fname = "Config_Xtra.txt", object = NULL)
fname |
Character, configuration file name. |
object |
any R object containing xtra model info - typically a list of stuff. The content and meaning of 'object' is completely model-specific. |
An object of class 'xtraModelInfo'.
x <- xtraModelInfo()x <- xtraModelInfo()