Last data update: 2014.03.03

 R: Running the BayesFit optimisation
runBayesFitR Documentation

Running the BayesFit optimisation

Description

"runBayesFit" defines the prior function and runs the BayesFit estimation

Usage

runBayesFit(opts)

Arguments

opts

list with entries as explained below. Options set – defines the problem and sets some parameters to control the MCMC algorithm. model: List of model parameters - to estimate. The parameter objects must each have a 'value' attribute containing the parameter's numerical value. estimate_params: list. List of parameters to estimate, all of which must also be listed in 'options$model$parameters'. initial_values: list of float, optional. Starting values for parameters to estimate. If omitted, will use the nominal values from 'options$model$parameters' step_fn: callable f(output), optional. User callback, called on every MCMC iteration. likelihood_fn: callable f(output, position). User likelihood function. prior_fn: callable f(output, position), optional. User prior function. If omitted, a flat prior will be used. nsteps: int. Number of MCMC iterations to perform. use_hessian: logical, optional. Wheter to use the Hessian to guide the walk. Defaults to FALSE. rtol: float or list of float, optional. Relative tolerance for ode solver. atol: float or list of float, optional. Absolute tolerance for ode solver. norm_step_size: float, optional. MCMC step size. Defaults to a reasonable value. hessian_period: int, optional. Number of MCMC steps between Hessian recalculations. Defaults to a reasonable but fairly large value, as Hessian calculation is expensive. hessian_scale: float, optional. Scaling factor used in generating Hessian-guided steps. Defaults to a reasonable value. sigma_adj_interval: int, optional. How often to adjust 'output$sig_value' while annealing to meet 'accept_rate_target'. Defaults to a reasonable value. anneal_length: int, optional. Length of initial "burn-in" annealing period. Defaults to 10 'nsteps', or if 'use_hessian' is TRUE, to 'hessian_period' (i.e. anneal until first hessian is calculated) T_init: float, optional. Initial temperature for annealing. Defaults to a resonable value. accept_rate_target: float, optional. Desired acceptance rate during annealing. Defaults to a reasonable value. See also 'sigma_adj_interval' above. sigma_max: float, optional. Maximum value for 'output$sig_value'. Defaults to a resonable value. sigma_min: float, optional. Minimum value for 'output$sig_value'. Defaults to a resonable value. sigma_step: float, optional. Increment for 'output$sig_value' adjustments. Defaults to a resonable value. To eliminate adaptive step size, set sigma_step to 1. thermo_temp: float in the range [0,1], optional. Temperature for thermodynamic integration support. Used to scale likelihood when calculating the posterior value. Defaults to 1, i.e. no effect.

Value

The output after the optimisation is finished - a list with entries as explained in 'Arguments'.

Examples


data("simpleExample", package="MEIGOR")
initial_pars = createLBodeContPars(model, LB_n=1, LB_k=0.1, LB_tau=0.01, UB_n=5, UB_k=0.9, UB_tau=10, random=TRUE)
simData = plotLBodeFitness(cnolist, model, initial_pars, reltol=1e-05, atol=1e-03, maxStepSize=0.01)

f_bayesFit <- function(position, params=initial_pars, exp_var=opts$exp_var) {
	# convert from log
	params$parValues = 10^position
	ysim = getLBodeDataSim(cnolist=cnolist, model=model,
	ode_parameters=params)
	data_as_vec = unlist(cnolist$valueSignals)
	sim_as_vec = unlist(ysim)
	# set nan (NAs) to 0
	sim_as_vec[is.na(sim_as_vec)] = 0
	sim_as_vec[is.nan(sim_as_vec)]= 0
	return(sum((data_as_vec-sim_as_vec)^2/(2*exp_var^2)))
}

prior_mean = log10(initial_pars$parValues)
prior_var = 10

opts <- list("model"=NULL, "estimate_params"=NULL,"initial_values"=NULL,
"tspan"=NULL, "step_fn"=NULL, "likelihood_fn"=NULL,
"prior_fn"=NULL, "nsteps"=NULL, "use_hessian"=FALSE,
"rtol"=NULL, "atol"=NULL, "norm_step_size"=0.75,
"hessian_period"=25000, "hessian_scale"=0.085,
"sigma_adj_interval"=NULL, "anneal_length"=NULL,
"T_init"=10, "accept_rate_target"=0.3, "sigma_max"=1,
"sigma_min"=0.25, "sigma_step"=0.125, "thermo_temp"=1, "seed"=NULL)
opts$nsteps = 2000
opts$likelihood_fn = f_bayesFit
opts$use_hessian = TRUE
opts$hessian_period = opts$nsteps/10
opts$model = list(parameters=list(name=initial_pars$parNames,
value=initial_pars$parValues))
opts$estimate_params = initial_pars$parValues
opts$exp_var = 0.01

