Abel 1068 Event-2D :
Fitting and MCMC

Back to Event-2D Demos page

10/4/08
dd@space.mit.edu

This page shows the effects of the commands in a1068_cg_mcmc.sl.
The command window text is given below along with select plots or images that are created.
The user-typed input is shown in bold and is assumed to be input at the "hydra>" prompt.
Commentary is shown in green.

% - - - Read in the data:

.source a1068_data
 Reading file:  ../Data/obs_1652/evt1_s3.fits
 - - -
 Reading file:  ../Data/obs_1652/evt1_s3.fits
 - - -
 - Data:                                     e2d_meta[id].;  e2d_data[id].;  e2d_view[id].
  id  Mission-grat-instrument    Date         2D-axes     ignore  prt  m src   events(k)   exp(ks)            name  
   1        CHANDRA-NONE-ACIS  2001.095            X,Y       0     -   -  -       74.24      26.8         ABELL 1068
>  2        CHANDRA-NONE-ACIS  2001.095  RADIUS,ENERGY       0     -   -  -       74.24      26.8         ABELL 1068
 - - -

% Includes both X-Y and RADIUS-ENERGY data sets.
We'll fit just the X,Y Image, so
ignore the RADIUS-ENERGY data:

e2d_ignore(2);

% Look at the X,Y data:

e2d_view_data(1); 
   e2d_view_data ( [ -90 , 90 , 2.5 ] , [ -90 , 90 , 2.5 ] ); 
     total events displayed: 53420     e2d_view[e2di].scale = 555%
     bins w/events have 50%, 90%, and max numbers of:  4,  22,  626

% - - - Setup the model:

.source a1068_model_simple
 Spectrum created from: temp_a1068_T32A05.par
 - - -
      o check/update foreground NH ...
      o check/update spectra ...
      o check/update the component geometry ...
      o check/update opaque types ...

       --- Source-3D Model ---

 Model name      : A1068 simple model
 Distance        : 505000 kpc.               Date : 2001.000
 Foreground NH   : 0.022 x10^22 /cm^2
 Cube radius     : 100 arc sec  with 37 cells/radius, cell size of 2.7027 arc sec.
 - - -
 s3d_NHarf_method = 1  and 'slope = 3
 Nominal spectral grid : 1.45864 A  to 35.4241 A,  with delta-log of 0.0003
 - - -
 - Spectra:                                                                  s3d_spec[id].
  id    min - max es   ph/ks/100cm^2    file
   2     0.35 - 8.50     52130.7        temp_a1068_T32A05.par
 - - -
 - Components:                                                               s3d_comp[id].
  id         name   ignore   type       norm   ispec   vtype  velval   vturb    rndint          rndclr
   1    Beta model     0     xray        0.01    2      -1         0       0        1          [1,1,0]
 - - -
 The s3d_ps user parameter structure values:
{xc=0, yc=0, ratio=1.4, ang=-42.5, beta=0.5, rbeta=12.6}
 - - -

% Set some rough parameter values:
first look at the available parameters and their values:

print(s3d_ps);
{xc=0, yc=0, ratio=1.4, ang=-42.5, beta=0.5, rbeta=12.6}

% Set these to values we might have guessed:

s3d_ps.rbeta = 15.0;
s3d_ps.ratio = 1.2;
s3d_ps.ang = -30.0;
s3d_ps.beta = 0.50;

% and look at a projection view of the components:

s3d_update;
      o check/update foreground NH ...
      o check/update spectra ...
      o check/update the component geometry ...
      o check/update opaque types ...
s3d_proj_comps;

% - - - Setup the binning and science sigma % Evaluate the model

e2d_pars2fom;
         Event2D:   parameters to FOM
   - update parameters                                                  [e2d_update_pars]
   - updating the model ...                                                  [s3d_update]
      o check/update foreground NH ...
      o check/update spectra ...
      o check/update the component geometry ...
      o check/update opaque types ...
   - Dataset 1: generating photons from the model...                   [s3d_xray_photons]
      o Beta model: generating 291724 MC photons...
      o Total generated =  291724
      o Total kept after NH/arf =  221539     ( 75.94 % )
             1: pass photons through instrument simulation...             [e2d_inst_xray]
             1: load simulation into model arrays                          [e2d_inst2mod]
   - calculate and show the figure of merit...                              [e2d_list_fom]
      o Data-Model 1,  chi^2 = 8603.27  from 4712 cells;   chi^2/cell = 1.82582;  sum(D)/sum(M) = 1.03418  +/- 0.00447451
      o Total FOM = 8603.27
 - - -

