[ADMB Users] comparing highly parameterized models
Steve Martell
SteveM at iphc.int
Fri May 3 06:55:45 PDT 2013
Hi Mark,
If you paste the following code in your GLOBALS_SECTION and run your model with the -mceval option after you've done your -mcmc run, then the on screen output will spit out the DIC and effective number of parameters.
Note that I've modified the mcmc_eval code from the source code.
void function_minimizer::mcmc_eval(void)
{
// |---------------------------------------------------------------------------|
// | Added DIC calculation. Martell, Jan 29, 2013 |
// |---------------------------------------------------------------------------|
// | DIC = pd + dbar
// | pd = dbar - dtheta (Effective number of parameters)
// | dbar = expectation of the likelihood function (average f)
// | dtheta = expectation of the parameter sample (average y)
gradient_structure::set_NO_DERIVATIVES();
initial_params::current_phase=initial_params::max_number_phases;
uistream * pifs_psave = NULL;
#if defined(USE_LAPLACE)
#endif
#if defined(USE_LAPLACE)
initial_params::set_active_random_effects();
int nvar1=initial_params::nvarcalc();
#else
int nvar1=initial_params::nvarcalc(); // get the number of active parameters
#endif
int nvar;
pifs_psave= new
uistream((char*)(ad_comm::adprogram_name + adstring(".psv")));
if (!pifs_psave || !(*pifs_psave))
{
cerr << "Error opening file "
<< (char*)(ad_comm::adprogram_name + adstring(".psv"))
<< endl;
if (pifs_psave)
{
delete pifs_psave;
pifs_psave=NULL;
return;
}
}
else
{
(*pifs_psave) >> nvar;
if (nvar!=nvar1)
{
cout << "Incorrect value for nvar in file "
<< "should be " << nvar1 << " but read " << nvar << endl;
if (pifs_psave)
{
delete pifs_psave;
pifs_psave=NULL;
}
return;
}
}
int nsamp = 0;
double sumll = 0;
independent_variables y(1,nvar);
independent_variables sumy(1,nvar);
do
{
if (pifs_psave->eof())
{
break;
}
else
{
(*pifs_psave) >> y;
sumy = sumy + y;
if (pifs_psave->eof())
{
double dbar = sumll/nsamp;
int ii=1;
y = sumy/nsamp;
initial_params::restore_all_values(y,ii);
initial_params::xinit(y);
double dtheta = 2.0 * get_monte_carlo_value(nvar,y);
double pd = dbar - dtheta;
double dic = pd + dbar;
dicValue = dic;
dicNoPar = pd;
cout<<"Number of posterior samples = "<<nsamp <<endl;
cout<<"Expectation of log-likelihood = "<<dbar <<endl;
cout<<"Expectation of theta = "<<dtheta <<endl;
cout<<"Number of estimated parameters = "<<nvar1 <<endl;
cout<<"Effective number of parameters = "<<dicNoPar <<endl;
cout<<"DIC = "<<dicValue <<endl;
break;
}
int ii=1;
initial_params::restore_all_values(y,ii);
initial_params::xinit(y);
double ll = 2.0 * get_monte_carlo_value(nvar,y);
sumll += ll;
nsamp++;
// cout<<sumy(1,3)/nsamp<<" "<<get_monte_carlo_value(nvar,y)<<endl;
}
}
while(1);
if (pifs_psave)
{
delete pifs_psave;
pifs_psave=NULL;
}
return;
}
On 2013-05-03, at 5:06 AM, "Fowler, Mark" <Mark.Fowler at dfo-mpo.gc.ca>
wrote:
> Thanks, John, I've got it. Slight misunderstanding, my subject header
> was inadequate. I'm okay for comparing models where the number of
> parameters are easily determined. My question concerned determination of
> the number of 'real' parameters represented by deviation parameters. Ian
> Taylor developed an approach to approximate the number of parameters for
> comparing Stock Synthesis models (the excerpt in my email), and I was
> asking for details on how a couple of scalars used by Ian were derived.
> Both Ian and Steve Martell responded, giving me a handle on them. I'll
> need to compute DIC's to estimate the effective number of parameters.
>
>> Mark Fowler
> Population Ecology Division
>> Bedford Inst of Oceanography
>> Dept Fisheries & Oceans
>> Dartmouth NS Canada
> B2Y 4A2
> Tel. (902) 426-3529
> Fax (902) 426-9710
> Email Mark.Fowler at dfo-mpo.gc.ca
> Home Tel. (902) 461-0708
> Home Email mark.fowler at ns.sympatico.ca
>
>
> -----Original Message-----
> From: users-bounces at admb-project.org
> [mailto:users-bounces at admb-project.org] On Behalf Of John Sibert
> Sent: May 2, 2013 4:03 PM
> To: users at admb-project.org
> Subject: Re: [ADMB Users] comparing highly parameterized models
>
> A less garbled link perhaps
> http://www.soest.hawaii.edu/PFRP/reprints/multifan.pdf
>
> John Sibert
> Emeritus Researcher, SOEST
> University of Hawaii at Manoa
> Honolulu HI (GMT-10)
> 808-294-3842
>
> Visit the ADMB project http://admb-project.org/
>
> On 05/02/2013 08:15 AM, dave fournier wrote:
>>
>> This is far too long to read, but we approached model selection in
>> the multifan cl model using approximate Bayes factors whihc enabled us
>
>> to copmoare non nested models.
>>
>> This link might have survived being copied.
>>
>>
>>
>> http://www.google.ca/url?sa=t&rct=j&q=&esrc=s&source=web&cd=4&cad=rja&
>> sqi=2&ved=0CEAQFjAD&url=http%3A%2F%2Fwww.soest.hawaii.edu%2Fpfrp%2Frep
>> rints%2Fmultifan.pdf&ei=lKyCUc-lMYPHiwKk_YDQDA&usg=AFQjCNEW-tR28NpujfA
>> kDNC5BDVIumdQhA&sig2=49idpQ7NXkxhrLbIYcULfg&bvm=bv.45960087,d.cGE
>>
>> _______________________________________________
>> Users mailing list
>> Users at admb-project.org
>> http://lists.admb-project.org/mailman/listinfo/users
>>
>
> _______________________________________________
> Users mailing list
> Users at admb-project.org
> http://lists.admb-project.org/mailman/listinfo/users
> _______________________________________________
> Users mailing list
> Users at admb-project.org
> http://lists.admb-project.org/mailman/listinfo/users
| ---------------------- |
| steven martell |
| stevem at iphc.int |
| (206) 552-7683 |
| ---------------------- |
________________________________
This internet e-mail message, and any files transmitted with it, contains confidential, privileged information that is intended only for the addressee. If you have received this e-mail message in error, please call us at (206) 634-1838 collect if necessary) and ask to speak to the message sender. Nothing in this e-mail or the act of transmitting it, is to be construed as a waiver of any rights or privileges enjoyed by the sender or the International Pacific Halibut Commission pursuant to the International Organizations Immunities Act, 22 U.S.C. Sec. 288 et seq.
More information about the Users
mailing list