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