News aggregator

Uploading to hackage fails: bad file names in tarball

haskell-cafe - Mon, 06/16/2014 - 9:44am
Since some time, I have been unable to upload packages to Hackage, via either the 'cabal upload' command, or using the web form. The error message it gives is this: ==== ✂ ==== Hackage username: AriePeterson Hackage password: Uploading dist/np-linear-0.1.1.1.tar.gz... Error: dist/np-linear-0.1.1.1.tar.gz: 400 Bad Request Error: Invalid package Invalid windows file name in tar archive: "np-linear-0.1.1.1\\src\\Aux.hs". For portability, hackage requires that file names be valid on both Unix and Windows systems, and not refer outside of the tarball. ==== ✂ ==== The tarball is created by 'cabal sdist': cabal-install version 1.18.0.2 using version 1.18.1.1 of the Cabal library. I also installed the newest cabal-install, on another machine, but this did not help (same error). I also tried to create a tarball by hand, using 'tar --format=ustar', but this again resulted in the same error message. By the way, I'm on linux, not Windows, so it is not clear how the backslashes get in the file names. Wh
Categories: Offsite Discussion

Using newtype to "tag" values with a context

Haskell on Reddit - Mon, 06/16/2014 - 8:22am

I'm not really sure if this pattern has a name, but it occured to me a few days ago and I've just been thinking about this, but haven't really used it.

Say that you have a function magic which accepts a CreditCard type and does some magic with it. At that point you really expect the card to be already validated (which happens up in the call stack). But since the CreditCard type is just a value, it doesn't carry around the context having passed validation.

Just to make this clear, I'm not talking about simple string validation, but something like checking with a payment gateway and making a $0 transaction, to actually check that the card is valid. Think of this example as something you only want to validate once.

Now we've established that the magic function should never be called with a raw CreditCard that hasn't been through the validation, a program that does such thing shouldn't even compile, right?

The thing I came up with is to use a newtype wrapper, something like newtype Valid a = Valid a, so that I can then have the function type defined as magic :: Valid CreditCard -> IO ().

This way when I get a new CreditCard from a user or something, I can run it through a validation which would return either a Valid CreditCard, or some error. For example validate :: CreditCard -> Either ValidationError (Valid CreditCard).

What do you guys think about this? Do you do this? Does it make sense? Does it have a name?

submitted by progfu
[link] [9 comments]
Categories: Incoming News

MonadReactive, embedded monadic reactiveness

Haskell on Reddit - Mon, 06/16/2014 - 7:55am

Hi! Today I decided to present something I’ve been wandering around for some days now.

The initial problem is this:

Given a context with values, how could we react to a change in that context?

I came up with a very small solution, but that works well for me. Here’s a lpast link to my idea – I’ll make a package later if we find it cool and sexy.

The idea: you have two concepts: the monad you want to react to (m), and the monad you want to react in (r). When something happens in m, something should happen in r. We can state that m is observed by r. It’s then quite simple: my library provides a way to wrap monadic functions in the observed monad into the observer so that when such a function is called, the observer reacts.

I know FRP, and I know it’s great, but I wanted to present something way simpler, for not-that simpler problems.

What do you think folks? Should I make a monad-reactive package, or do you think it’s worthless?

submitted by _skp
[link] [13 comments]
Categories: Incoming News

BNFC-meta being maintained? Is there public src repo?

haskell-cafe - Mon, 06/16/2014 - 5:33am
Dear BNFC-meta developers and users, I've been using BNFC-meta for the parser of my experimental language Nax language implementation ( https://github.com/kyagrd/mininax ). It's very cool 'cause you don't need to write makefile or cabal entries to call the command line tools (bnfc, alex, happy), and you get quasiquoters for free, which really makes the early development pleasant. But due to the bug fixed in recent versions of BNFC, probably after latest BNFC-meta release, I'll have to switch back to using plain BNFC at some point. (FYI,the bug is that block comments arn't working). If there is a source repository that has ported more recent versions of BNFC into BNFC-meta, it'll be helpful for people like me. Or, if the source repository has open access and easy to participate (e.g., putting on github or something like that), we can draw community support maintaining BNFC-meta to keep up with the recent BNFC (and also happy and alex) versions. -- Ki Yung Ahn
Categories: Offsite Discussion

