Notes from 4th Bayesian Mixer Meetup
The agenda had two feature talks by Dominic Steinitz and Volodymyr Kazantsev and a lightning talk by Jon Sedar.
 |
| Dominic Steinitz: Hamiltonian and Sequential MC samplers to model ecosystems |
 |
| Volodymyr Kazantsev: Bayesian Model Averaging |
 |
| Jon Sedar: Easier Plate Notation in Python using Daft |
Next meeting
The next Bayesian Mixer Meetup meeting is already scheduled for 21 October. We will be back at Cass Business School, with two talks:- Darren Wilkinson: Hierarchical Bayesian Modelling of Growth Curves inc Stochastic Processes
- Peadar Coyle: Advanced PyMC3
Phase plane analysis in R
As an example I use here the FitzHugh-Nagumo system introduced earlier:
\begin{align}
\dot{v}=&2 (w + v - \frac{1}{3}v^3) + I_0 \\
\dot{w}=&\frac{1}{2}(1 - v - w)\\
\end{align}
The FitzHugh-Nagumo system is a simplification of the Hodgkin-Huxley model of spike generation in squid giant axon. Here \(I_0\) is a bifurcation parameter. As I decrease \(I_0\) from 0 the system dynamics change (Hopf-bifurcation): a stable equilibrium solution transform into a limit cycle.
Following Michael's paper, I can use
phaseR to plot the velocity field, add nullclines and plot trajectories from different starting points.Here I plot the FitzHugh-Nagumo system for four different parameters of \(I_0\) and three different initial starting values. The blue line show the nullcline of \(w\) i.e. \(\dot{w}=0\), while the red line shows the nullcline of \(v\). For \(I_0=-2\) I can observe the limit cycle.
Yet, I was a little surprised that the paper didn't make any references to the XPPAUT software by Bard Ermentrout, which has been around for many years as tool to analyse dynamical systems.
The GUI to the software itself gives many more options to analyse dynamical systems, including an interface to the popular bifurcation program AUTO. A good tutorial with the FitzHugh-Nagumo model was given by Mathieu Desroches at the ICS summer school 2012.
Of course I could use XPPAUT as a pure integration engine from R as well:
Considering that R started as a tool for statisticians it has made a remarkable journey; here competing with more traditional engineering tools like Matlab with MatCont or special software like XPPAUT. If someday, someone would find the time and motivation to write an interface to AUTO then R would indeed be a very good environment for the analysis of dynamical systems.
Session Info
R version 3.1.2 (2014-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
other attached packages:
[1] phaseR_1.3 deSolve_1.10-9
loaded via a namespace (and not attached):
[1] tools_3.1.2Simulating neurons or how to solve delay differential equations in R
My model is based on the paper: Epileptiform activity in a neocortical network: a mathematical model by F. Giannakopoulos, U. Bihler, C. Hauptmann and H. J. Luhmann. The article presents a flexible and efficient modelling framework for:
- large populations with arbitrary geometry
- different synaptic connections with individual dynamic characteristics
- cell specific axonal dynamics
Hodgkin-Huxley model in R
One of the great research papers of the 20th century celebrates its 60th anniversary in a few weeks time: A quantitative description of membrane current and its application to conduction and excitation in nerve by Alan Hodgkin and Andrew Huxley. Only shortly after Andrew Huxley died, 30th May 2012, aged 94.
In 1952 Hodgkin and Huxley published a series of papers, describing the basic processes underlying the nervous mechanisms of control and the communication between nerve cells, for which they received the Nobel prize in physiology and medicine, together with John Eccles in 1963.
Their research was based on electrophysiological experiments carried out in the late 1940s and early 1950 on a giant squid axon to understand how action potentials in neurons are initiated and propagated.
Dynamical systems in R with simecol
This evening I will talk about Dynamical systems in R with simecol at the LondonR meeting.
Thanks to the work by Thomas Petzoldt, Karsten Rinke, Karline Soetaert and R. Woodrow Setzer it is really straight forward to model and analyse dynamical systems in R with their deSolve and simecol packages.
I will give a brief overview of the functionality using a predator-prey model as an example.
This is of course a repeat of my presentation given at the Köln R user group meeting in March.
For a further example of a dynamical system with simecol see my post about the Hodgkin-Huxley model, which describes the action potential of a giant squid axon.
I shouldn't forget to mention the other talks tonight as well:
- Writing R for Dummies - Andrie De Vries
- News from data.table 1.6, 1.7 and 1.8 - Matthew Dowle
- Converting S Plus Applications into R - Andy Nicholls (postponed to 18 September 2012)
Logistic map: Feigenbaum diagram in R
The other day I found some old basic code I had written about 15 years ago on a Mac Classic II to plot the Feigenbaum diagram for the logistic map. I remember, it took the little computer the whole night to produce the bifurcation chart.
logistic.map <- function(r, x, N, M){
## r: bifurcation parameter
## x: initial value
## N: number of iteration
## M: number of iteration points to be returned
z <- 1:N
z[1] <- x
for(i in c(1:(N-1))){
z[i+1] <- r *z[i] * (1 - z[i])
}
## Return the last M iterations
z[c((N-M):N)]
}
## Set scanning range for bifurcation parameter r
my.r <- seq(2.5, 4, by=0.003)
system.time(Orbit <- sapply(my.r, logistic.map, x=0.1, N=1000, M=300))
## user system elapsed (on a 2.4GHz Core2Duo)
## 2.910 0.018 2.919
Orbit <- as.vector(Orbit)
r <- sort(rep(my.r, 301))
plot(Orbit ~ r, pch=".")
Let's not forget when Mitchell Feigenbaum started this work in 1975 he did this on his little calculator!
Update, 18 March 2012
The comment from Berend has helped to speedup the code by a factor of about four, thanks to byte compiling (using the same parameters as above), and Owe got me thinking about the alpha value of the plotting colour. Here is the updated result, with the R code below:
library(compiler) ## requires R >= 2.13.0
logistic.map <- cmpfun(logistic.map) # same function as above
my.r <- seq(2.5, 4, by=0.001)
N <- 2000; M <- 500; start.x <- 0.1
orbit <- sapply(my.r, logistic.map, x=start.x, N=N, M=M)
Orbit <- as.vector(orbit)
r <- sort(rep(my.r, (M+1)))
plot(Orbit ~ r, pch=".", col=rgb(0,0,0,0.05))
