Re: [Help-gsl] quasi-Monte Carlo integration

[David Zaslavsky]

Hi Rhys,

Sorry, it completely slipped my mind to respond to your email earlier.

I definitely considered that option before writing my own code. It would
be a little tricky at best, though, because a quasirandom number
generator produces a tuple of N values at a time (where N is the
dimension of the integration region), whereas the Monte Carlo methods
request random numbers individually. See around line 77 of monte/plain.c
compared to around line 75 of my quasimontecarlo.c.

You could, in principle, get around that by having the wrapper store the
values produced from the gsl_qrng it wraps and return them from
gsl_rng_uniform_pos one at a time, only calling the underlying gsl_qrng
when necessary to get new values. But then I'm not sure what would be
sensible behavior if any of the other sampling functions (gsl_rng_get,
gsl_rng_uniform, gsl_rng_uniform_int) were called. Possibly some kind of
error would be appropriate, because consuming a value and continuing
would throw off the QMC sampling, but that would make the gsl_qrng
wrapper fairly useless for anything other than plain QMC integration.
Now, sure, it's a valid solution to just say so in the documentation,
but it seems inelegant to me, and since this is a constraint that can be
enforced by the API (i.e. that a gsl_qrng, even wrapped, can only be
sensibly used for plain Monte Carlo), I think it makes sense to do so.

Someday when I have the time (not any time soon, unfortunately), I'll
look into implementing adaptive QMC integration, and perhaps that will
suggest some better way that the QMC and MC code could be merged. I
definitely agree that *if* there is a sensible way to do so, it's
preferable to having separate QMC and MC code. But it's a more subtle
issue than it appears at first.

:) David

On 11/11/2013 11:13 AM, Rhys Ulerich wrote:
> Hi David,
> 
>> It occurred to me today that I should have sent this out some time ago,
>> but anyway... unless I'm missing something, I don't believe GSL natively
>> includes any way to do quasi-Monte Carlo integration, and neither do any
>> of the listed extensions allow for it. So I wrote my own, at
>>  https://github.com/diazona/quasimontecarlo
>> It is just a copy of the (plain) GSL Monte Carlo integrator that uses a
>> quasirandom number generator instead of a pseudorandom number generator.
> 
> Thank you.  Is there some way that we could wrap a gsl_qrng into
> something providing a gsl_rng interface?  This would permit handing
> the entire Quasi MC problem off to the normal MC code.
> 
> - Rhys
>

signature.asc
Description: OpenPGP digital signature