type families and type classes

haskell-cafe - Sun, 06/15/2014 - 11:25pm
Hi all, I have a bunch of data types which are "divided" in other data types as following: type FileName = String data FileWish = File FileName deriving (Show, Eq) data FileFact = FileExist | FileDontExist deriving (Show, Eq) data FileAction = FileCreate deriving (Show, Eq) observe :: FileWish -> [FileFact] observe :: undefined plan :: FileWish -> FileFact -> [FileAction] plan (File _) FileDontExist = [FileCreate] plan (File _) FileExist = [] perform :: FileWish -> [FileAction] -> IO () perform :: undefined I want to have a model as extensible as possible, even if I'm new to Haskell, I think that it's preferable to base parts of system on behavior (type classes) rather than data types. So I have File that give a Wish, a Fact (via observe) and a Action (via plan). I want to put observe, plan and perform in different type classes (let's say Observable, Plannable and Performable). When I create a new thing (like File) I want, fo
Categories: Offsite Discussion

What are the problems with instances for polymorphictypes?

haskell-cafe - Sun, 06/15/2014 - 8:47pm
In other words instances for forall-types, such as: instance Foo (forall a. [a]) where ... It feels obvious to me that there *would* be problems with this, but I'm curious about what, exactly, they are. Could someone familiar with the matter either elaborate on them, or refer me to an existing explanation, a previous discussion, or something of the sort? I *don't* have any kind of use case in mind, I'm merely seeking a better understanding of the type-system issues involved. (I attempted Google, but didn't have much success.) Thanks in advance. _______________________________________________ Haskell-Cafe mailing list Haskell-Cafe< at >haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Categories: Offsite Discussion

Can't get satisfaction from tying knots

Haskell on Reddit - Sun, 06/15/2014 - 6:20pm

I'm currently working on language-gcl, a library that currently doesn't do much more than implement Dijkstra's guarded command language (GCL). GCL only defines two 'building blocks' of a programming language, because the manuscript it's defined in (EWD 418) focuses mostly on formal methods.

The building blocks of GCL are the alternative construct (if ... fi) and the repetitive construct (do ... od), which contain guarded commands. Dijkstra doesn't define the expression language (except that the guards must be decidable predicates) and explicitly leaves room for 'other statements'; he uses examples of (concurrent) assignment, e.g.:

if true -> x := Y; y := Y; do x > y -> x := x - y [] y > x -> x, y := y, x od fi

So anything if...fi or do...od is defined by GCL, but everything else (the guard expressions, the assignment statements) is up to whoever uses GCL.

So I decided to start with this:

data GCL stmt guard = Alternative [GC stmt guard] | Repetitive [GC stmt guard] | Statement stmt data GC stmt guard = GC guard [GCL stmt guard]

This allows for some extensibility, such as concurrent assignment, assert and assume statements, and many other things.

I ran into problems with the extensibility on two distinct occasions. One involved adding loop invariants (type parameter explosion!), but I believe that my current solution (adding a [guard] field to Repetitive) is reasonable.

The other is a little more problematic, which (initially) involves declaring variables. If I want to extend GCL with any constructs that may contain more GCL, I run into knot-tying problems. Ideally, I want to simply define something like this:

data Variables name ty expr = Declare [(name, ty)] gcl | Assign [(name, expr)]

-edit-

I could end this right here, with (gcl ~ GCL (MoreStmts name ty expr) expr), but this doesn't really sit right with me, DRY being one problem:

data Exceptions gcl = Raise | Rescue gcl gcl

It would also be really nice to be able to spice things up à la carte, e.g. GCL (Exceptions :+: Variables name ty expr) expr, but I would lose that by using this solution.

So, to tie the knot, I could make gcl the last type variable of MoreStmts and introduce the horrific nice auxiliary type

newtype GclKnot stmt1 guard = GclKnot (GCL (stmt1 (GclKnot stmt1 guard)) guard)

