Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: MCMC in Mathematica

  • To: mathgroup at smc.vnet.net
  • Subject: [mg96012] Re: MCMC in Mathematica
  • From: Mark Fisher <particlefilter at gmail.com>
  • Date: Sun, 1 Feb 2009 04:41:06 -0500 (EST)
  • References: <gm0q4e$rmq$1@smc.vnet.net>

On Jan 31, 1:12 am, Ramiro <ramiro.barran... at gmail.com> wrote:
> Hello,
>
> I use Mathematica a lot and now I would like to do MCMC to work on
> some of my problems.  Is this something that people with experience
> would recommend?  Or should I move to a different programming language
> such as C++ or R?  I have done a little samples of MCMC on mathematica
> and it does not take much code to get it running, but I am wondering
> if I would be better off in the long-run to work on something that
> would be faster.
>
> Thank you,
> Ramiro

Hi Ramiro,

I wrote a package to do random-walk Metropolis draws and I use it to
do pretty much all of my MCMC. (When I want to do Gibbs sampling, I
write specific code for each problem.) You can get it here:

http://www.markfisher.net/~mefisher/mma/mathematica.html

Here's a simple illustration of how to use it. (Please read the
instructions at the web site on where to put the file.) BTW, I'm using
Version 7, although the code works in Version 6 too.

<< RandomWalkMetropolis`

data = RandomReal[NormalDistribution[], 30];

expr = Product[
     PDF[NormalDistribution[\[Mu], \[Sigma]], data[[i]]], {i, 30}] //
    Log // PowerExpand // Simplify

mc = RandomWalkMetropolis[expr, {{\[Mu], \[Sigma]}, {0, 1}}, .1, 10^3,
    BurnInNumber -> 10^3, InBoundsQ -> (#[[2]] > 0 &)];

AcceptanceRate[mc]

Table[ListLinePlot[mc[[All, i]], PlotRange -> All], {i, 2}]

Table[Histogram[mc[[All, i]], {.05}], {i, 2}]

At its heart, the code runs NestList on a Metropolis step. If you use
the option ThinningNumber -> n (with n > 1), then the code acutally
runs NestList on Nest!

The code is pseudo-compiled, so its fast-ish. Is it fast enough? It's
fine for me. But all of my colleagues (economists) use either C
or some other econometrics-related software). 

Good luck.


--Mark


  • Prev by Date: Re: NEWBIE: How do I use the results of a Solve?
  • Next by Date: Re: Dynamic Module
  • Previous by thread: Re: MCMC in Mathematica
  • Next by thread: Re: MCMC in Mathematica