This is a short example of taking a stochastic system, using Fokker-Planck to convert it to a PDES and using SUNDIALS to solve a 2D partial differential equation in Haskell via the hmatrix-sundials library. The example is taken from the C examples that come with the SUNDIALS source. Here’s the full blog.
Idontgetoutmuch's Weblog
https://idontgetoutmuch.wordpress.com/ · 85 posts · history since 2008 · active
26 Apr 2020
22 Oct 2019
This is a short example of how to use SUNDIALS to solve a simple partial differential equation in Haskell via the hmatrix-sundials library. The example is taken from the C examples that come with the SUNDIALS source. Here’s the full blog. I’ll give a better URL soonish.
18 Apr 2019
Blog post on my new site.
4 Apr 2019
Blog post on my new site.
15 Jul 2018
Post available from my new site. Sadly WordPress doesn’t allow me to render the html exported by a Jupyter notebook.
19 May 2018
Introduction This blog started off life as a blog post on how I use nix but somehow transformed itself into a “how I do data visualisation” blog post. The nix is still here though quietly doing its work in the background. Suppose you want to analyze your local election results and visualize them using a … Continue reading Cartography in…
23 Apr 2018
I’m the chair this year for the first(!) ACM SIGPLAN International Workshop on Numerical Programming in Functional Languages (NPFL), which will be co-located with ICFP this September in St. Louis, Missouri, USA. Please consider submitting something! All you have to do is submit between half a page and a page describing your talk. There will … Continue reading Workshop on…
4 Mar 2018
Introduction These are some very hasty notes on Runge-Kutta methods and IRK2 in particular. I make no apologies for missing lots of details. I may try and put these in a more digestible form but not today. Some Uncomprehensive Theory In general, an implicit Runge-Kutta method is given by where and Traditionally this is written … Continue reading Implicit Runge…
25 Feb 2018
Introduction For the blog post still being written on variatonal methods, I referred to the still excellent Bishop (2006) who uses as his example data, the data available in R for the geyser in Yellowstone National Park called “Old Faithful”. While explaining this to another statistician, they started to ask about the dataset. Since I … Continue reading Reproducibility and…
2 Jun 2017
Introduction Summary Back in January, a colleague pointed out to me that GHC did not produce very efficient code for performing floating point abs. I have yet to produce a write-up of my notes about hacking on GHC: in summary it wasn’t as difficult as I had feared and the #ghc folks were extremely helpful. … Continue reading Haskell for…
18 Apr 2017
Introduction Tribbles originate from the planet Iota Geminorum IV and, according to Dr. McCoy, are born pregnant. No further details are given but we can follow Gurtin and MacCamy (1974) and perhaps recover some of what happens on the Enterprise. Of course, age-dependent population models are of more than fictional use and can be applied, for … Continue reading Trouble…
27 Jan 2017
I have been thinking about writing a blog on why the no u-turn sampler (NUTS) works rather than describing the actual algorithm. This led me to look at Jared Tobin’s Haskell implementation. His example tries to explore the Himmelblau function but only finds one local minima. This is not unexpected; as the excellent Stan manual … Continue reading Warming up…
14 Jan 2017
As part of improving the random number generation story for Haskell, I want to be able to use the testu01 library with the minimal amount of Haskell wrapping. testu01 assumes that there is a C function which returns the random number. The ghc manual gives an example but does not give all the specifics. These … Continue reading Calling Haskell…
5 Jan 2017
Introduction I was intrigued by a tweet by the UK Chancellor of the Exchequer stating "exports [to South Korea] have doubled over the last year. Now worth nearly £11bn” and a tweet by a Member of the UK Parliament stating South Korea "our second fastest growing trading partner". Although I have never paid much attention … Continue reading UK /…
24 Nov 2016
Introduction In most presentations of Riemannian geometry, e.g. O’Neill (1983) and Wikipedia, the fundamental theorem of Riemannian geometry (“the miracle of Riemannian geometry”) is given: that for any semi-Riemannian manifold there is a unique torsion-free metric connection. I assume partly because of this and partly because the major application of Riemannian geometry is General Relativity, … Continue reading Mercator: A…
4 Jul 2016
Introduction Recall from the previous post that the Hare growth parameter undergoes Brownian motion so that the further into the future we go, the less certain we are about it. In order to ensure that this parameter remains positive, let’s model the log of it to be Brownian motion. where the final equation is a … Continue reading Modelling an…
26 Jun 2016
Introduction In the 1920s, Lotka (1909) and Volterra (1926) developed a model of a very simple predator-prey ecosystem. Although simple, it turns out that the Canadian lynx and showshoe hare are well represented by such a model. Furthermore, the Hudson Bay Company kept records of how many pelts of each species were trapped for almost … Continue reading Ecology, Dynamical…
1 May 2016
Introduction This is a bit different from my usual posts (well apart from my write up of hacking at Odessa) in that it is a log of how I managed to get LibBi (Library for Bayesian Inference) to run on my MacBook and then not totally satisfactorily (as you will see if you read on). … Continue reading Fun with…
17 Apr 2016
Introduction In their paper Betancourt et al. (2014), the authors give a corollary which starts with the phrase “Because the manifold is paracompact”. It wasn’t immediately clear why the manifold was paracompact or indeed what paracompactness meant although it was clearly something like compactness which means that every cover has a finite sub-cover. It turns … Continue reading Every Manifold…
14 Feb 2016
Introduction In proposition 58 Chapter 1 in the excellent book O’Neill (1983), the author demonstrates that the Lie derivative of one vector field with respect to another is the same as the Lie bracket (of the two vector fields) although he calls the Lie bracket just bracket and does not define the Lie derivative preferring … Continue reading The Lie…
20 Jan 2016
Introduction The equation of motion for a pendulum of unit length subject to Gaussian white noise is We can discretize this via the usual Euler method where and The explanation of the precise form of the covariance matrix will be the subject of another blog post; for the purpose of exposition of forward filtering / … Continue reading Particle Smoothing
15 Jan 2016
Introduction The equation of motion for a pendulum of unit length subject to Gaussian white noise is We can discretize this via the usual Euler method where and The explanation of the precise form of the covariance matrix will be the subject of another blog post; for the purpose of exposition of using Stan and, … Continue reading Inferring Parameters…
6 Dec 2015
Introduction Let be a (hidden) Markov process. By hidden, we mean that we are not able to observe it. And let be an observable Markov process such that That is the observations are conditionally independent given the state of the hidden process. As an example let us take the one given in Särkkä (2013) where … Continue reading Naive Particle…
12 Nov 2015
Every so often, someone bitten by floating point arithmetic behaving in an unexpected way is tempted to suggest that a calculation should be done be precisely and rounding done at the end. With floating point rounding is done at every step. Here’s an example of why floating point might really be the best option for … Continue reading Floating Point:…
31 Oct 2015
Introduction We previously used importance sampling in the case where we did not have a sampler available for the distribution from which we wished to sample. There is an even more compelling case for using importance sampling. Suppose we wish to estimate the probability of a rare event. For example, suppose we wish to estimate … Continue reading Girsanov’s Theorem
12 Aug 2015
Introduction Suppose we wish to model a process described by a differential equation and initial condition But we wish to do this in the presence of noise. It’s not clear how do to this but maybe we can model the process discretely, add noise and somehow take limits. Let be a partition of then we … Continue reading Stochastic Integration
13 Jul 2015
Theorem Let and be measures on with , a sub -algebra and an integrable random variable () then Proof Thus Hence We note that is -measurable (it is the result of a projection) and that Hence as required.
20 Jun 2015
Introduction If you look at the wikipedia article on Hidden Markov Models (HMMs) then you might be forgiven for concluding that these deal only with discrete time and finite state spaces. In fact, HMMs are much more general. Furthermore, a better understanding of such models can be helped by putting them into context. Before actually … Continue reading Some Background…
7 Jun 2015
Here’s the same analysis of estimating population growth using Stan. data { int<lower=0> N; // number of observations vector[N] y; // observed population } parameters { real r; } model { real k; real p0; real deltaT; real sigma; real mu0; real sigma0; vector[N] p; k <- 1.0; p0 <- 0.1; deltaT <- 0.0005; sigma … Continue reading Population Growth…
6 Jun 2015
Introduction Let us see if we can estimate the parameter for population growth using MCMC in the example in which we used Kalman filtering. We recall the model. And we are allowed to sample at regular intervals In other words , where is known so the likelihood is Let us assume a prior of then … Continue reading Population Growth…
31 May 2015
Thames Flux It is roughly 150 miles from the source of the Thames to Kingston Bridge. If we assume that it flows at about 2 miles per hour then the water at Thames Head will have reached Kingston very roughly at days. The Environmental Agency measure the flux at Kingston Bridge on a twice daily … Continue reading The Flow…
27 Apr 2015
Introduction Suppose you want to sample from the truncated normal distribution. One way to do this is to use rejection sampling. But if you do this naïvely then you will run into performance problems. The excellent Devroye (1986) who references Marsaglia (1964) gives an efficient rejection sampling scheme using the Rayleigh distribution. The random-fu package … Continue reading Rejection Sampling
11 Mar 2015
Introduction Simple models for e.g. financial option pricing assume that the volatility of an index or a stock is constant, see here for example. However, simple observation of time series show that this is not the case; if it were then the log returns would be white noise One approach which addresses this, GARCH (Generalised AutoRegressive … Continue reading Stochastic…
25 Jan 2015
Conor McBride was not joking when he and his co-author entitled their paper about dependent typing in Haskell “Hasochism”: Lindley and McBride (2013). In trying to resurrect the Haskell package yarr, it seemed that a dependently typed reverse function needed to be written. Writing such a function turns out to be far from straightforward. How … Continue reading A Type…
25 Dec 2014
Let an be stopping times and let the filtration on which they are defined be right continuous. Then , , and are stopping times where . For the first we have and both the latter are in by the definition of a stopping time. Similarly for the second . For the fourth we have since … Continue reading Stopping Times
9 Sept 2014
Summary An extended Kalman filter in Haskell using type level literals and automatic differentiation to provide some guarantees of correctness. Population Growth Suppose we wish to model population growth of bees via the logistic equation We assume the growth rate is unknown and drawn from a normal distribution but the carrying capacity is known and … Continue reading Fun with…
26 Aug 2014
Introduction Suppose we have a vector of weights which sum to 1.0 and we wish to sample n samples randomly according to these weights. There is a well known trick in Matlab / Octave using sampling from a uniform distribution. num_particles = 2*10^7 likelihood = zeros(num_particles,1); likelihood(:,1) = 1/num_particles; [_,index] = histc(rand(num_particles,1),[0;cumsum(likelihood/sum(likelihood))]); s = sum(index); … Continue reading Haskell Vectors…
23 Aug 2014
Importance Sampling Suppose we have an random variable with pdf and we wish to find its second moment numerically. However, the random-fu package does not support sampling from such as distribution. We notice that So we can sample from and evaluate > {-# OPTIONS_GHC -Wall #-} > {-# OPTIONS_GHC -fno-warn-name-shadowing #-} > {-# OPTIONS_GHC -fno-warn-type-defaults … Continue reading Importance Sampling
6 Aug 2014
Introduction Suppose we have particle moving in at constant velocity in 1 dimension, where the velocity is sampled from a distribution. We can observe the position of the particle at fixed intervals and we wish to estimate its initial velocity. For generality, let us assume that the positions and the velocities can be perturbed at … Continue reading Fun with…
19 Jul 2014
Suppose we wish to estimate the mean of a sample drawn from a normal distribution. In the Bayesian approach, we know the prior distribution for the mean (it could be a non-informative prior) and then we update this with our observations to create the posterior, the latter giving us improved information about the distribution of … Continue reading Fun with…
15 Jun 2014
This is really intended as a draft chapter for our book. Given the diverse natures of the intended intended audiences, it is probably a bit light on explanation of the Haskell (use of monad transformers) for those with a background in numerical methods. It is hoped that the explanation of the mathematics is adequate for … Continue reading Gibbs Sampling…
24 May 2014
I have never felt entirely comfortable with Haskell’s arrows and skimming the literature for their categorical basis didn’t reveal anything as straightforward as monads or applicatives. It did however lead me to start thinking about monoidal categories and since I always want an example, I thought I would write up Hilbert spaces. Let and be … Continue reading A Monoidal…
11 May 2014
Introduction I have seen Hölder’s inequality and Minkowski’s inequality proved in several ways but this seems the most perspicuous (to me at any rate). Young’s Inequality If and such that then A and satisfying the premise are known as conjugate indices. Proof Since is convex we have Substituting in appropriate values gives or Now take … Continue reading Hölder’s and…
9 Apr 2014
Introduction It’s possible to Gibbs sampling in most languages and since I am doing some work in R and some work in Haskell, I thought I’d present a simple example in both languages: estimating the mean from a normal distribution with unknown mean and variance. Although one can do Gibbs sampling directly in R, it … Continue reading Gibbs Sampling…
2 Apr 2014
Introduction The other speaker at the Machine Learning Meetup at which I gave my talk on automatic differentiation gave a very interesting talk on A/B testing. Apparently this is big business these days as attested by the fact I got 3 ads above the wikipedia entry when I googled for it. It seems that people … Continue reading Student’s T…
20 Mar 2014
Introduction This is meant to be shorter blog post than normal with the expectation that the material will be developed further in future blog posts. A Bayesian will have a prior view of the distribution of some data and then based on data, update that view. Mostly the updated distribution, the posterior, will not be … Continue reading Bayesian Analysis:…
10 Feb 2014
Introduction Suppose we have a square thin plate of metal and we hold each of edges at a temperature which may vary along the edge but is fixed for all time. After some period depending on the conductivity of the metal, the temperature at every point on the plate will have stabilised. What is the … Continue reading Laplace’s Equation…
10 Jan 2014
I have recently started providing consultancy to a hedge fund and as far as I can see, R looks like it has a good set of libraries for this domain. In my previous job I used an embedded domain specific language in Haskell (Frankau et al. 2009). I’d like to be able to use Haskell … Continue reading Getting Financial…
2 Jan 2014
While on holiday, I started reading a couple of on-line texts on algebraic topology on my Kindle (Hatcher 2002), (May 1999). Kindle’s aren’t great for reading maths but they save having to carry around heavy books of which one might only read a fraction. Along with the usual holiday activities (ok not that usual as … Continue reading Fundamental Group…
7 Dec 2013
Introduction About a year ago there was a reddit post on the Ising Model in Haskell. The discussion seems to have fizzled out but Ising models looked like a perfect fit for Haskell using repa. In the end it turns out that they are not a good fit for repa, at least not using the … Continue reading Haskell, Ising,…
23 Oct 2013
I had a fun weekend analysing car parking data in Westminster at the Future Cities Hackathon along with Amit Nandi Bart Baddeley Jackie Steinitz Ian Ozsvald Mateusz Łapsa-Malawski Apparently in the world of car parking where Westminster leads the rest of UK follows. For example Westminster is rolling out individual parking bay monitors. Our analysis … Continue reading Parking in…
13 Oct 2013
Preface The intended audience of this article is someone who knows something about Machine Learning and Artifical Neural Networks (ANNs) in particular and who recalls that fitting an ANN required a technique called backpropagation. The goal of this post is to refresh the reader’s knowledge of ANNs and backpropagation and to show that the latter … Continue reading Backpropogation is…
12 Sept 2013
Introduction The planet Mercury has a highly elliptical orbit with a perihelion of about 0.31 AU and an aphelion of about 0.47 AU. This ellipse is not stationary but itself rotates about the Sun, a phenomenon known as the precession of the perihelion. A calculation carried out using Newtonian mechanics gives a value at variance … Continue reading The Precession…
6 Aug 2013
This article attempts to show that Haskell [@Hudak:2007:HHL:1238844.1238856] performs reasonably well on numerical problems. When I started to do this, it seemed straightforward enough: pick a problem which admitted a numerical solution, find an algorithm and code it up. I chose the problem of orbital dynamics as I had always been fascinated by the precession … Continue reading Planetary Simulation…
31 May 2013
Introduction Neural networks are a method for classifying data based on a theory of how biological systems operate. They can also be viewed as a generalization of logistic regression. A method for determining the coefficients of a given model, backpropagation, was developed in the 1970’s and rediscovered in the 1980’s. The article “A Functional Approach … Continue reading Neural Networks…
5 May 2013
This is a very informal blog of the Cabal hacking team at OdHac: Sasha Manzyuk, Bram Schuur and me (Dominic Steinitz). I really enjoyed pair programming with both of them and I certainly wouldn’t have tracked down the bug I investigated without their help. Also a big thank you to Roman Cheplyaka for organising a … Continue reading Cabal Hacking…
30 Apr 2013
Introduction Having shown how to use automated differentiation to estimate parameters in the case of linear regression let us now turn our attention to the problem of classification. For example, we might have some data about people’s social networking such as volume of twitter interactions and number of twitter followers together with a label which … Continue reading Logistic Regression…
26 Apr 2013
Introduction Automated differentiation was developed in the 1960’s but even now does not seem to be that widely used. Even experienced and knowledgeable practitioners often assume it is either a finite difference method or symbolic computation when it is neither. This article gives a very simple application of it in a machine learning / statistics … Continue reading Regression and…
23 Feb 2013
It’s part of Haskell folklore that the archetypal example for comonads is Conway’s game of life. Here’s an implementation using arrays. > {-# OPTIONS_GHC -Wall #-} > {-# OPTIONS_GHC -fno-warn-name-shadowing #-} > {-# OPTIONS_GHC -fno-warn-type-defaults #-} > import Diagrams.Prelude > import Diagrams.Backend.Cairo.CmdLine > import Data.Array > import Data.List The usual comonad class: > class Comonad … Continue reading Comonads, Life…
10 Feb 2013
Path Dependency Now we are able to price options using parallelism, let us consider a more exotic financial option. Let us suppose that we wish to price an Asian call option. The payoff at time is Thus the payoff depends not just on the value of the underlying at time but also on the path … Continue reading Parallelising Path…
5 Jan 2013
Can we get better performance for pricing financial options in Haskell? First let us translate our pricer using the Explicit Euler Method to use repa. {-# LANGUAGE FlexibleContexts, TypeOperators #-} {-# OPTIONS_GHC -Wall -fno-warn-name-shadowing -fno-warn-type-defaults #-} import Data.Array.Repa as Repa import Data.Array.Repa.Eval import Control.Monad r, sigma, k, t, xMax, deltaX, deltaT :: Double m, n, … Continue reading Option Pricing…
4 Nov 2012
One way of visualizing a distribution is to take lots of samples for it and plot the resulting histogram. To do this we use the diagrams package. {-# LANGUAGE TupleSections #-} import Diagrams.Prelude import Diagrams.Backend.Cairo.CmdLine import Data.Colour (withOpacity) To generate the actual samples (in this case from the beta distribution) we use the random-fu and … Continue reading Plotting a…
15 Aug 2012
A Simple Example In Section 7.1 of Doing Bayesian Data Analysis, John Kruschke gives a simple example of the Metropolis algorithm, in which we generate samples from a distribution without knowing the distribution itself. Of course the example is contrived as we really do know the distribtion. In the particular example, for i = 1…n where n … Continue reading…
7 May 2012
A few weeks ago, two of my colleagues and I had a birthday in the same week. What are the odds of that? There’s about 25 people in the department in which I work. We consider every outcome where pairs of people share birthdays but no 3 or more people share birthdays. Then we can … Continue reading A Triple…
22 Apr 2012
We can solve our differential equation using the Implicit Euler method which is unconditionally stable. We can also take this opportunity to use the Vector Package rather than Arrays as it has a richer set of combinators and to tidy up the code to make the payoff explicit (thanks to suggestions by Ben Moseley). First … Continue reading The Implicit…
1 Apr 2012
Suppose we want to find the price of a European call option. Then we need to solve the Black-Scholes equation: Although this particular equation can be solved explicitly, under more realistic assumptions we have to rely on numerical methods. We can approximate the partial differential equation by a difference equation (the minus sign on the … Continue reading Solving a…
17 Mar 2012
Protter(Protter 2004) notes that if some assumptions are made on how fast the mesh (see below for a definition) converges then the calculation of the almost sure value of the quadratic variation of Brownian Motion can be done using Chebyshev’s Inequality and the Borel-Cantelli Lemma. Williams and Rogers(Rogers and Williams 2000) also set this as … Continue reading The Quadratic…
28 Jan 2012
This is fairly standard. Let Xi be a symmetric random walk and let Ei be the event that Xi = 0. Then and . The diagram shows the number of ways of getting to a particular node. The dashed arrows are meant to represent that the random walk then continues. In general, at time 2n, the number … Continue reading…
8 Jan 2012
The title says it all. > {-# LANGUAGE > DeriveFunctor, > DeriveFoldable, > DeriveTraversable, > RankNTypes, > FlexibleContexts, > NoMonomorphismRestriction, > UndecidableInstances, > ScopedTypeVariables #-} > > import Prelude hiding (mapM) > import Data.Foldable > import Data.Traversable > import Data.Monoid > import Control.Monad.Writer hiding (mapM) > import Control.Monad.State hiding (mapM) > import qualified Data.Map as … Continue reading Catamorphisms are…
31 Dec 2011
One of my colleagues remarked that sets are an example of a topos; in other words that set theory can be generalized. Another of my colleagues immediately asked for another example of a topos. Here is one. A topos is a cartesian closed category with extra structure. Consider the category Sets. A subset S ⊆ X can … Continue reading…
12 Nov 2011
There don’t seem to be many examples of anamorphisms around. Here’s one to make up the lack. Let’s start off with our target language. > {-# LANGUAGE > DeriveFunctor, > DeriveFoldable, > DeriveTraversable #-} > > import Data.Foldable > import Data.Traversable > > data TermF a = PlusF a a > | MultF a a … Continue reading Anamorphism Example
17 Jul 2011
Doob’s Inequality for Right (Left) Continuous Martingales In all the proofs I have seen, some of the details are elided; here is my proof with all the details. Proof: Let be an increasing set of finite sets whose union is dense in and note that by right (left) continuity, we have: Thus: By monotone convergence … Continue reading Doob’s Inequality…
3 Jul 2011
Doob’s Inequality Let be a (discrete) non-negative sub-martingale, bounded in i.e. with for some fixed . Define by , and . Then . Further, . Proof: Fix and define the stopping time . Then which is -measurable and thus by The Stopping Time Lemma since Also by Fubini Using both of these facts and Fubini … Continue reading Doob’s Inequality
15 May 2011
One of my colleagues (Roland Zumkeller) posted some nifty functions to count the number of expressions in an AST for the DSL we work on. This led to an email and chat discussion that I have summarised in this post. Any errors are entirely mine. Let’s start off with our target language. It’s easy to … Continue reading Monadic Caching…
13 Mar 2011
is a martingale. We need to show: or alternatively that But is independent of and so
26 Jul 2009
O’Neill 2.6: from 2.13, we know that . Suppose that . Thus . And so we must have as required.
27 Dec 2008
See 3.12 of Williams. Let where and . and . is right continuous. Then by the proof . Thus and is therefore an upper bound (if not then such that and but by monotonicity ). Therefore . On the other hand suppose then (if not then but then is a lower bound for all such … Continue reading Skorokhod Representation…
29 Oct 2008
On page 19, O’Neill comments that the proof of Lemma 33 is a mild generalization of the proof of proposition 28. I think (2) (3) requires spelling out. Let and let be a co-ordinate system at . Let be a co-ordinate system at . Then by (2) has rank . Thus by exercise 7 and … Continue reading Immersions
18 May 2008
Let be a curve and be an element of the set of all smooth vector fields on . Let re-parameterize the curve. Then and . Thus which is smooth since the composition of smooth functions is a smooth function and therefore . By definition So Also by definition So, using the chain rule Which gives … Continue reading (Induced) Covariant…
12 May 2008
Nothing to do with differential geometry but there was a question on the haskell-cafe mailing list asking for “a proof that initial algebras & final coalgebras of a CPO coincide”. I presume that means a CPO-category. A category is a CPO-category iff There is a terminator, 1. Each hom-set, , is a CPO with a … Continue reading Isomorphic Types
5 May 2008
2 May 2008
Define . This is well defined since if is another bump function then and reversing and gives
26 Apr 2008
Let and be sub-manifolds of . Then and there are charts: for about adapted to for about adpapted to such that Now let be the injection map. Consider the commutative diagram below. Then which is smooth. Hence is smooth.
22 Apr 2008
Let be a vector space with linear isomorphisms to as charts. Let be a curve . Then and . But is linear (see the exercise) and so as was required to be shown.
17 Apr 2008
Let be integral curves in a manifold then is closed. Proof: Let be have a limit in . Let be a chart then converges to by continuity. But also converges to . Hence and so . Thus is closed.