but I fear the doubling of type variables will get me into all sorts of trouble once I want to use these things with type classes. I can also forget about deriving any instances.

Does anyone have an idea for improvement? Or am I trying too hard?

submitted by Rhymoid
[link] [6 comments]
Categories: Incoming News

I'm trying my hardest to learn this language and only get more and more discouraged.

Haskell on Reddit - Sun, 06/15/2014 - 4:28pm

Edit: I am blown away by how many really helpful and encouraging responses I've received here. I really wasn't expecting that. Thanks to everyone. Some thing's I've learned here: I need to read more about typeclasses and spend significantly more time studying the language than I have in the past. I should definitely check out Real World Haskell. I should get on #haskell on freenode as the people there are very knowledgable. Thanks again everyone. Not only do I now have a proper understanding of what was wrong with my very simple code I also have a new plan for learning. Thank you.

Every 6 months or so I read through LYAH (and feel utterly confused by the end) and try to do something super super super basic, I mean so trivial it's embarrassing, and I get absolutely nowhere. I spend most of my time writing in Ruby and Javascript and have recently gotten pretty comfortable with Clojure but there is something about Haskell that just wont stick in my brain.

The latest bit of code to drive me crazy is this:

mean :: [a] -> Float mean x = sum x / length x

This seems pretty straight forward to me. Sum up the list of any type and then divide by the size of the list. I chose a because it could be a list of integers or floats (probably incorrect). This is what I get when I attempt to load the file in ghci...

