X 1822 Event-2D :
C-G and MCMC in 1D


Back to Event-2D Demos page
 10/9/08
dd@space.mit.edu

This page shows the effects of the commands in x1822_cgmc_play.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. 


%   Apply the event2d cg fitting and mcmcm exploring to
%   the x1822 spectral modeling as an example.

%   Read in the x1822 data by doing:

.source x1822_cgmc_data

Current Spectrum List:
 id    instrument   m prt src    use/nbins   A   R     totcts   exp(ksec)  target
  1     HETG-ACIS  -1  2   1     502/ 1024   1   1  1.1262e+05    80.250  X1822-371
  2     HETG-ACIS   1  1   1     361/ 1024   2   2  8.0612e+04    80.250  X1822-371
       ...
Parameters[Variable] = 9[5]
            Data bins = 863
           Chi-square = 842.8604    (sigma = gehrels)
   Reduced chi-square = 0.9823548
phabs(1)*(cutoffpl(1)+bbody(1)+gauss(1))
 idx  param             tie-to  freeze         value         min         max
  1  phabs(1).nH            0     0       0.05956062           0      100000  10^22
  2  cutoffpl(1).norm       0     0       0.01538339           0       1e+10  
  3  cutoffpl(1).PhoIndex   0     0       -0.5294097          -2           9  
  4  cutoffpl(1).HighECut   0     1             3.64           1         500  keV
  5  bbody(1).norm          0     0     0.0002600437           0        0.01  
  6  bbody(1).kT            0     0        0.2203471        0.05           3  keV
  7  gauss(1).area          0     1          0.00033           0        0.01  photons/s/cm^2
  8  gauss(1).center        0     1         1.939947           0           0  A
  9  gauss(1).sigma         0     1      0.007719233           0           1  A

Here's a plot of the spectral data loaded into ISIS with selections and binning set as desired.
The best-fit model is also over plotted for each data set (MEG m=-1: black/red, HEG m=+1: blue/purple). 

% - - -
%   Make standard 2-parameter contour plots
%   

   thaw([1,3,6]);
   
   ()=fit_counts;
   px = conf_grid("phabs(1).nH", 0.003,0.25, 15);
   py = conf_grid("bbody(1).kT", 0.15,0.30, 15);
   cmap = conf_map_counts (px, py);
   save_conf(cmap, "cmap_nH_kT_v1.fits");

   xlin,xrange();
   ylin,yrange();
   xlabel(cmap.px.index); ylabel(cmap.py.index);
   plot_conf(cmap);
   

   ()=fit_counts;
   px = conf_grid("phabs(1).nH", 0.003,0.25, 15);
   py = conf_grid("cutoffpl(1).PhoIndex", -0.70,-0.40, 15);
   cmap = conf_map_counts (px, py);
   save_conf(cmap, "cmap_nH_PhoIndex_v1.fits");

   xlin,xrange();
   ylin,yrange();
   xlabel(cmap.px.index); ylabel(cmap.py.index);
   plot_conf(cmap);

Using standard ISIS commands, contour plots of the confidence intervals for pairs of parameters (nH & kT, nH & PhoIndex) are generated.
The plots below were re-created from the saved fits files and set to the desired x and y range values for comparison to C-G and MCMC output.



The C-G fitting and MCMC exploring were carried out (code below) and these plots show the FOM-log outputs in parameter space.
The scales are the same as the confidence plot above each of them.

The highlighted MCMC points (red "o"s) have chi-squared values within 2.7 of the minumum
and they do seem to roughly agree with the ranges of the contours above.



The C-G fitting trajectory is shown by the light blue and dark blue colors. The black "X" on yellow marks the C-G best fit result
(click plot for larger version.)
The MCMC exploration points are in red (accepted steps) and yellow (rejected trial points).
The symbols are changed from "/" to "o" when the chi-squared at the point is below a level.
The plot below shows these chi-squared ("FOM") values as a function of iteration in the fitting and MCMC exploration. 

% - - -
% 
%   Setup to do CG fitting and MCMC exploration
%   of this model-data over the three interesting parameters.
%   (The no-argument help message for e2d_cg_fit reminds of these steps.)
% i)
public define e2d_norm_optimum (toleran) { return 1.0; }
% ii)
s3d_ps = struct {nH, kT, PhoIndex };
% iii)
public define e2d_pars2fom ()
{
   e2d_update_pars;  % used for scanning and fitting
   %
   % Set the three interesting parameters from s3d_ps
   % and then do a fit with the 2 norms and return chi2:
   set_par("phabs(1).nH",s3d_ps.nH,1);
   set_par("bbody(1).kT",s3d_ps.kT,1);
   set_par("cutoffpl(1).PhoIndex",s3d_ps.PhoIndex,1);
   variable fit_out;
   ()=fit_counts(Assigned_RMFARF, &fit_out);
   % Set the chi-squared into the e2d variable where it is expected,
   % also a little noise to it since it is expected by cg fit:
   e2d_tot_fom = fit_out.statistic + 0.1*grand();
}
% iv)  For cg_fit and mcmc_explore:
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.nH = pvals[0]; 
   if (pvals[0] < 0.0) s3d_ps.nH = 0.0;
   s3d_ps.kT = pvals[1]; 
   s3d_ps.PhoIndex = pvals[2];
} 

% Starting point and smand values:
cg_p0 = [0.07, 0.20, -0.50];
cg_smand = [0.02, 0.01, 0.015];


% Do CG-fit
cg_p0 = [0.07, 0.20, -0.50];
e2d_cg_fit(1);


% Do MCMC
cg_p0 = [0.07, 0.20, -0.50];
e2d_mcmc_explore(300);


% Look at FOM-log info:
e2d_fomlog_plot("nH","kT");
e2d_fomlog_plot("nH","PhoIndex");

% - - - - -