% look at the residuals in 2D Images are: Data --- Model --- Residual

e2d_view_resid(1);

% change to a coarser binning for fitting:

e2d_view_data (1, [ -90 , 90 , 7 ] , [ -90 , 90 , 7 ] );
   e2d_view_data ( [ -90 , 90 , 7 ] , [ -90 , 90 , 7 ] ); 
     total events displayed: 53420     e2d_view[e2di].scale = 1538%
     bins w/events have 50%, 90%, and max numbers of:  29,  159,  3126
 - - -
e2d_list_fom;
      o Data-Model 1,  chi^2 = 7047.34  from 676 cells;   chi^2/cell = 10.4251;  sum(D)/sum(M) = 1.03418  +/- 0.00447451
      o Total FOM = 7047.34

% Look at the Data-Model residuals in 2D

e2d_view_resid(1);

% Check the "chi-percentage" plot

e2d_chipc_plot;

% Change to a 10% science-sigma limit:

e2d_meta[1].scisigma = 10.0e-2;
e2d_list_fom;
      o Data-Model 1,  chi^2 = 2948.09  from 676 cells;   chi^2/cell = 4.36108;  sum(D)/sum(M) = 1.03418  +/- 0.00447451
      o Total FOM = 2948.09
 - - -

% And look at the "chi-percentage" plot now

e2d_chipc_plot;

% Check that the norm is roughly correct and adjust it:

norm_adj = e2d_norm_optimum(NULL);
    - Minimum FOM of 1847.57 at norm factor = 0.796832
 - - -

% Scale the norm for better fit:

s3d_comp[1].norm = s3d_comp[1].norm * norm_adj;

% update the model and look at our current residuals:

e2d_pars2fom;
         Event2D:   parameters to FOM
   - update parameters                                                  [e2d_update_pars]
   - updating the model ...                                                  [s3d_update]
      o check/update foreground NH ...
      o check/update spectra ...
      o check/update the component geometry ...
      o check/update opaque types ...
   - Dataset 1: generating photons from the model...                   [s3d_xray_photons]
      o Beta model: generating 232455 MC photons...
      o Total generated =  232455
      o Total kept after NH/arf =  176404     ( 75.88 % )
             1: pass photons through instrument simulation...             [e2d_inst_xray]
             1: load simulation into model arrays                          [e2d_inst2mod]
   - calculate and show the figure of merit...                              [e2d_list_fom]
      o Data-Model 1,  chi^2 = 1778.95  from 676 cells;   chi^2/cell = 2.63158;  sum(D)/sum(M) = 1.29902  +/- 0.00562036
      o Total FOM = 1778.95
 - - -
e2d_view_resid(1);

% - - - Scan the parameters
To get some idea of the scale over which each parameter changes,
do 1-D scans of each parameter with the others fixed.

Set the parameters close to where the scan has its
minimum or low point and make a note of the
"should make a noticable difference" scale for each parameter.

Each scan requires the "e2d_update_pars" routine be
redefined for the appropriate scan, and then
the e2d_scan_par routine can be run.

reset the FOM-log values:
(The scanned parameters and chi-squareds are saved for
viewing - this resets the save list.)


e2d_fomlog_i = 0;

% Scan over values of rbeta:

public define e2d_update_pars () 
{         
   s3d_ps.rbeta =  e2d_scan_vals[e2d_iscan];
}                 
e2d_scan_par([8.0:18.0:#9],3);
 ====================
         Event2D:   parameters to FOM
   - update parameters                                                  [e2d_update_pars]
   - updating the model ...                                                  [s3d_update]
      o check/update foreground NH ...
      o check/update spectra ...
      o check/update the component geometry ...
      o check/update opaque types ...
   - Dataset 1: generating photons from the model...                   [s3d_xray_photons]
      o Beta model: generating 232455 MC photons...
      o Total generated =  232455
      o Total kept after NH/arf =  176387     ( 75.88 % )
             1: pass photons through instrument simulation...             [e2d_inst_xray]
             1: load simulation into model arrays                          [e2d_inst2mod]
   - calculate and show the figure of merit...                              [e2d_list_fom]
      o Data-Model 1,  chi^2 = 2050.84  from 676 cells;   chi^2/cell = 3.03378;  sum(D)/sum(M) = 1.28308  +/- 0.00555138
      o Total FOM = 2050.84
 - - -
   - Minimum FOM of 1709.27 at norm factor = 1.15549
 - - -
         ...
     This takes a couple of minutes !
         ...
 - - -
   - FOM = 2069.66
   - Scan parameter value was: 18
======   ======
 
 ====================
 Minimum FOM is at scan value = 10.5   giving a FOM of 1603.82  +/- 23.0599
    with a best norm adjust of 1.09113
    *** Parameter(s) set to this best-FOM value, BUT model not re-evaluated. ***
 Fitting gives minimum at scan value = 9.90345 +/- 0.683506,  with a FOM of 1625.85

 ====================

% Set it to ~ min and note that smand ~ 3.0

s3d_ps.rbeta=10.5;

% Scan over values of ratio:

public define e2d_update_pars () 
{         
   s3d_ps.ratio = e2d_scan_vals[e2d_iscan]; 
}                 
e2d_scan_par([1.0:1.8:#9],3);
         ...
     This takes a couple of minutes !
         ...

% Set it to ~ min and note that smand ~ 0.2

s3d_ps.ratio = 1.35;

% Scan over values of angle: * Note: nominal value is negative ! *

public define e2d_update_pars () 
{         
   s3d_ps.ang =  e2d_scan_vals[e2d_iscan];
}                 
e2d_scan_par([-60.0:0.0:#9],3);
         ...
     This takes a couple of minutes !
         ...

% Set it to ~ min and note that smand ~ 15.0

s3d_ps.ang = -40.0;

% Scan over values of beta:

public define e2d_update_pars () 
{         
   s3d_ps.beta =  e2d_scan_vals[e2d_iscan];
}                 
e2d_scan_par([0.40:0.60:#9],3);
         ...
     This takes a couple of minutes !
         ...

% Set it to ~ min and note that smand ~ 0.04

s3d_ps.beta = 0.50;

% Reset/cancel the definition of e2d_update_pars :
( It is good to be sure this is done habitually
after scanning or fitting to cancel the
e2d_update_pars definition. )

e2d_update_pars_NULL;

% Look at our new scan-fit parameter values
and evaluate the model etc as before:

s3d_ps.rbeta=10.5;
s3d_ps.ratio = 1.35;
s3d_ps.ang = -40.0;
s3d_ps.beta = 0.50;

e2d_pars2fom;
         Event2D:   parameters to FOM
   - update parameters                                                  [e2d_update_pars]
   - updating the model ...                                                  [s3d_update]
      o check/update foreground NH ...
      o check/update spectra ...
      o check/update the component geometry ...
      o check/update opaque types ...
   - Dataset 1: generating photons from the model...                   [s3d_xray_photons]
      o Beta model: generating 232455 MC photons...
      o Total generated =  232455
      o Total kept after NH/arf =  176417     ( 75.89 % )
             1: pass photons through instrument simulation...             [e2d_inst_xray]
             1: load simulation into model arrays                          [e2d_inst2mod]
   - calculate and show the figure of merit...                              [e2d_list_fom]
      o Data-Model 1,  chi^2 = 1419.11  from 676 cells;   chi^2/cell = 2.09927;  sum(D)/sum(M) = 1.28916  +/- 0.00557771
      o Total FOM = 1419.11
 - - -
norm_adj = e2d_norm_optimum(NULL);
   - Minimum FOM of 1275.11 at norm factor = 1.09611
 - - -

% Scale the norm (by ~ 1.09) for better fit:

s3d_comp[1].norm = s3d_comp[1].norm * norm_adj;
e2d_pars2fom;
         Event2D:   parameters to FOM
   - update parameters                                                  [e2d_update_pars]
   - updating the model ...                                                  [s3d_update]
      o check/update foreground NH ...
      o check/update spectra ...
      o check/update the component geometry ...
      o check/update opaque types ...
   - Dataset 1: generating photons from the model...                   [s3d_xray_photons]
      o Beta model: generating 254796 MC photons...
      o Total generated =  254796
      o Total kept after NH/arf =  193434     ( 75.91 % )
             1: pass photons through instrument simulation...             [e2d_inst_xray]
             1: load simulation into model arrays                          [e2d_inst2mod]
   - calculate and show the figure of merit...                              [e2d_list_fom]
      o Data-Model 1,  chi^2 = 1287  from 676 cells;   chi^2/cell = 1.90385;  sum(D)/sum(M) = 1.17807  +/- 0.00509706
      o Total FOM = 1287
 - - -

% These scan-set values are an improvement over the rough values, giving chi-squared (FOM) of ~ 1300 instead of ~ 1800.

e2d_view_resid(1);

% From these residuals it's clear that the center is
unlikely to be well-fit with a simple beta model only,
e.g., we should add a central point source, etc.
However, let's go with this model and fit it and
MCMC it for demonstration purposes.

- - - Fit the paramters

Each fit requires the "e2d_update_pars" routine be
redefined for the appropriate fit, and then
the e2d_cg_fit routine can be run.
(this looks similar to the scan case because the
cg fit routine calls the scan routine.)

Set up to fit 4 parameters:
Assign them to the cg_p0 array in the order that
we scanned them above for convenience:

public define e2d_update_pars () 
{ 
   % this line is the same for all cg fits: 
   variable pvals = cg_p0 + cg_smand*cg_dir*e2d_scan_vals[e2d_iscan]; 
   % map the fitting parameter values to the desired model params: 
   s3d_ps.rbeta = pvals[0]; 
   s3d_ps.ratio = pvals[1]; 
   s3d_ps.ang = pvals[2]; 
   s3d_ps.beta = pvals[3]; 
   % can also clamp parameters, etc. here:
   if (s3d_ps.rbeta < 3.0) s3d_ps.rbeta = 3.0;
}

% Of course it would be better to start the fitting at the scan
values we found above, but let's make the cg routine
show us how it works. (We could have skipped the scanning
if we had some idea of the "smand" scales to use.)

Start at the same very rough guess set of values as
we used above before the scanning.

cg_p0 = [15.0, 1.20, -30.0, 0.50];

% and set the "should make a noticable difference" values
for the parameters according to what the scans suggested:

cg_smand = [3.0, 0.2, 15.0, 0.04];

% === do the fit: (the "1" is a dummy argument)

e2d_cg_fit(1);

This fitting takes ~ 30 minutes ...

The fitting process is recorded in an HTML log file
with links to plots of the scans made while fitting:

cg_log_Fri_Oct_3_13_48_25_2008.html

% Plot the FOM-log values
Show the delta chi-squared that corresponds to
Data-Model agreement = sqrt(2*N_dof)

variable delta_chi = sqrt(2.0*e2d_data[1].nfom);
delta_chi;
e2d_fomlog_plot("rbeta","beta", delta_chi);

% Make the same kind of plot for the ang[le] and ratio parameters:

e2d_fomlog_plot("ang","ratio", delta_chi);

% Look at the data and model and residual for these "best fit" values:

e2d_update_pars_NULL;
s3d_ps.rbeta = 7.832;
s3d_ps.ratio = 1.419;
s3d_ps.ang = -41.73;
s3d_ps.beta = 0.4897;
e2d_pars2fom;
         Event2D:   parameters to FOM
   ... etc ...
   - calculate and show the figure of merit...                              [e2d_list_fom]
      o Data-Model 1,  chi^2 = 1288.96  from 676 cells;   chi^2/cell = 1.90675;  sum(D)/sum(M) = 1.17748  +/- 0.00509451
      o Total FOM = 1288.96
 - - -
norm_adj = e2d_norm_optimum(NULL);
   - Minimum FOM of 1249.79 at norm factor = 1.04778
 - - -
s3d_comp[1].norm = s3d_comp[1].norm * norm_adj;
e2d_pars2fom;
         Event2D:   parameters to FOM
   ... etc ...
   - calculate and show the figure of merit...                              [e2d_list_fom]
      o Data-Model 1,  chi^2 = 1208.38  from 676 cells;   chi^2/cell = 1.78754;  sum(D)/sum(M) = 1.11793  +/- 0.00483685
      o Total FOM = 1208.38
 - - -

% And look at the residuals for this case:

e2d_view_resid(1);

% For fun make a plot of the model projection for these parameters (on left below)
as compared to the model from the rough parameter values (on right, repeated from above):

s3d_proj_comps;

% - - - Do MCMC exploration

Uses the same e2d_update_pars routine and the
cg_p0 and cg_smand parameters already defined above -
this makes it easy to do MCMC and CG fitting together/alternately.


public define e2d_update_pars () 
{ 
   % this line is the same for all cg fits: 
   variable pvals = cg_p0 + cg_smand*cg_dir*e2d_scan_vals[e2d_iscan]; 
   % map the fitting parameter values to the desired model params: 
   s3d_ps.rbeta = pvals[0]; 
   s3d_ps.ratio = pvals[1]; 
   s3d_ps.ang = pvals[2]; 
   s3d_ps.beta = pvals[3]; 
   % can also clamp parameters, etc. here:
   if (s3d_ps.rbeta < 3.0) s3d_ps.rbeta = 3.0;
}

% Start the MCMC at the same point as the cg_fit did:

cg_p0 = [15.0, 1.20, -30.0, 0.50];
cg_smand = [3.0, 0.2, 15.0, 0.04];

% Put that same delta-chi here to help scale delta-chisquared
to a probability, like simulated annealing with "T" = sqrt(2*N_dof)

e2d_mcmc_explore(100, 1.0/delta_chi);
 ====================
   ... etc ...
 ====================

% start at the fit min:

cg_p0 = [7.832, 1.419, -41.73, 0.4897];
e2d_mcmc_explore(100, 1.0/delta_chi);
 ====================
   ... etc ...
 ====================

% start between them: :-)

cg_p0 = 0.5*( [7.832, 1.419, -41.73, 0.4897] + [15.0, 1.20, -30.0, 0.50] );

% and "lower the temp" by half to emphasize the minimums:

e2d_mcmc_explore(100, 2.0/delta_chi);
 ====================
   ... etc ...
 ====================

% start at the fit min: and T/3

cg_p0 = [7.832, 1.419, -41.73, 0.4897];
e2d_mcmc_explore(100, 3.0/delta_chi);
 ====================
   ... etc ...
 ====================

% Make FOM-log plots with the cg fit and MCMC points together:

e2d_fomlog_plot("rbeta","beta", delta_chi);
()=system("cp fomlog_plot.ps fomlog_rbeta_beta_wMCMC_plot.ps");

e2d_fomlog_plot("ang","ratio",delta_chi);
()=system("cp fomlog_plot.ps fomlog_ang_ratio_wMCMC_plot.ps");

% Make a plot of the parameter values vs iteration:

PSwin=open_plot("fomlog_ps_vs_iter.ps/VCPS",1,4);
title(""); xlabel("Iteration");
variable str_tag = "rbeta";
ylabel(str_tag); yrange();
plot( [0:e2d_fomlog_i-1], array_struct_field(e2d_fomlog_ps[[0:e2d_fomlog_i-1]], str_tag));
variable str_tag = "ratio";
ylabel(str_tag); yrange();
plot( [0:e2d_fomlog_i-1], array_struct_field(e2d_fomlog_ps[[0:e2d_fomlog_i-1]], str_tag));
variable str_tag = "ang";
ylabel(str_tag); yrange();
plot( [0:e2d_fomlog_i-1], array_struct_field(e2d_fomlog_ps[[0:e2d_fomlog_i-1]], str_tag));
variable str_tag = "beta";
ylabel(str_tag); yrange();
plot( [0:e2d_fomlog_i-1], array_struct_field(e2d_fomlog_ps[[0:e2d_fomlog_i-1]], str_tag));
close_plot(PSwin);

% The rbeta-beta plot makes it clear that rbeta and beta are very
correlated. The parameters vs iter plot show that rbeta is
heading toward smaller values. So, let's fix rbeta and do MCMC...

Reset the FOM-log to iter 196, just after the cg fitting:

e2d_fomlog_i=196;

% And do MCMC with rbeta fixed:

public define e2d_update_pars () 
{ 
   % this line is the same for all cg fits: 
   variable pvals = cg_p0 + cg_smand*cg_dir*e2d_scan_vals[e2d_iscan]; 
   % map the fitting parameter values to the desired model params: 
   s3d_ps.ratio = pvals[0]; 
   s3d_ps.ang = pvals[1]; 
   s3d_ps.beta = pvals[2]; 
}

% fix it at the cg fit value:

s3d_ps.rbeta=7.832;

% and set the other three parameters:

cg_p0 = [1.20, -30.0, 0.50];
cg_smand = [0.2, 15.0, 0.04];
variable delta_chi = sqrt(2.0*e2d_data[1].nfom);

% Now do MCMC explore with just these 3 parameters:

e2d_mcmc_explore(400, 1.0/delta_chi);
 ====================
   ... etc ...
 ====================
e2d_fomlog_plot("ang","ratio",delta_chi);
()=system("cp fomlog_plot.ps fomlog_ang_ratio_wMCMCrbfixed_plot.ps");

e2d_fomlog_plot("ang","beta", delta_chi);
()=system("cp fomlog_plot.ps fomlog_ang_beta_wMCMCrbfixed_plot.ps");

% Remake the ps vs iter plots for this case too (using code above and different file name for PSwin)

Abel 1068 Event-2D : Fitting and MCMC

Abel 1068 Event-2D :
Fitting and MCMC