res = runBayesFit(opts)

initial_pars$parValues = cur_params(output=res, options=opts)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(MEIGOR)
Loading required package: Rsolnp
Loading required package: snowfall
Loading required package: snow
Loading required package: CNORode
Loading required package: CellNOptR
Loading required package: RBGL
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'parallel'

The following objects are masked from 'package:snow':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
    parCapply, parLapply, parRapply, parSapply, splitIndices,
    stopCluster


Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:snow':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parRapply, parSapply

The following objects are masked from 'package:stats':

    IQR, mad, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, cbind, colnames, do.call, duplicated, eval, evalq,
    get, grep, grepl, intersect, is.unsorted, lapply, lengths, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, rank,
    rbind, rownames, sapply, setdiff, sort, table, tapply, union,
    unique, unsplit

Loading required package: hash
hash-2.2.6 provided by Decision Patterns

Loading required package: ggplot2
Loading required package: RCurl
Loading required package: bitops
Loading required package: Rgraphviz
Loading required package: grid
Loading required package: XML

Attaching package: 'XML'

The following object is masked from 'package:graph':

    addNode

Loading required package: genalg

Attaching package: 'CNORode'

The following object is masked from 'package:stats':

    simulate

Loading required package: deSolve

Attaching package: 'deSolve'

The following object is masked from 'package:graphics':

    matplot

> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/MEIGOR/runBayesFit.Rd_%03d_medium.png", width=480, height=480)
> ### Name: runBayesFit
> ### Title: Running the BayesFit optimisation
> ### Aliases: runBayesFit
> 
> ### ** Examples
> 
> 
> data("simpleExample", package="MEIGOR")
> initial_pars = createLBodeContPars(model, LB_n=1, LB_k=0.1, LB_tau=0.01, UB_n=5, UB_k=0.9, UB_tau=10, random=TRUE)
> simData = plotLBodeFitness(cnolist, model, initial_pars, reltol=1e-05, atol=1e-03, maxStepSize=0.01)
> 
> f_bayesFit <- function(position, params=initial_pars, exp_var=opts$exp_var) {
+ 	# convert from log
+ 	params$parValues = 10^position
+ 	ysim = getLBodeDataSim(cnolist=cnolist, model=model,
+ 	ode_parameters=params)
+ 	data_as_vec = unlist(cnolist$valueSignals)
+ 	sim_as_vec = unlist(ysim)
+ 	# set nan (NAs) to 0
+ 	sim_as_vec[is.na(sim_as_vec)] = 0
+ 	sim_as_vec[is.nan(sim_as_vec)]= 0
+ 	return(sum((data_as_vec-sim_as_vec)^2/(2*exp_var^2)))
+ }
> 
> prior_mean = log10(initial_pars$parValues)
> prior_var = 10
> 
> opts <- list("model"=NULL, "estimate_params"=NULL,"initial_values"=NULL,
+ "tspan"=NULL, "step_fn"=NULL, "likelihood_fn"=NULL,
+ "prior_fn"=NULL, "nsteps"=NULL, "use_hessian"=FALSE,
+ "rtol"=NULL, "atol"=NULL, "norm_step_size"=0.75,
+ "hessian_period"=25000, "hessian_scale"=0.085,
+ "sigma_adj_interval"=NULL, "anneal_length"=NULL,
+ "T_init"=10, "accept_rate_target"=0.3, "sigma_max"=1,
+ "sigma_min"=0.25, "sigma_step"=0.125, "thermo_temp"=1, "seed"=NULL)
> opts$nsteps = 2000
> opts$likelihood_fn = f_bayesFit
> opts$use_hessian = TRUE
> opts$hessian_period = opts$nsteps/10
> opts$model = list(parameters=list(name=initial_pars$parNames,
+ value=initial_pars$parValues))
> opts$estimate_params = initial_pars$parValues
> opts$exp_var = 0.01
> 
> res = runBayesFit(opts)
Running BayesFit 
========================= 
iter=0 sigma=1 T=10 acc=1 lkl=13062.5 prior=0.028 post=13062.53
iter=20 sigma=1 T=5.939305 acc=0.952 lkl=13062.93 prior=0.644 post=13063.58
iter=40 sigma=1 T=3.710748 acc=0.902 lkl=13062.5 prior=0.885 post=13063.39
iter=60 sigma=1 T=2.48769 acc=0.934 lkl=13062.59 prior=0.833 post=13063.43
iter=80 sigma=1 T=1.816462 acc=0.951 lkl=13062.48 prior=2.011 post=13064.49
iter=100 sigma=1 T=1.448084 acc=0.871 lkl=13059.76 prior=1.907 post=13061.66
iter=120 sigma=1 T=1.245914 acc=0.802 lkl=9046.33 prior=2.011 post=9048.34
iter=140 sigma=1 T=1.13496 acc=0.709 lkl=8300.11 prior=2.13 post=8302.24
iter=160 sigma=1 T=1.074068 acc=0.634 lkl=8201.99 prior=2.377 post=8204.37
iter=180 sigma=1 T=1.040649 acc=0.564 lkl=8201.99 prior=2.377 post=8204.37
iter=200 sigma=1 T=1.022988 acc=0.507 lkl=8201.99 prior=2.377 post=8204.37
iter=220 sigma=1 T=1.022988 acc=0.466 lkl=8199.86 prior=2.376 post=8202.23
iter=240 sigma=1 T=1.022988 acc=0.436 lkl=8192.55 prior=2.372 post=8194.92
iter=260 sigma=1 T=1.022988 acc=0.402 lkl=8192.55 prior=2.372 post=8194.92
iter=280 sigma=1 T=1.022988 acc=0.374 lkl=8192.55 prior=2.372 post=8194.92
iter=300 sigma=1 T=1.022988 acc=0.349 lkl=8192.55 prior=2.372 post=8194.92
iter=320 sigma=1 T=1.022988 acc=0.327 lkl=8192.55 prior=2.372 post=8194.92
iter=340 sigma=1 T=1.022988 acc=0.311 lkl=8191.6 prior=2.362 post=8193.96
iter=360 sigma=0.875 T=1.022988 acc=0.294 lkl=8191.6 prior=2.362 post=8193.96
iter=380 sigma=0.75 T=1.022988 acc=0.278 lkl=8191.6 prior=2.362 post=8193.96
iter=400 sigma=0.625 T=1.022988 acc=0.267 lkl=8189.87 prior=2.355 post=8192.22
iter=420 sigma=0.5 T=1.022988 acc=0.257 lkl=8190.76 prior=2.352 post=8193.12
iter=440 sigma=0.375 T=1.022988 acc=0.247 lkl=8188.97 prior=2.276 post=8191.24
iter=460 sigma=0.25 T=1.022988 acc=0.243 lkl=8187.72 prior=2.187 post=8189.91
iter=480 sigma=0.25 T=1.022988 acc=0.239 lkl=8174.36 prior=2.239 post=8176.6
iter=500 sigma=0.25 T=1.022988 acc=0.232 lkl=8151.11 prior=2.253 post=8153.36
iter=520 sigma=0.25 T=1.022988 acc=0.223 lkl=8151.11 prior=2.253 post=8153.36
iter=540 sigma=0.25 T=1.022988 acc=0.214 lkl=8151.11 prior=2.253 post=8153.36
iter=560 sigma=0.25 T=1.022988 acc=0.207 lkl=8151.11 prior=2.253 post=8153.36
iter=580 sigma=0.25 T=1.022988 acc=0.2 lkl=8151.11 prior=2.253 post=8153.36
iter=600 sigma=0.25 T=1.022988 acc=0.193 lkl=8151.11 prior=2.253 post=8153.36
iter=620 sigma=0.25 T=1.022988 acc=0.187 lkl=8151.11 prior=2.253 post=8153.36
iter=640 sigma=0.25 T=1.022988 acc=0.181 lkl=8151.11 prior=2.253 post=8153.36
iter=660 sigma=0.25 T=1.022988 acc=0.177 lkl=8151.33 prior=2.199 post=8153.53
iter=680 sigma=0.25 T=1.022988 acc=0.173 lkl=8151.26 prior=2.126 post=8153.39
iter=700 sigma=0.25 T=1.022988 acc=0.168 lkl=8151.26 prior=2.126 post=8153.39
iter=720 sigma=0.25 T=1.022988 acc=0.165 lkl=8151.58 prior=2.169 post=8153.75
iter=740 sigma=0.25 T=1.022988 acc=0.161 lkl=8151.58 prior=2.169 post=8153.75
iter=760 sigma=0.25 T=1.022988 acc=0.156 lkl=8151.58 prior=2.169 post=8153.75
iter=780 sigma=0.25 T=1.022988 acc=0.152 lkl=8151.58 prior=2.169 post=8153.75
iter=800 sigma=0.25 T=1.022988 acc=0.149 lkl=8151.58 prior=2.169 post=8153.75
iter=820 sigma=0.25 T=1.022988 acc=0.146 lkl=8151.57 prior=2.087 post=8153.66
iter=840 sigma=0.25 T=1.022988 acc=0.143 lkl=8151.57 prior=2.087 post=8153.66
iter=860 sigma=0.25 T=1.022988 acc=0.139 lkl=8151.57 prior=2.087 post=8153.66
iter=880 sigma=0.25 T=1.022988 acc=0.137 lkl=8154.42 prior=2.061 post=8156.48
iter=900 sigma=0.25 T=1.022988 acc=0.134 lkl=8154.42 prior=2.061 post=8156.48
iter=920 sigma=0.25 T=1.022988 acc=0.131 lkl=8154.42 prior=2.061 post=8156.48
iter=940 sigma=0.25 T=1.022988 acc=0.129 lkl=8154.42 prior=2.061 post=8156.48
iter=960 sigma=0.25 T=1.022988 acc=0.127 lkl=8151.22 prior=2.037 post=8153.25
iter=980 sigma=0.25 T=1.022988 acc=0.124 lkl=8151.22 prior=2.037 post=8153.25
iter=1000 sigma=0.25 T=1.022988 acc=0.122 lkl=8151.22 prior=2.037 post=8153.25
iter=1020 sigma=0.25 T=1.022988 acc=0.119 lkl=8151.22 prior=2.037 post=8153.25
iter=1040 sigma=0.25 T=1.022988 acc=0.117 lkl=8151.22 prior=2.037 post=8153.25
iter=1060 sigma=0.25 T=1.022988 acc=0.115 lkl=8151.22 prior=2.037 post=8153.25
iter=1080 sigma=0.25 T=1.022988 acc=0.113 lkl=8151.22 prior=2.037 post=8153.25
iter=1100 sigma=0.25 T=1.022988 acc=0.111 lkl=8151.22 prior=2.037 post=8153.25
iter=1120 sigma=0.25 T=1.022988 acc=0.109 lkl=8151.22 prior=2.037 post=8153.25
iter=1140 sigma=0.25 T=1.022988 acc=0.107 lkl=8151.22 prior=2.037 post=8153.25
iter=1160 sigma=0.25 T=1.022988 acc=0.105 lkl=8151.22 prior=2.037 post=8153.25
iter=1180 sigma=0.25 T=1.022988 acc=0.104 lkl=8151.77 prior=2.164 post=8153.93
iter=1200 sigma=0.25 T=1.022988 acc=0.103 lkl=8154.97 prior=2.071 post=8157.04
iter=1220 sigma=0.25 T=1.022988 acc=0.102 lkl=8154.97 prior=2.071 post=8157.04
iter=1240 sigma=0.25 T=1.022988 acc=0.1 lkl=8154.97 prior=2.071 post=8157.04
iter=1260 sigma=0.25 T=1.022988 acc=0.099 lkl=8152.51 prior=2.015 post=8154.52
iter=1280 sigma=0.25 T=1.022988 acc=0.098 lkl=8152.51 prior=2.015 post=8154.52
iter=1300 sigma=0.25 T=1.022988 acc=0.097 lkl=8153.74 prior=1.916 post=8155.66
iter=1320 sigma=0.25 T=1.022988 acc=0.095 lkl=8153.74 prior=1.916 post=8155.66
iter=1340 sigma=0.25 T=1.022988 acc=0.094 lkl=8153.74 prior=1.916 post=8155.66
iter=1360 sigma=0.25 T=1.022988 acc=0.093 lkl=8153.74 prior=1.916 post=8155.66
iter=1380 sigma=0.25 T=1.022988 acc=0.091 lkl=8153.74 prior=1.916 post=8155.66
iter=1400 sigma=0.25 T=1.022988 acc=0.09 lkl=8153.74 prior=1.916 post=8155.66
iter=1420 sigma=0.25 T=1.022988 acc=0.089 lkl=8153.74 prior=1.916 post=8155.66
iter=1440 sigma=0.25 T=1.022988 acc=0.087 lkl=8153.74 prior=1.916 post=8155.66
iter=1460 sigma=0.25 T=1.022988 acc=0.086 lkl=8153.74 prior=1.916 post=8155.66
iter=1480 sigma=0.25 T=1.022988 acc=0.085 lkl=8153.74 prior=1.916 post=8155.66
iter=1500 sigma=0.25 T=1.022988 acc=0.085 lkl=8153.87 prior=1.938 post=8155.81
iter=1520 sigma=0.25 T=1.022988 acc=0.083 lkl=8153.87 prior=1.938 post=8155.81
iter=1540 sigma=0.25 T=1.022988 acc=0.082 lkl=8153.87 prior=1.938 post=8155.81
iter=1560 sigma=0.25 T=1.022988 acc=0.081 lkl=8153.87 prior=1.938 post=8155.81
iter=1580 sigma=0.25 T=1.022988 acc=0.081 lkl=8153.65 prior=1.913 post=8155.57
iter=1600 sigma=0.25 T=1.022988 acc=0.08 lkl=8153.65 prior=1.913 post=8155.57
iter=1620 sigma=0.25 T=1.022988 acc=0.079 lkl=8153.65 prior=1.913 post=8155.57
iter=1640 sigma=0.25 T=1.022988 acc=0.078 lkl=8153.65 prior=1.913 post=8155.57
iter=1660 sigma=0.25 T=1.022988 acc=0.077 lkl=8153.65 prior=1.913 post=8155.57
iter=1680 sigma=0.25 T=1.022988 acc=0.076 lkl=8153.65 prior=1.913 post=8155.57
iter=1700 sigma=0.25 T=1.022988 acc=0.075 lkl=8153.65 prior=1.913 post=8155.57
iter=1720 sigma=0.25 T=1.022988 acc=0.074 lkl=8153.65 prior=1.913 post=8155.57
iter=1740 sigma=0.25 T=1.022988 acc=0.074 lkl=8153.43 prior=1.912 post=8155.34
iter=1760 sigma=0.25 T=1.022988 acc=0.073 lkl=8153.43 prior=1.912 post=8155.34
iter=1780 sigma=0.25 T=1.022988 acc=0.072 lkl=8153.43 prior=1.912 post=8155.34
iter=1800 sigma=0.25 T=1.022988 acc=0.072 lkl=8153.43 prior=1.912 post=8155.34
iter=1820 sigma=0.25 T=1.022988 acc=0.071 lkl=8153.43 prior=1.912 post=8155.34
iter=1840 sigma=0.25 T=1.022988 acc=0.07 lkl=8153.43 prior=1.912 post=8155.34
iter=1860 sigma=0.25 T=1.022988 acc=0.069 lkl=8153.43 prior=1.912 post=8155.34
iter=1880 sigma=0.25 T=1.022988 acc=0.069 lkl=8153.43 prior=1.912 post=8155.34
iter=1900 sigma=0.25 T=1.022988 acc=0.068 lkl=8153.43 prior=1.912 post=8155.34
iter=1920 sigma=0.25 T=1.022988 acc=0.068 lkl=8146.7 prior=1.953 post=8148.65
iter=1940 sigma=0.25 T=1.022988 acc=0.067 lkl=8146.7 prior=1.953 post=8148.65
iter=1960 sigma=0.25 T=1.022988 acc=0.066 lkl=8146.7 prior=1.953 post=8148.65
iter=1980 sigma=0.25 T=1.022988 acc=0.066 lkl=8146.7 prior=1.953 post=8148.65
param	 actual	 fitted	 % error	 
A_n_B 	 3.84 	 0 	 200406.15 
A_k_B 	 0.6 	 0 	 32859.33 
tau_B 	 5.05 	 0.03 	 14778.32 
B_n_C 	 2.33 	 0.13 	 1761.73 
B_k_C 	 0.49 	 0 	 1402074.75 
tau_C 	 3.66 	 18.52 	 80.21 
> 
> initial_pars$parValues = cur_params(output=res, options=opts)
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>


Created & Maintained by Osamu Ogasawara (osamu.ogasawara@gmail.com) and IMSBIO