[ADMB Users] comparing highly parameterized models

Fowler, Mark Mark.Fowler at dfo-mpo.gc.ca
Fri May 3 07:17:30 PDT 2013


Thank you, Steve. I'll tackle that shortly (I just crashed - I don't
think ADMB and WinBUGS share toys very well). 

Is this the code on the TODO list for ISCAM (safe to add to ISCAM as
well)? 

>	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: Steve Martell [mailto:SteveM at iphc.int] 
Sent: May 3, 2013 10:56 AM
To: Fowler, Mark
Cc: John Sibert; <users at admb-project.org>
Subject: Re: [ADMB Users] comparing highly parameterized models

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%2Fre
>> p 
>> rints%2Fmultifan.pdf&ei=lKyCUc-lMYPHiwKk_YDQDA&usg=AFQjCNEW-tR28Npujf
>> A 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