Probabilistic programming in Haskell and Python (Langevin again).

[Dilawar Singh]

In Haskell, any computation which creates side effect (modification of state) is done in a Monad. Monad (and Arrows et. al.) are the way to do "stateful computation" in Haskell. Writing to a file, printing string to console, generating random number are all stateful computations. Understanding and using Monad is perhaps the hardest part of learning Haskell. This note is not a tutorial about Monads, but a working code snippet which solves a similar problem of Langevin simulation posted few days ago.

To do a Langevin simulation, one needs to generate a Poisson distribution with mean 0 and variance 1. I've used Haskell package `normaldistribuition` to do this. Let's write down the code.

-- This function creates an infinite list of normally distributed numbers with mean 0 and variance 1.
-- To take n numbers out of this infinite list, one uses `take n alphas`. 
alphas :: (Random a, Floating a) => IO [a]
alphas = do
    g <- newStdGen
    return $ normals g

-- A normal distribution. Draw a number from this distribution by repeated calling this function.
-- alpha :: (Random a, Floating a) => IO a
alpha :: IO Double
alpha = normalIO

Rest is now straight-forward, multiply alpha with Random Walk term (Wiener process) and add to the differential equation. Now simulate it using step by step. The full code with commentary should be available here. For the curious souls, this simulates a Genetic Switch modified in this paper  and first explored by this guy. Paper should be worth going through for "bistability in system dynamics" people.

The idea behind this implementation was to get a comparison b/w Python and Haskell for numerical simulations.

When it comes to stateful computations, nothing is straightforward in Haskell. The type system of Haskell is just amazing and lives up to its expectation; often frustrating. Speed is bit faster than Python in interpreter; 3-5 times faster with compiled code. It can easily be improved further by using Strict evaluation but I did not bother much about it. Code in Haskell more readable, and much easier to parallelize. Plotting is a headache in Haskell.

What a post about simulation without a figure. First in in Haskell (used Chart-cairo library) and second in Python (matplotlib).



​A note about performance b/w Data.List and Data.Sequence is added to this stackoverflow question.
--

Dilawar
NCBS Bangalore

[Saket Choudhary]

Hey Dilawar,

Cool stuff!

I think you meant 'Normal' with mean 0, var=1 instead of a
'poisson'?(mean and var are otherwise contradictory for a poisson)

(https://github.com/dilawar/Courses/blob/master/RandomnessInBiology_Monsoon_2014/Assignment07/solve.py#L46)

> A note about performance b/w Data.List and Data.Sequence is added to this stackoverflow question.
> --
> Dilawar
> NCBS Bangalore
>

> --
> --
> The website for the club is http://wncc-iitb.org/
> To post to this group, send email to wncc...@googlegroups.com
> ---
> You received this message because you are subscribed to the Google Groups "Web and Coding Club IIT Bombay" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to wncc_iitb+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Dilawar Singh]

Oh yeah; good catch about mean == var for Poisson. Thanks.

best,
Dilawar