[Developers] MCMC not working in ADMB10

[anders at nielsensweb.org]

Hi Allan,

Great job finding this, and narrowing it down. Kasper and I put in
the replacement code for the random number generator (Mersenne
Twister), and did some testing, but clearly not enough.

Kind of strange that nobody has had trouble with this before now, but
I guess relatively few run the code in the svn repository, and of those
very few use the mcmc stuff.

One thing that makes it difficult to fix this in an optimal way is that
the previous rng in ADMB could summarize its current state in a single
number, but the MT cannot it needs an array.

So the thing that is happening in the bet_rand.cpp file, where a new rng
is started at the state saved in idum is difficult to duplicate, I guess
one would have to use the state-array in place of idum:

double better_rand(long int& idum)
{
  random_number_generator rng(idum);
  double rr = ((random_number_generator&) rng).better_rand();
  return rr;
}

One thing that we could do easily is to simply draw a new random seed
within this version of better_rand, so that the next rng defined would
be different, but I worry about the "random" properties, when not using
the rng in sequence. The above could be done by:

double better_rand(long int& idum)
{
  random_number_generator rng(idum);
  idum=(long int)((((random_number_generator&)
rng).better_rand())*2147483647.0);
  double rr = ((random_number_generator&) rng).better_rand();
  return rr;
}

(the idum=... ; should all be one line)

I tested it on a simple problem and it seemed to fix the problem,
but before submitting it I would like to hear Dave, or some other
expert's opinion.

We could also switch back to numerical recipes generator.

Cheers,

Anders.

PS: Another possibly related issue is the speed of MCMC. To me it
seems that it starts out fast and then gets slower and slower. To
my memory that was not always the case.




> Thanks Dave! Yes, the random number generator is the issue.  The
> function better_rand was replaced due to the issues with Numerical
> Recipes copyright, I believe.  I tested this by replacing the
> bet_rand.cpp file with the ADMB9 bet_rand.cpp file and the mcmc appears
> to work (at least the acceptance rate is reasonable). The random number
> generated in ADMB 10 was always constant at 0.937887 for my example
> (changed with -mcseed, though).
>
> I looked into the random number generator (in rngen.cpp) and it matches
> the Matsumota & Nishimura paper. But when that function is called from
> the bet_rand.cpp file, it is always using the same seed because a new
> instance of the random_number_generator class is created with the same
> seed (iseed) in the call from xxmcmc.cpp (line 744). It seems to me that
> what needs to be done is the class "rng" needs to be declared once,
> outside of the better_rand function in bet_rand.cpp.
>
> This is as far as I got, and I think it is best left up to the people
> that understand the code best to decide how the seed should be set only
> once.  But, at least we know our results with random numbers can be
> reproduced with the same seed.
>
> Great job with the new random number generators, and thanks again to
> Dave for instantly seeing what this problem likely was.
>
> Allan
>
>
>
> ----- Original Message -----
> From: dave fournier <davef at otter-rsch.com>
> Date: Thursday, November 18, 2010 8:52 am
> Subject: Re: [Developers] MCMC not working in ADMB10
>
>> Allan Hicks wrote:
>>
>> I think there is a problem with the random number generator. Maybe
>> it
>> got changed?
>> > Hi All,
>> >
>> > I wish that I had the time to present a solution to this problem,
>> but
>> > unfortunately I am only identifying the issue.
>> >
>> > It appears that the MCMC algorithm is not working in ADMB 10
>> beta.
>> > When an MCMC is run, the acceptance rate goes to 0 or 1, which
>> seems
>> > to depend where it left off during the scaling phase.  Attached
>> are
>> > screen outputs from the catage example run in ADMB10 and ADMB9.
>> You
>> > can see that ADMB10 also runs out of memory after 21800
>> iterations.
>> > ADMB9 appears to have run the MCMC fine.
>> >
>> > This brings up a point I have been thinking about. Should we
>> include
>> > testing of MCMC results with new versions of ADMB? This is
>> something
>> > we are just beginning to discuss with regard to testing new
>> versions
>> > of SS3.
>> >
>> > Thanks for all your help,
>> > Allan
>> > ------------------------------------------------------------------
>> ------
>> >
>> > _______________________________________________
>> > Developers mailing list
>> > Developers at admb-project.org
>> > http://lists.admb-project.org/mailman/listinfo/developers
>>
>>
> _______________________________________________
> Developers mailing list
> Developers at admb-project.org
> http://lists.admb-project.org/mailman/listinfo/developers
>