Couldn't match expected type `Float' with actual type `Int' In the return type of a call of `length' In the second argument of `(/)', namely `length x' In the expression: sum x / length x

If I leave out the type declaration I get a different error but an error nonetheless. Does anyone have any advice as to how to learn this language? I consider myself to be of slightly above average intelligence and I'm fairly determined to figure this out. I don't want to just give up.

submitted by dunnowins
[link] [110 comments]
Categories: Incoming News

Taking a number to a power and taking its mod

haskell-cafe - Sun, 06/15/2014 - 4:10pm
Hello, I've been trying to do a problem from SPOJ which requires me to find prime numbers in a range. I've been battling with several possibility, so far I have had no luck. I remember primality tests from college. So I looked up the one I knew: fermat's little theorem. It requires me to do this check: a^(p-1) `mod` p == 1 However, taking a number to a large power is very slow. So it was even worse than the solution I was already trying. I'm not here for your guys to give me the solution. But I'm here to understand how can I do that operation quickly. This wiki entry implements powMod http://www.haskell.org/haskellwiki/Testing_primality Calling that funcion with powMod 999999527 2 (999999527-1) is very fast. What's going on? Can someone shed me some light? []'s Rafael _______________________________________________ Haskell-Cafe mailing list Haskell-Cafe< at >haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Categories: Offsite Discussion

Noob here: How can I define a list length function non-recursively?

Haskell on Reddit - Sun, 06/15/2014 - 12:47pm

It's a thought-exercise from 'the Haskell School of Music.' I was staring at it for a long time to no avail. I realize the answer's probably obvious.

submitted by hoomei
[link] [19 comments]
Categories: Incoming News

Support of multi-constructor types by UNPACK pragma

haskell-cafe - Sun, 06/15/2014 - 12:43pm
Currently the pragma only supports single-constructor types. So in the following example it will simply be ignored: data A = A1 Char | A2 {-# UNPACK #-} !B data B = B1 Int | B2 Bool However it seems to be easily solvable by changing the type A to something like the following during unpacking: data A = A1 Char | A2_1 Int | -- from B1 A2_2 Bool -- from B2 Am I missing something? Why is it not implemented? ​ _______________________________________________ Haskell-Cafe mailing list Haskell-Cafe< at >haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Categories: Offsite Discussion

Yesod Web Framework: Three evil ways to avoid conditional compilation

Planet Haskell - Sun, 06/15/2014 - 8:15am

Let's suppose that you're using a library at version 1. It exposes a function:

someFunc :: Int -> String someFunc i = 'x' : show i

In version 2 of the library, someone realizes that this library isn't general-purpose enough: why is the 'x' character hardcoded? So version 2 exposes a more powerful version of the function:

someFunc :: Int -> Char -> String someFunc i c = c : show i

In your current codebase, you have:

main = putStrLn $ someFunc 5

Changing that code to work with version 2 is trivial:

main = putStrLn $ someFunc 5 'x'

But what if you want to make your code work with both versions? The real, proper answer, that I hope everyone actually uses, is to use Cabal CPP macros:

#if MIN_VERSION_somelib(2, 0, 0) main = putStrLn $ someFunc 5 'x' #else main = putStrLn $ someFunc 5 #endif

And sure, you should do that... but let's have some fun. I'm going to present three evil techniques to accomplish the same conditional compilation result, and try to point out their relative merits. I encourage others to come up with other ridiculous ideas of their own.

If anyone's curious where I came up with the idea to do this, it was thinking about prerelease GHC patch level releases that break backwards compatibility. And if you want to play around with the code, either open it in FP Haskell Center or clone the Github repo. In all the examples, simply switch whether V1 or V2 is imported to simulated upgrading/downgrading dependencies.

Typeclasses

A well known truth in Haskell is, "for every problem, there's an abuse of typeclasses waiting to be discovered." As far as typeclass abuses go, this one is pretty benign.

{-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE TypeSynonymInstances #-} module Typeclass where import V1 -- import V2 class Evil a where evil :: a -> String instance Evil String where evil = id instance Evil a => Evil (Char -> a) where evil f = evil (f 'x') main :: IO () main = putStrLn $ evil $ someFunc 5

What's "nice" about this approach (if anything here is nice) is that everything is compile-time checked, and the code is actually pretty readable. However, as the different cases you want to support get more complicated, you'll need to add in ever harrier language extensions.

Typeable

This approach is for the Pythonista in you. Next time someone proudly states that Haskell is a statically typed language, just pull this one out:

module Typeable where import Data.Typeable -- import V1 import V2 evil :: Typeable a => a -> String evil a | Just s <- cast a = s | Just f <- cast a = f 'x' | otherwise = error "Yay, runtime type errors!" main :: IO () main = putStrLn $ evil $ someFunc 5

Advantage: it's so incredibly trivial to add more cases. Downsides: it's runtime type checking, and the dispatch is performed at runtime, not compile time.

Template Haskell

Of course, no blog post about Haskell abuse would be complete without some usage of Template Haskell. Due to the stage restriction, we have to split this into two files. First, THHelper:

{-# LANGUAGE TemplateHaskell #-} module THHelper where import Language.Haskell.TH evil :: Name -> Q Exp evil name = do VarI _ typ _ _ <- reify name case typ of AppT (AppT ArrowT (ConT _)) (ConT _) -> return $ VarE name AppT (AppT ArrowT (ConT _)) (AppT (AppT ArrowT (ConT _)) (ConT _)) -> [|flip $(return $ VarE name) 'x'|] _ -> error $ "Unknown type: " ++ show typ

Notice how beautiful our pattern matching is. This combines the best (worst) of both worlds from above: we get full compile time checking, and can easily (hah!) pattern match on all possible signatures for the function at hand.

And calling this beast is equally elegant:

{-# LANGUAGE TemplateHaskell #-} module TH where import THHelper import V1 -- import V2 main :: IO () main = putStrLn $ $(evil 'someFunc) 5
Categories: Offsite Blogs

Open source projects using postgresql-simple?

Haskell on Reddit - Sun, 06/15/2014 - 7:54am

A while ago, I made a post asking for advice on how to use databases in Haskell. I've played around a bit, and decided that postgresql-simple might be the best way to go. However, I was wondering if there were any examples floating around of projects that have used it in production? Googling didn't turn anything up.

Writing the interface between Haskell and database types seems like it is the sort of thing that you would have to be quite disciplined in doing. If anyone could point me at some code that shows examples of doing so, I would be very grateful!

submitted by parrotchute
[link] [9 comments]
Categories: Incoming News

Dominic Steinitz: Gibbs Sampling in Haskell

Planet Haskell - Sun, 06/15/2014 - 6:13am

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 those with a background in Haskell but not necessarily in numerical methods. As always, any feedback is gratefully accepted.

Introduction

Imagine an insect, a grasshopper, trapped on the face of a clock which wants to visit each hour an equal number of times. However, there is a snag: it can only see the value of the hour it is on and the value of the hours immediately anti-clockwise and immediately clockwise. For example, if it is standing on 5 then it can see the 5, the 4, and the 6 but no others.

It can adopt the following strategy: toss a fair coin and move anti-clockwise for a head and move clockwise for a tail. Intuition tells us that over a large set of moves the grasshopper will visit each hour (approximately) the same number of times.

Can we confirm our intuition somehow? Suppose that the strategy has worked and the grasshopper is now to be found with equal probability on any hour. Then at the last jump, the grasshopper must either have been at the hour before the one it is now on or it must have been at the hour after the one it is now on. Let us denote the probability that the grasshopper is on hour by and the (conditional) probability that the grasshopper jumps to state given it was in state by . Then we have

Substituting in where is a normalising constant (12 in this case) we obtain

This tells us that the required distribution is a fixed point of the grasshopper’s strategy. But does the strategy actually converge to the fixed point? Let us perform an experiment.

First we import some modules from hmatrix.

> {-# LANGUAGE FlexibleContexts #-} > module Chapter1 where > import Data.Packed.Matrix > import Numeric.LinearAlgebra.Algorithms > import Numeric.Container > import Data.Random > import Control.Monad.State > import qualified Control.Monad.Writer as W > import qualified Control.Monad.Loops as ML > import Data.Random.Source.PureMT

Let us use a clock with 5 hours to make the matrices sufficiently small to fit on one page.

Here is the strategy encoded as a matrix. For example the first row says jump to position 1 with probablity 0.5 or jump to position 5 with probability 0.5.

> eqProbsMat :: Matrix Double > eqProbsMat = (5 >< 5) > [ 0.0, 0.5, 0.0, 0.0, 0.5 > , 0.5, 0.0, 0.5, 0.0, 0.0 > , 0.0, 0.5, 0.0, 0.5, 0.0 > , 0.0, 0.0, 0.5, 0.0, 0.5 > , 0.5, 0.0, 0.0, 0.5, 0.0 > ]

We suppose the grasshopper starts at 1 o’clock.

> startOnOne :: Matrix Double > startOnOne = ((1 >< 5) [1.0, 0.0, 0.0, 0.0, 0.0])

If we allow the grasshopper to hop 1000 times then we see that it is equally likely to be found on any hour hand with a 20% probability.

ghci> eqProbsMat (5><5) [ 0.0, 0.5, 0.0, 0.0, 0.5 , 0.5, 0.0, 0.5, 0.0, 0.0 , 0.0, 0.5, 0.0, 0.5, 0.0 , 0.0, 0.0, 0.5, 0.0, 0.5 , 0.5, 0.0, 0.0, 0.5, 0.0 ] ghci> take 1 $ drop 1000 $ iterate (<> eqProbsMat) startOnOne [(1><5) [ 0.20000000000000007, 0.2, 0.20000000000000004, 0.20000000000000004, 0.2 ]]

In this particular case, the strategy does indeed converge.

Now suppose the grasshopper wants to visit each hour in proportion the value of the number on the hour. Lacking pen and paper (and indeed opposable thumbs), it decides to adopt the following strategy: toss a fair coin as in the previous strategy but only move if the number is larger than the one it is standing on; if, on the other hand, the number is smaller then choose a number at random from between 0 and 1 and move if this value is smaller than the ratio of the proposed hour and the hour on which it is standing otherwise stay put. For example, if the grasshopper is standing on 5 and gets a tail then it will move to 6 but if it gets a head then four fifths of the time it will move to 4 but one fifth of the time it will stay where it is.

Suppose that the strategy has worked (it is not clear that is has) and the grasshopper is now to be found at 12 o’clock 12 times as often as at 1 o’clock, at 11 o’clock 11 times as often as at 1 o’clock, etc. Then at the last jump, the grasshopper must either have been at the hour before the one it is now on, the hour after the one it is now on or the same hour it is now on. Let us denote the probability that the grasshopper is on hour by .

Substituting in at 4 say

The reader can check that this relationship holds for all other hours. This tells us that the required distribution is a fixed point of the grasshopper’s strategy. But does this strategy actually converge to the fixed point?

Again, let us use a clock with 5 hours to make the matrices sufficiently small to fit on one page.

Here is the strategy encoded as a matrix. For example the first row says jump to position 1 with probablity 0.5 or jump to position 5 with probability 0.5.

> incProbsMat :: Matrix Double > incProbsMat = scale 0.5 $ > (5 >< 5) > [ 0.0, 1.0, 0.0, 0.0, 1.0 > , 1.0/2.0, 1.0/2.0, 1.0, 0.0, 0.0 > , 0.0, 2.0/3.0, 1.0/3.0, 1.0, 0.0 > , 0.0, 0.0, 3.0/4.0, 1.0/4.0, 1.0 > , 1.0/5.0, 0.0, 0.0, 4.0/5.0, 1.0/5.0 + 4.0/5.0 > ]

We suppose the grasshopper starts at 1 o’clock.

If we allow the grasshopper to hop 1000 times then we see that it is equally likely to be found on any hour hand with a probability of times the probability of being found on 1.

ghci> incProbsMat (5><5) [ 0.0, 0.5, 0.0, 0.0, 0.5 , 0.25, 0.25, 0.5, 0.0, 0.0 , 0.0, 0.3333333333333333, 0.16666666666666666, 0.5, 0.0 , 0.0, 0.0, 0.375, 0.125, 0.5 , 0.1, 0.0, 0.0, 0.4, 0.5 ] ghci> take 1 $ drop 1000 $ iterate (<> incProbsMat) startOnOne [(1><5) [ 6.666666666666665e-2, 0.1333333333333333, 0.19999999999999996, 0.2666666666666666, 0.33333333333333326 ]]

In this particular case, the strategy does indeed converge.

Surprisingly, this strategy produces the desired result and is known as the Metropolis Algorithm. What the grasshopper has done is to construct a (discrete) Markov Process which has a limiting distribution (the stationary distribution) with the desired feature: sampling from this process will result in each hour being sampled in proportion to its value.

Markov Chain Theory

Let us examine what is happening in a bit more detail.

The grasshopper has started with a very simple Markov Chain: one which jumps clockwise or anti-clockwise with equal probability and then modified it. But what is a Markov Chain?

A time homogeneous Markov chain is a countable sequence of random variables
such that

We sometimes say that a Markov Chain is discrete time stochastic process with the above property.

So the very simple Markov Chain can be described by

The grasshopper knows that so it can calculate without knowing . This is important because now, without knowing , the grasshopper can evaluate

where takes the maximum of its arguments. Simplifying the above by substituing in the grasshopper’s probabilities and noting that is somewhat obscure way of saying jump clockwise or anti-clockwise we obtain

The Ergodic Theorem

In most studies of Markov chains, one is interested in whether a chain has a stationary distribution. What we wish to do is take a distribution and create a chain with this distribution as its stationary distribution. We will still need to show that our chain does indeed have the correct stationary distribution and we state the relevant theorem somewhat informally and with no proof.

Theorem

An irreducible, aperiodic and positive recurrent Markov chain has a unique stationary distribution.

Roughly speaking

  • Irreducible means it is possible to get from any state to any other state.

  • Aperiodic means that returning to a state having started at that state occurs at irregular times.

  • Positive recurrent means that the first time to hit a state is finite (for every state and more pedantically except on sets of null measure).

Note that the last condition is required when the state space is infinite – see Skrikant‘s lecture notes for an example and also for a more formal definition of the theorem and its proof.

Algorithm

Let be a probability distribution on the state space with for all and let be an ergodic Markov chain on with transition probabilities (the latter condition is slightly stronger than it need be but we will not need fully general conditions).

Create a new (ergodic) Markov chain with transition probabilities

where takes the maximum of its arguments.

Calculate the value of interest on the state space e.g. the total magnetization for each step produced by this new chain.

Repeat a sufficiently large number of times and take the average. This gives the estimate of the value of interest.

Convergence

Let us first note that the Markov chain produced by this algorithm almost trivially satisfies the detailed balance condition, for example,

Secondly since we have specified that is ergodic then clearly is also ergodic (all the transition probabilities are ).

So we know the algorithm will converge to the unique distribution we specified to provide estimates of values of interest.

Gibbs Sampling Random Scan

For simplicity let us consider a model with two parameters and that we sample from either parameter with equal probability. In this sampler, We update the parameters in a single step.

The transition density kernel is then given by

where is the Dirac delta function.

Detailed balance

This sampling scheme satisifies the detailed balance condition. We have

In other words

Hand waving slightly, we can see that this scheme satisfies the premises of the ergodic theorem and so we can conclude that there is a unique stationary distribution and must be that distribution.

Systematic Scan

Most references on Gibbs sampling do not describe the random scan but instead something called a systematic scan.

Again for simplicity let us consider a model with two parameters. In this sampler, we update the parameters in two steps.

We observe that this is not time-homegeneous; at each step the transition matrix flips between the two transition matrices given by the individual steps. Thus although, as we show below, each individual transtion satisifies the detailed balance condition, we cannot apply the ergodic theorem as it only applies to time-homogeneous processes.

The transition density kernel is then given by

where .

Thus

Detailed balance

Suppose that we have two states and and that . Then . Trivially we have

Now suppose that

So again we have

Similarly we can show

But note that

whereas

and these are not necessarily equal.

So the detailed balance equation is not satisfied, another sign that we cannot appeal to the ergodic theorem.

An Example: The Bivariate Normal

Let us demonstrate the Gibbs sampler with a distribution which we actually know: the bivariate normal.

The conditional distributions are easily calculated to be

Let’s take a correlation of 0.8, a data point of (0.0, 0.0) and start the chain at (2.5, 2.5).

> rho :: Double > rho = 0.8 > > y :: (Double, Double) > y = (0.0, 0.0) > > y1, y2 :: Double > y1 = fst y > y2 = snd y > > initTheta :: (Double, Double) > initTheta = (2.5, 2.5)

We pre-calculate the variance needed for the sampler.

> var :: Double > var = 1.0 - rho^2

In Haskell and in the random-fu package, sampling from probability distributions is implemented as a monad. We sample from the relevant normal distributions and keep the trajectory using a writer monad.

> gibbsSampler :: Double -> RVarT (W.Writer [(Double,Double)]) Double > gibbsSampler oldTheta2 = do > newTheta1 <- rvarT (Normal (y1 + rho * (oldTheta2 - y2)) var) > lift $ W.tell [(newTheta1, oldTheta2)] > newTheta2 <- rvarT (Normal (y2 + rho * (newTheta1 - y1)) var) > lift $ W.tell [(newTheta1, newTheta2)] > return $ newTheta2

It is common to allow the chain to “burn in” so as to “forget” its starting position. We arbitrarily burn in for 10,000 steps.

> burnIn :: Int > burnIn = 10000

We sample repeatedly from the sampler using the monadic form of iterate. Running the monadic stack is slightly noisy but nonetheless straightforward. We use mersenne-random-pure64 (albeit indirectly via random-source) as our source of entropy.

> runMCMC :: Int -> [(Double, Double)] > runMCMC n = > take n $ > drop burnIn $ > snd $ > W.runWriter (evalStateT (sample (ML.iterateM_ gibbsSampler (snd initTheta))) (pureMT 2))

We can look at the trajectory of our sampler for various run lengths.

For bigger sample sizes, plotting the distribution sampled re-assures us that we are indeed sampling from a bivariate normal distribution as the theory predicted.

Applications to Bayesian Statistics

Some of what is here and here excluding JAGS and STAN (after all this is a book about Haskell).

Applications to Physics

Applications to Physics

Most of what is here.


Categories: Offsite Blogs