[ADMB Users] big(ger) models
dave fournier
davef at otter-rsch.com
Fri May 9 13:27:37 PDT 2014
I noticed that someone had posted on the SAS list about a moderate
sized model
requiring something like 80 hours to fit using GLIMMIX. I thought this
might be
an interesting application of some of the sparse matrix stuff. The
poster had
posted a simple example of the kind of problem he was having trouble with.
http://listserv.uga.edu/cgi-bin/wa?A2=ind1405b&L=sas-l&P=4001
He also later said that his real problem has an ar(1) structure which we
know how to deal with.
(I think). I wrote the attached TPL file to test the ideas.
I ran the model with the following.
time ./ryan -mno 1000000 -ams 10000000 -ilmn 5 -shess -ndi 100000
-noinit -l1 100000000 -nl1 100000000 -iprint 1 -crit 1.e-3 -noreders
the -noreders is to prevent the calculation of the estimates of the std
devs of the random effects via the delta method.
This saves a lot of time and space since there about 75,000 random
effects. Fitting the model required less than 4 minutes.
I don't know how many ADMB people are interested in this kind of
larger model.
It appears that the -noreders option has never been included in the
code. It is only about 8 or so lines of code to add.
Maybe someone is interested. Per usual I will hold my breath.
-------------- next part --------------
DATA_SECTION
int nloc
int nsub
int nobs
int iseed
LOC_CALCS
nloc=250;
nsub=300;
nobs=4;
END_CALCS
3darray y(1,nloc,1,nsub,1,nobs)
matrix x(1,nloc,1,nsub)
LOC_CALCS
iseed=12345;
random_number_generator rng(iseed);
double lre;
double sre;
double prob;
double eta;
for (int i=1;i<=nloc;i++)
{
lre=sqrt(.3)*randn(rng);
for (int j=1;j<=nsub;j++)
{
sre=sqrt(.5)*randn(rng);
x(i,j)=randn(rng);
eta=1.6+2.5*x(i,j)+lre+sre;
prob=exp(eta)/(1.0+exp(eta));
for (int k=1;k<=nobs;k++)
{
double z=randu(rng);
if (z<=prob)
y(i,j,k)=1.0;
else
y(i,j,k)=0.0;
}
}
}
PARAMETER_SECTION
objective_function_value f
init_number a
init_number b
init_bounded_number lsd(0.0,5.0,2)
init_bounded_number ssd(0.0,5.0,2)
!! lsd.set_scalefactor(100.);
!! ssd.set_scalefactor(100.);
random_effects_vector eloc(1,nloc,2)
random_effects_matrix esub(1,nloc,1,nsub,2)
PROCEDURE_SECTION
for (int i=1;i<=nloc;i++)
{
for (int j=1;j<=nsub;j++)
{
f1(i,j,a,b,lsd,ssd,eloc(i),esub(i,j));
}
}
SEPARABLE_FUNCTION void f1(int i,int j,const prevariable& a,const prevariable& b,const prevariable& lsd,const prevariable& ssd,const prevariable& el,const prevariable& es)
dvariable eta=a+b*x(i,j)+lsd*el+ssd*es;
dvariable prob;
if (value(prob)<20.)
prob=exp(eta)/(1.0+exp(eta));
else
prob=1.0/(1.0+exp(-eta));
for (int k=1;k<=nobs;k++)
{
if (y(i,j,k)>0)
f-=log(1.e-10+prob);
else
f-=log(1.e-10+(1.0-prob));
}
if (j==1)
{
f+=0.5*square(el);
}
f+=0.5*square(es);
More information about the Users
mailing list