[ADMB Users] Fwd: admb10 64bit has problem?
Ben Bolker
bbolker at gmail.com
Fri Feb 18 13:46:53 PST 2011
On 32-bit Ubuntu 10.04 with the latest SVN ADMB I see problems with
everything up to obs=11; thereafter the three approaches match.
*catch,mu,s 0 1 1 nll 0.5 0.5 1
*catch,mu,s 1 1 1 nll 0.25 1 1
*catch,mu,s 2 1 1 nll 0.125 0.25 1
*catch,mu,s 3 1 1 nll 0.0625 0.0833333 1
*catch,mu,s 4 1 1 nll 0.03125 12 0.0416667
*catch,mu,s 5 1 1 nll 0.015625 2.26508e+26 0.00833333
*catch,mu,s 6 1 1 nll 0.0078125 2.28229e+26 4.3816e-27
*catch,mu,s 7 1 1 nll 0.00390625 0.000539342 4.3787e-27
*catch,mu,s 8 1 1 nll 0.00195312 2.48016e-05 2.28399e+26
*catch,mu,s 9 1 1 nll 0.000976562 389749 2.28401e+26
*catch,mu,s 10 1 1 nll 0.000488281 2.28401e+26 2.28401e+26
*catch,mu,s 11 1 1 nll 0.000244141 2.28401e+26 2.28401e+26
As far as I can tell it seems to be something screwy in ADMB, not just
a numerically unstable calculation; I tried translating the first two
versions into plain C code (attached) and the results look just fine.
Looking more closely at what ADMB actually does explains why
log_negbinomial_density blows up (I'm surprised that it works reliably
on any platform for s=1!) -- it calculates an intermediate parameter r =
mu/(1e-120 + (s-1)). When tau happens to equal 1, the calculated
variable r is going to be 1e120, and the result one gets for tmp will
depend entirely on what order one adds the terms in (and hence which
significant digits fall off the table: 1e120 + any value less than about
1e104 is equal to 1e120 in double precision ...
This certainly seems like a nasty trap. I can imagine the results
would be bad (if not completely pathological) for values of s *near*
(not exactly equal to) 1. I tried s=1.001 in the example above and got
wacky answers from log_negbinomial_density for x= 9 to 11 .
At face value it seems like Weihai's function is much better, but of
course I can imagine it will fail in other extreme cases (e.g. when
s>>mu or mu>>s).
cheers
Ben Bolker
df1b2variable log_negbinomial_density(double x, const df1b2variable & _xmu,
const df1b2variable & _xtau)
{
ADUNCONST(df1b2variable, xmu)
ADUNCONST(df1b2variable, xtau) init_df3_two_variable mu(xmu);
init_df3_two_variable tau(xtau);
*mu.get_u_x() = 1.0;
*tau.get_u_y() = 1.0;
if (value(tau) - 1.0 < 0.0)
{
cerr << "tau <=1 in log_negbinomial_density " << endl;
ad_exit(1);
}
df3_two_variable r = mu / (1.e-120 + (tau - 1.0));
df3_two_variable tmp;
tmp = gammln(x + r) - gammln(r) - gammln(x + 1)
+ r * log(r) + x * log(mu) - (r + x) * log(r + mu);
df1b2variable tmp1;
tmp1 = tmp;
return tmp1;
}
On 02/18/2011 03:15 PM, Johnoel Ancheta wrote:
>
> On Fri, 18 Feb 2011, Weihai Liu wrote:
>
>>
>
>>
> Hi Arni,
>
>>
>>
> Currently I found out there are some unconsistent about using
> log_negbinomial_density() function for ADMB 10 64 bit.
>
>>
>>
> Following is my testing tpl, it runs ok on my windows admb9 32bit. 3
> variant function can agree. But not consistent on ubuntu admb10 64
> bit, I am quickly test on your new release windows 64 bit, it has the
> same issue, not consistent.
>
>
>>
>>
> The testing function is nll for NB(m,s), which use gammln() built in
> function, why I build my own version, because at first I found the
> ADMB builtin log_negbinomial_density() failed on linux admb10 64 bit,
> so to figure out what cause the problem, I set it up two other
> alternative functions,see below.
>
>
>>
>>
> From my own nllNegativeBinomial(), I have to separate some terms out,
> then the problem gone. I feel there are some issue on admb10 64 bit,
> or some issue come from gammln(), since I am not testing it on admb10
> 32 bit on windows,so not sure the problem is from different version
> of admb or different bits or combinations.
>
>
>>
>>
> you can play around by changing the mu,s in NB(mu,s), or get rid of
> mfexp in the cout to see nll.
>
>>
>
>>
>>
> GLOBALS_SECTION
>>
> #include "admodel.h"
>>
> const double EPS = 1.e-20;
>
>>
>>
> DATA_SECTION
>
>>
>>
> PARAMETER_SECTION
>>
> init_number junk
>>
> //random_effects_vector junk1(1,2)
>>
> objective_function_value nll
>
>>
>>
> PROCEDURE_SECTION
>>
> dvariable mu=1.;
>>
> dvariable s=1.;
>
>>
>>
> for(int n=0;n<=20;n++){
>>
> //full negative log likelihood for NB(mu,s)
>>
> cout<<"*catch,mu,s "<<n<<" "<<mu<<" "<<s<<" nll ";
>>
> cout<< mfexp(-nllNegativeBinomial(double(n),mu,s)); //works both on
>>
> admb9(win) 32bit,admb10(linux) 64bit
>>
> cout<<" "<<mfexp(-nllNegativeBinomial2(double(n),mu,s)); //good on
>>
> admb9(win) 32bit,fail on admb10(linux) 64bit
>>
> cout<<"
> "<<mfexp(log_negbinomial_density(double(n),mu,1.+mu/s))<<endl;
>>
> //good on admb9(win)32bit,fail on admb10(linux) 64bit
>>
> }
>>
> exit(8);
>
>>
>
>>
>>
> //following version get consistent results for windows admb9 32bit
> and linux
>>
> admb10 64bit
>>
> FUNCTION dvariable nllNegativeBinomial(double obs,const dvariable&
> m,const
>>
> dvariable& s)
>>
> dvariable nll1=0;dvariable nll2=0; dvariable nll3=0; dvariable
> nll=0;
>>
> //nll= -gammln(obs+s)+ gammln(s)- s*log(s/(m+s))-
> obs*log(m/(m+s)+EPS)+
>>
> gammln(obs+1.);
>>
> nll1= -gammln(obs+s)+ gammln(s)+gammln(obs+1.);
>>
> nll2= - s*log(s/(m+s)+EPS);
>>
> nll3= - obs*log(m/(m+s)+EPS);
>>
> nll=nll1+nll2+nll3;
>>
> return nll;
>
>>
>>
> //following version works for windows admb9 32bit, same as builtin
>>
> log_negbinomial_density()
>>
> //not good on linux admb10 64bit, failed as
> log_negbinomial_density(),but
>>
> they got somewhat different results
>>
> FUNCTION dvariable nllNegativeBinomial2(double obs,const dvariable&
> m,const
>>
> dvariable& s)
>>
> dvariable nll=0;
>>
> nll= -gammln(obs+s)+ gammln(s)- s*log(s/(m+s)+EPS)-
> obs*log(m/(m+s)+EPS)+
>>
> gammln(obs+1.);
>>
> //nll= -gammln(obs+s)+ gammln(s)- s*(log(s)-log(m+s))-
>>
> obs*(log(m)-log(m+s))+ gammln(obs+1.);
>>
> return nll;
>
>>
>>
> weihai
>
>>
>
>
>>
>>>
>> -- Weihai Liu ~~~~~~~~~~~~~~~~~~~~~~~~ Quantitative Fisheries
>> Center Department of Fisheries & Wildlife Michigan State
>> University 153 Giltner Hall East Lansing, MI 48824
>>
> Phone: 517-355-0126
>> Email: liuweih at msu.edu http://qfc.fw.msu.edu/weihai
>> ~~~~~~~~~~~~~~~~~~~~~~
>
>
> _______________________________________________ Users mailing list
> Users at admb-project.org
> http://lists.admb-project.org/mailman/listinfo/users
-------------- next part --------------
A non-text attachment was scrubbed...
Name: nbtest.c
Type: text/x-csrc
Size: 1387 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/users/attachments/20110218/6722216c/attachment.c>
More information about the Users
mailing list