5. ADI forward modeling of point sources

Authors: Valentin Christiaens and Carlos Alberto Gomez Gonzalez
Suitable for VIP v1.0.0 onwards
Last update: 2022/09/21

Table of contents

This tutorial shows:

  • how to generate and inject fake companions in a cube;
  • how to estimate the astrometry and photometry of a directly imaged companion, and associated uncertainties.

Let’s first import a couple of external packages needed in this tutorial:

[1]:
%matplotlib inline
from hciplot import plot_frames, plot_cubes
from matplotlib.pyplot import *
from matplotlib import pyplot as plt
import numpy as np
from packaging import version

In the following box we import all the VIP routines that will be used in this tutorial. The path to some routines has changed between versions 1.0.3 and 1.1.0, which saw a major revamp of the modular architecture, hence the if statements.

[2]:
import vip_hci as vip
vvip = vip.__version__
print("VIP version: ", vvip)
if version.parse(vvip) < version.parse("1.0.0"):
    msg = "Please upgrade your version of VIP"
    msg+= "It should be 1.0.0 or above to run this notebook."
    raise ValueError(msg)
elif version.parse(vvip) <= version.parse("1.0.3"):
    from vip_hci.conf import VLT_NACO
    from vip_hci.metrics import cube_inject_companions
    from vip_hci.negfc import (confidence, firstguess, mcmc_negfc_sampling, show_corner_plot, show_walk_plot,
                               speckle_noise_uncertainty)
    from vip_hci.pca import pca, pca_annular, pca_annulus, pca_grid
else:
    from vip_hci.config import VLT_NACO
    from vip_hci.fm import (confidence, cube_inject_companions, cube_planet_free, firstguess, mcmc_negfc_sampling,
                            normalize_psf, show_corner_plot, show_walk_plot, speckle_noise_uncertainty)
    from vip_hci.psfsub import median_sub, pca, pca_annular, pca_annulus, pca_grid

# common to all versions:
from vip_hci.fits import open_fits, write_fits, info_fits
from vip_hci.metrics import contrast_curve, detection, significance, snr, snrmap, throughput
from vip_hci.preproc import frame_crop
from vip_hci.var import fit_2dgaussian, frame_center
VIP version:  1.5.0

5.1. Loading ADI data

In the ‘dataset’ folder of the VIP_extras repository you can find a toy ADI (Angular Differential Imaging) cube and a NACO point spread function (PSF) to demonstrate the capabilities of VIP. This is an L’-band VLT/NACO dataset of beta Pictoris published in Absil et al. (2013) obtained using the Annular Groove Phase Mask (AGPM) Vortex coronagraph. The sequence has been heavily sub-sampled temporarily to make it smaller. The frames were also cropped to the central 101x101 area. In case you want to plug-in your cube just change the path of the following cells.

More info on this dataset, and on opening and visualizing fits files with VIP in general, is available in Tutorial 1. Quick start.

Let’s load the data:

[3]:
psfnaco = '../datasets/naco_betapic_psf.fits'
cubename = '../datasets/naco_betapic_cube_cen.fits'
angname = '../datasets/naco_betapic_pa.fits'

cube = open_fits(cubename)
psf = open_fits(psfnaco)
pa = open_fits(angname)
Fits HDU-0 data successfully loaded, header available. Data shape: (61, 101, 101)
Fits HDU-0 data successfully loaded, header available. Data shape: (39, 39)
Fits HDU-0 data successfully loaded, header available. Data shape: (61,)
[4]:
derot_off = 104.84 # NACO derotator offset for this observation (Absil et al. 2013)
TN = -0.45         # Position angle of true north for NACO at the epoch of observation (Absil et al. 2013)

angs = pa+derot_off+TN

Let’s measure the FWHM by fitting a 2D Gaussian to the core of the unsaturated non-coronagraphic PSF:

[5]:
%matplotlib inline
DF_fit = fit_2dgaussian(psf, crop=True, cropsize=9, debug=True, full_output=True)
../_images/tutorials_05_fm_planets_14_0.png
FWHM_y = 4.733218722257211
FWHM_x = 4.4736824050602895

centroid y = 19.006680059041177
centroid x = 18.999424475165444
centroid y subim = 4.006680059041176
centroid x subim = 3.9994244751654446

amplitude = 0.10413004853269539
theta = -34.08563676834199
[6]:
fwhm_naco = np.mean([DF_fit['fwhm_x'],DF_fit['fwhm_y']])
print(fwhm_naco)
4.60345056365875

Let’s normalize the flux to one in a 1xFWHM aperture and crop the PSF array:

[7]:
psfn = normalize_psf(psf, fwhm_naco, size=19, imlib='ndimage-fourier')
Flux in 1xFWHM aperture: 1.228
[8]:
plot_frames(psfn, grid=True, size_factor=4)
../_images/tutorials_05_fm_planets_18_0.png

Let’s finally define the pixel scale for NACO (L’ band), which we get from a dictionary stored in the config subpackage:

[9]:
pxscale_naco = VLT_NACO['plsc']
print(pxscale_naco, "arcsec/px")
0.02719 arcsec/px

5.2. Generating and injecting synthetic planets

We first select an image library imlib for image operations (shifts, rotations) and associated interpolation order. ‘vip-fft’ is more accurate, but ‘skimage’ is faster, and ‘opencv’ even faster - see Tutorial 7 for more details.

[10]:
imlib_rot = 'opencv' #'skimage'
interpolation='lanczos4' #'biquintic'

The cube_inject_companions function in the fm module (VIP versions >= 1.1.0) makes the injection of fake companions at arbitrary fluxes and locations very easy. The normalized non-coronagraphic PSF should be provided for the injection. If the user does not have access to an observed PSF, the create_synth_psf from the var module can be used to create synthetic ones (based on 2D Gaussian, Moffat or Airy models).

Some procedures, e.g. the negative fake companion technique and the contrast curve generation, heavily rely on the injection of fake companions. The coordinates for the injection should be provided in the derotated image, while the actual injection occurs in the images of the input cube, i.e. in the rotated field.

[11]:
rad_fc = 30.5
theta_fc = 240
flux_fc = 400.

gt = [rad_fc, theta_fc, flux_fc]
[12]:
cubefc = cube_inject_companions(cube, psf_template=psfn, angle_list=angs, flevel=flux_fc, plsc=pxscale_naco,
                                rad_dists=[rad_fc], theta=theta_fc, n_branches=1,
                                imlib=imlib_rot, interpolation=interpolation)

Let’s set the corresponding cartesian coordinates:

[13]:
cy, cx = frame_center(cube[0])
x_fc = cx + rad_fc*np.cos(np.deg2rad(theta_fc))
y_fc = cy + rad_fc*np.sin(np.deg2rad(theta_fc))

xy_test = (x_fc, y_fc)
print('({:.1f}, {:.1f})'.format(xy_test[0],xy_test[1]))
(34.7, 23.6)

Let’s double-check the fake companion was injected at the right location, by post-processing the cube and checking the final image. Let’s use PCA, and infer the optimal \(n_{\rm pc}\) while we are at it - this will be useful for the next section.

[14]:
res_ann_opt = pca_grid(cubefc, angs, fwhm=fwhm_naco, range_pcs=(1,41,1), source_xy=xy_test, mode='annular',
                       annulus_width=4*fwhm_naco, imlib=imlib_rot, interpolation=interpolation,
                       full_output=True, plot=True)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2023-08-01 01:21:44
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.258349
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 41
Optimal number of PCs = 10, for S/N=13.415
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 34.7, 23.6
Flux in a centered 1xFWHM circular aperture = 104.785
Central pixel S/N = 18.264
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 13.415
Max S/N (shifting the aperture center) = 18.472
stddev S/N (shifting the aperture center) = 4.088

../_images/tutorials_05_fm_planets_30_1.png

The grid search looking for the optimal number of principal components (npc) found that 10 principal components maximizes the S/N ratio of the injected fake companion.

[15]:
_, final_ann_opt, _, opt_npc_ann = res_ann_opt
[16]:
plot_frames(final_ann_opt, label='PCA-ADI annulus (opt. npc={:.0f})'.format(opt_npc_ann),
            dpi=100, vmin=-5, colorbar=True, circle=xy_test)
../_images/tutorials_05_fm_planets_33_0.png

We can see that the fake companion was indeed injected at the requested location.

5.3. Flux and position estimation with NEGFC

When a companion candidate is detected, the next step is to characterize it, i.e. infer its exact position (astrometry) and flux (photometry).

Question 5.1: Why would a simple 2D Gaussian fit (as performed e.g. for the stellar PSF in Section 5.1) be inappropriate to extract the astrometry and photometry of a candidate companion?

VIP implements the Negative fake companion (NEGFC) technique for robust extraction of the position and flux of detected point-like sources. The technique can be summarized as follow (see full description in Wertz et al. 2017):

  1. Estimate the position and flux of the planet, from either the visual inspection of reduced images or a previous estimator (see ABC below).
  2. Scale (in flux) and shift the normalized off-axis PSF to remove the estimate from the input data cube.
  3. Process the cube with PCA in a single annulus encompassing the point source.
  4. Measure residuals in an aperture centered on the approximate location of the companion candidate.
  5. Iterate on the position and flux of the injected negative PSF (steps 2-4), until the absolute residuals in the aperture are minimized (i.e. the injected negative companion flux and the position match exactly that of the true companion).

Iterations between steps 2-4 can be performed in one of 3 ways - sorted in increasing computation time and accuracy:

  1. a grid search on the flux only, provided a fixed estimate of the position (implemented in the firstguess function);
  2. a Nelder-Mead simplex algorithm (firstguess function with the simplex=True option);
  3. an MCMC sampler, which has the advantage to also yield uncertainties on each of the parameters of the point source (mcmc_negfc_sampling function).

Different figures of merit can be used for minimization of the residuals (Wertz et al. 2017; Christiaens et al. 2021):

\[\chi^2 = \sum_j^N \frac{(I_j-\mu)^2}{\sigma^2}~~({\rm default}), ~~~~~~~~\chi^2 = \sum_j^N |I_j|, ~~~~~~~~\chi^2 = N {\rm std}(I_j).\]

where \(j \in {1,...,N}\), \(N\) is the total number of pixels contained in the circular aperture around the companion candidate, \(\mu\) and \(\sigma\) are the mean and standard deviation (\(N_{\rm resel}\) degrees of freedom) of pixel intensities in a truncated annulus at the radius of the companion candidate, but avoiding the azimuthal region encompassing the negative side lobes.

5.3.1. Nelder-Mead based optimization

With the function firstguess, we can obtain a first estimation of the flux and position by running A) a naive grid minimization (grid of values for the flux through parameter f_range), and B) a Nelder-mead based minimization (if the parameter simplex is set to True). The latter is done based on the preliminary guess of the grid minimization. The maximum number of iterations and error can be set with the parameter simplex_options as a dicitionary (see scipy.minimize function for the Nelder-Mead options).

Fisrt we define the position of the sources by examining a flux frame or S/N map. planets_xy_coord takes a list or array of X,Y pairs like ((x1,y1),(x2,y2)…(x_n,y_n)). Let’s take the coordinates of the previously injected companion.

Let’s test the algorithm with different values for the # of PCs: 5 and 25.

[17]:
r_lo, theta_lo, f_lo = firstguess(cubefc, angs, psfn, ncomp=5, planets_xy_coord=[xy_test],
                                  fwhm=fwhm_naco, f_range=None, annulus_width=4*fwhm_naco,
                                  aperture_radius=2, simplex=True, imlib=imlib_rot,
                                  interpolation=interpolation, plot=True, verbose=True)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2023-08-01 01:21:49
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
             Planet 0
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Planet 0: flux estimation at the position [34.749999999999986,23.58622518457463], running ...
Step | flux    | chi2r
1/30   0.100   1.933
2/30   0.149   1.933
3/30   0.221   1.932
4/30   0.329   1.931
5/30   0.489   1.930
6/30   0.728   1.927
7/30   1.083   1.924
8/30   1.610   1.919
9/30   2.395   1.912
10/30   3.562   1.901
11/30   5.298   1.886
12/30   7.880   1.862
13/30   11.721   1.827
14/30   17.433   1.776
15/30   25.929   1.698
16/30   38.566   1.581
17/30   57.362   1.416
18/30   85.317   1.191
19/30   126.896   0.907
20/30   188.739   0.543
21/30   280.722   0.175
22/30   417.532   0.043
23/30   621.017   0.751
24/30   923.671   3.713
25/30   1373.824   11.487
../_images/tutorials_05_fm_planets_40_1.png
Planet 0: preliminary position guess: (r, theta)=(30.5, 240.0)
Planet 0: preliminary flux guess: 417.5
Planet 0: Simplex Nelder-Mead minimization, running ...
Planet 0: Success: True, nit: 88, nfev: 208, chi2r: 0.027626969432102817
message: Optimization terminated successfully.
Planet 0: simplex result: (r, theta, f)=(30.547, 239.971, 384.840) at
          (X,Y)=(34.71, 23.55)

 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
DONE !
 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Running time:  0:01:15.364625
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
[18]:
r_hi, theta_hi, f_hi = firstguess(cubefc, angs, psfn, ncomp=25, planets_xy_coord=[xy_test],
                                  fwhm=fwhm_naco, f_range=None, annulus_width=4*fwhm_naco,
                                  aperture_radius=2, imlib=imlib_rot, interpolation=interpolation,
                                  simplex=True, plot=True, verbose=True)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2023-08-01 01:23:05
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
             Planet 0
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Planet 0: flux estimation at the position [34.749999999999986,23.58622518457463], running ...
Step | flux    | chi2r
1/30   0.100   0.322
2/30   0.149   0.322
3/30   0.221   0.322
4/30   0.329   0.322
5/30   0.489   0.322
6/30   0.728   0.322
7/30   1.083   0.322
8/30   1.610   0.322
9/30   2.395   0.322
10/30   3.562   0.322
11/30   5.298   0.321
12/30   7.880   0.320
13/30   11.721   0.319
14/30   17.433   0.317
15/30   25.929   0.313
16/30   38.566   0.310
17/30   57.362   0.301
18/30   85.317   0.272
19/30   126.896   0.244
20/30   188.739   0.197
21/30   280.722   0.117
22/30   417.532   0.064
23/30   621.017   0.210
24/30   923.671   0.464
25/30   1373.824   0.559
../_images/tutorials_05_fm_planets_41_1.png
Planet 0: preliminary position guess: (r, theta)=(30.5, 240.0)
Planet 0: preliminary flux guess: 417.5
Planet 0: Simplex Nelder-Mead minimization, running ...
Planet 0: Success: True, nit: 58, nfev: 168, chi2r: 0.05751191932736462
message: Optimization terminated successfully.
Planet 0: simplex result: (r, theta, f)=(30.186, 239.897, 424.143) at
          (X,Y)=(34.86, 23.89)

 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
DONE !
 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Running time:  0:00:49.105108
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

For both the \(n_{\rm pc} = 5\) and \(n_{\rm pc} = 25\) cases, the parameters estimated by the Nelder-Mead optimization are not exactly equal to the original values (radius=30.5, theta=240, flux=400), which reflects:

  • the limitations of this heuristic minization procedure (depending on the initial guess the minimization may get trapped in a different local minimum);
  • the higher residual speckle noise level in images obtained with low \(n_{\rm pc}\) values;
  • the higher self-subtraction for high \(n_{\rm pc}\) values.

These estimates are provided without uncertainties (error bars). We will come back to this question later on.

For comparison, let’s use the optimal \(n_{\rm pc} = 10\) found in Section 5.2:

[19]:
r_0, theta_0, f_0 = firstguess(cubefc, angs, psfn, ncomp=opt_npc_ann, planets_xy_coord=[xy_test],
                               fwhm=fwhm_naco, f_range=None, annulus_width=4*fwhm_naco,
                               aperture_radius=2, imlib=imlib_rot, interpolation=interpolation,
                               simplex=True, plot=True, verbose=True)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2023-08-01 01:23:54
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
             Planet 0
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Planet 0: flux estimation at the position [34.749999999999986,23.58622518457463], running ...
Step | flux    | chi2r
1/30   0.100   1.951
2/30   0.149   1.951
3/30   0.221   1.950
4/30   0.329   1.950
5/30   0.489   1.948
6/30   0.728   1.946
7/30   1.083   1.944
8/30   1.610   1.939
9/30   2.395   1.933
10/30   3.562   1.924
11/30   5.298   1.911
12/30   7.880   1.891
13/30   11.721   1.861
14/30   17.433   1.818
15/30   25.929   1.749
16/30   38.566   1.661
17/30   57.362   1.523
18/30   85.317   1.288
19/30   126.896   1.010
20/30   188.739   0.637
21/30   280.722   0.218
22/30   417.532   0.052
23/30   621.017   0.634
24/30   923.671   2.632
25/30   1373.824   6.486
../_images/tutorials_05_fm_planets_43_1.png
Planet 0: preliminary position guess: (r, theta)=(30.5, 240.0)
Planet 0: preliminary flux guess: 417.5
Planet 0: Simplex Nelder-Mead minimization, running ...
Planet 0: Success: True, nit: 68, nfev: 165, chi2r: 0.04395945745570059
message: Optimization terminated successfully.
Planet 0: simplex result: (r, theta, f)=(30.534, 239.933, 390.873) at
          (X,Y)=(34.70, 23.57)

 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
DONE !
 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Running time:  0:00:26.698135
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

We see that using the optimal \(n_{\rm pc}\) leads to a closer parameter estimates to the ground truth.

Answer 5.1: If relying on a 2D Gaussian fit, both the photometry and astrometry would be biased by self-subtraction and the negative side lobes.

5.3.2. Planet subtraction

Let’s use the values obtained with the simplex optimization to subtract the planet with the function cube_planet_free.

First we define a list with the parameters (r, theta, flux) is each companion that we obtained via the NegFC, in this case one:

[20]:
plpar_fc = [(r_0[0], theta_0[0], f_0[0])]

Note: r_0, theta_0 and f_0 have the same length as the number of planet coordinates planets_xy_coord provided to firstguess. Here there is only one planet, so we take the zeroth index. The number of tuples in (i.e the length of) plpar_fc should match the number of planets.

[21]:
cube_emp = cube_planet_free(plpar_fc, cubefc, angs, psfn, imlib=imlib_rot, interpolation=interpolation)

Let’s double-check the fake companion was well removed by computing a PCA post-processed image:

[22]:
from vip_hci.preproc import cube_derotate
from vip_hci.config import time_ini, timing
t0 = time_ini()
for i in range(100):
    fr_pca_emp = pca_annulus(cube_emp, angs, ncomp=opt_npc_ann, annulus_width=4*fwhm_naco,
                             r_guess=rad_fc, imlib=imlib_rot, interpolation=interpolation)
    #cube_tmp = cube_derotate(cube_emp, angs, imlib=imlib_rot, interpolation=interpolation, edge_blend='')
timing(t0)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2023-08-01 01:24:20
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Running time:  0:00:11.764804
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Let’s take a look at the PSF of the planet in the full-frame PCA final image and the same PSF in the frame resulting of processing the planet-subtracted cube:

[23]:
cropped_frame1 = frame_crop(final_ann_opt, cenxy=xy_test, size=15)
New shape: (15, 15)
[24]:
cropped_frame2 = frame_crop(fr_pca_emp, cenxy=xy_test, size=15)
New shape: (15, 15)

Let’s use both 'mode=surface' and the default image mode of plot_frames to show the residuals in the vicinity of the companion:

[25]:
plot_frames((cropped_frame1, cropped_frame2), mode='surface', vmax=8)
../_images/tutorials_05_fm_planets_58_0.png
[26]:
plot_frames((final_ann_opt, fr_pca_emp), vmin = float(np.amin(final_ann_opt)),
            vmax= float(np.amax(final_ann_opt)), axis=False)
../_images/tutorials_05_fm_planets_59_0.png

Not only the bright point-like signal is subtracted, but so are the negative side lobes. A subtraction not leaving any significant artifact/defect is a good sign that the inferred parameters are correct. However, keep in mind that even for slightly inaccurate parameters, the final image can still look relatively clean. Let’s take for example the parameters inferred with non-optimal \(n_{\rm pc}\):

[27]:
# planet parameters inferred from npc=5 used to make empty cube, then processed with PCA (npc=10)
plpar_fc = [(r_lo[0], theta_lo[0], f_lo[0])]
print(plpar_fc)
cube_emp = cube_planet_free(plpar_fc, cubefc, angs, psfn, imlib=imlib_rot, interpolation=interpolation)
final_ann_5 = pca_annulus(cubefc, angs, ncomp=5, annulus_width=4*fwhm_naco, r_guess=30.5,
                                  imlib=imlib_rot, interpolation=interpolation)
fr_pca_emp5 = pca_annulus(cube_emp, angs, ncomp=5, annulus_width=4*fwhm_naco, r_guess=30.5,
                                  imlib=imlib_rot, interpolation=interpolation)
plot_frames((final_ann_5, fr_pca_emp5), vmin = float(np.amin(final_ann_5)),
            vmax= float(np.amax(final_ann_5)), axis=False)
[(30.547288438104047, 239.9708851445477, 384.8400913238594)]
../_images/tutorials_05_fm_planets_61_1.png
[28]:
# parameters inferred from npc=25 used to make empty cube, then processed with PCA (npc=10)
plpar_fc = [(r_hi[0], theta_hi[0], f_hi[0])]
print(plpar_fc)
cube_emp = cube_planet_free(plpar_fc, cubefc, angs, psfn, imlib=imlib_rot, interpolation=interpolation)
final_ann_25 = pca_annulus(cubefc, angs, ncomp=25, annulus_width=4*fwhm_naco, r_guess=30.5,
                                   imlib=imlib_rot, interpolation=interpolation)
fr_pca_emp25 = pca_annulus(cube_emp, angs, ncomp=25, annulus_width=4*fwhm_naco, r_guess=30.5,
                                   imlib=imlib_rot, interpolation=interpolation)
plot_frames((final_ann_25, fr_pca_emp25), vmin = float(np.amin(final_ann_25)),
            vmax= float(np.amax(final_ann_25)), axis=False)
[(30.185654487910202, 239.89660151207698, 424.14272218644646)]
../_images/tutorials_05_fm_planets_62_1.png

Inaccurate parameters still leading to an apparently good removal of the companion brings the question of the uncertainties on each of the three parameters characterizing the companion. The next sections are dedicated to this question.

5.3.3. NEGFC technique coupled with MCMC

5.3.3.1. Running the MCMC sampler

MCMC is a more robust way of obtaining the flux and position. It samples the posterior distributions of the parameters and from them we can infer both the most likely parameter values and uncertainties on each parameter. The relevant function is mcmc_negfc_sampling, which can accept a number of parameters. Let’s define them in the next few boxes:

Let’s first define observation-related parameters, such as the non-coronagraphic psf, its FWHM and the pixel scale od the detector:

[29]:
obs_params = {'psfn': psfn,
              'fwhm': fwhm_naco}

In NEGFC, PCA in a single annulus is used by default to speed up the algorithm - although other algorithms can be used through the algo parameter. Let’s set the \(n_{\rm pc}\) to the optimal \(n_{\rm pc}\) found in Section 5.2. We set the width of the annulus on which PCA is performed (in pixels) with the annulus_width parameter. We also set a few other algorithm-related parameters in the following box:

[30]:
annulus_width = 4*fwhm_naco

algo_params = {'algo': pca_annulus,
               'ncomp': opt_npc_ann,
               'annulus_width': annulus_width,
               'svd_mode': 'lapack',
               'imlib': imlib_rot,
               'interpolation': interpolation}

The choice of log-likelihood expression to be used is determined by mu_sigma and fm. If the former is True (default; mu and sigma calculated automatically) or a tuple of 2 values, the following figure of merit will be used: \(\chi^2 = \sum_j^N \frac{(I_j-\mu)^2}{\sigma^2}\). Otherwise, the choice will be determined by fm: ‘sum’ for the sum of absolute residuals, or ‘stddev’ for the standard deviation of residuals (which can be useful if the point source is contained within a more extended signal). aperture_radius is the radius of the aperture (in fwhm units) in which the residual intensities \(I_j\) are considered.

[31]:
mu_sigma=True
aperture_radius=2

negfc_params = {'mu_sigma': mu_sigma,
                'aperture_radius': aperture_radius}

Parameter initial_state corresponds to the initial first estimation of the planets parameters (r, theta, flux). We set it to the result of the simplex optimization, obtained with optimal \(n_{\rm pc}\). Note that the MCMC minimization can only run for one companion candidate at a time, hence the first dimension of init should always be 3 (not the number of planets, as opposed to planets_xy_coord in firstguess).

[32]:
initial_state = np.array([r_0[0], theta_0[0], f_0[0]])

Beware that the MCMC procedure is a CPU intensive procedure and can take several hours when run properly on a standard laptop. We use the affine invariant sampler from emcee which can be run in parallel (nproc sets the number of CPUs to be used). At least 100 walkers are recommended for our MCMC chain, although both the number of walkers and iterations will depend on your dataset. For the sake of preventing this tutorial to take too long to run, we set the maximum number of iterations to 500, although feel free to significantly increase it in case of non-convergence.

[33]:
from multiprocessing import cpu_count

nwalkers, itermin, itermax = (100, 200, 500)

mcmc_params = {'nwalkers': nwalkers,
               'niteration_min': itermin,
               'niteration_limit': itermax,
               'bounds': None,
               'nproc': cpu_count()//2}

Another update from Christiaens et al. (2021) is that the convergence can now be evaluated based on the auto-correlation time (see more details in the documentation of ``emcee` <https://emcee.readthedocs.io/en/stable/tutorials/autocorr/>`__), instead of the Gelman-Rubin test, which is inappropriate for non-independent samples (such as an MCMC chain). We set the convergence test to the autocorrelation time based criterion using conv_test='ac' (instead of Gelman-Rubin ‘gb’). We also set the autocorrelation criterion \(N/\tau >= a_c\), where \(N\) is the number of samples and \(\tau\) the autocorrelation time, to ac_c=50 (the value recommended in the emcee documentation). Finally, we set the number of consecutive times the criterion must be met to: ac_count_thr=1, and the maximum gap in number of steps between 2 checks of the convergence criterion to: check_maxgap=50. If the maximum number of iterations niteration_limit is large enough, the chain will stop upon meeting the convergence criterion (spoiler: it won’t be the case here since we chose a small value for niteration_limit).

[34]:
conv_test, ac_c, ac_count_thr, check_maxgap = ('ac', 50, 1, 50)

conv_params = {'conv_test': conv_test,
               'ac_c': ac_c,
               'ac_count_thr': ac_count_thr,
               'check_maxgap': check_maxgap}

Setting bounds=None does not mean no bounds are used for parameter exploration, but rather that they are set automatically to:

  • \(r \in [r_0-w_{ann}/2, r_0+w_{ann}/2]\), where \(w_{ann}\) is the annulus_width,
  • \(\theta \in [\theta_0-\Delta {\rm rot}, \theta_0+\Delta {\rm rot}]\), where \(\Delta {\rm rot}\) is the angle subtended by min(aperture_radius/2,``fwhm``) at \(r_0\),
  • \(f \in [0.1*f_0, 2*f_0]\),

where (\(r_0, \theta_0, f_0\)) = initial_state. If the bounds are provided manually (as a tuple of tuples), they will supersede the automatic setting above.

Now let’s start the sampler. Note that this step is computer intensive and may take a long time to run depending on your machine. Feel free to skip the next box if you do not wish to run the MCMC or can’t run it in a reasonable time. The results are already saved as ‘MCMC_results’ in the ‘datasets’ folder and will be loaded in the subsequent boxes.

[35]:
chain = mcmc_negfc_sampling(cubefc, angs, **obs_params, **algo_params, **negfc_params,
                            initial_state=initial_state, **mcmc_params, **conv_params,
                            display=True, verbosity=2, save=False, output_dir='./')
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2023-08-01 01:24:33
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
        MCMC sampler for the NEGFC technique
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
The mean and stddev in the annulus at the radius of the companion (excluding the PA area directly adjacent to it) are -0.05 and 1.50 respectively.
Beginning emcee Ensemble sampler...
emcee Ensemble sampler successful

Start of the MCMC run ...
Step  |  Duration/step (sec)  |  Remaining Estimated Time (sec)
0               3.34020                 1666.76205
1               2.19021                 1090.72508
2               2.19359                 1090.21224
3               2.31976                 1150.59898
4               2.69719                 1335.11152
5               2.51392                 1241.87450
6               2.24456                 1106.56561
7               2.41634                 1188.83928
8               2.23449                 1097.13704
9               2.26186                 1108.31140
10              2.21542                 1083.33940
11              2.50433                 1222.11304
12              2.53846                 1236.23051
13              2.31661                 1125.87246
14              2.24655                 1089.57917
15              2.26821                 1097.81364
16              2.14571                 1036.37890
17              2.13355                 1028.37158
18              2.22640                 1070.89936
19              2.12646                 1020.70032
20              2.16859                 1038.75653
21              2.65686                 1269.98099
22              2.41074                 1149.92107
23              2.37324                 1129.66224
24              2.14155                 1017.23720
25              2.14439                 1016.44276
26              2.14505                 1014.60818
27              2.11589                 998.69866
28              2.27338                 1070.76057
29              2.24540                 1055.33706
30              2.25070                 1055.58065
31              2.17123                 1016.13470
32              2.22591                 1039.49857
33              2.14495                 999.54577
34              2.34227                 1089.15741
35              2.25026                 1044.12203
36              2.13896                 990.33755
37              2.13012                 984.11405
38              2.13743                 985.35477
39              2.13073                 980.13534
40              2.14593                 984.98141
41              2.13185                 976.38867
42              2.14525                 980.37879
43              2.22047                 1012.53478
44              2.31905                 1055.16957
45              2.19298                 995.61292
46              2.37515                 1075.94521
47              2.21718                 1002.16581
48              2.13358                 962.24503
49              2.18660                 983.96955
50              2.08625                 936.72535

 ac convergence test in progress...
../_images/tutorials_05_fm_planets_81_1.png
51              2.54530                 1140.29306
52              2.08067                 930.05815
53              2.11232                 942.09561
54              2.09718                 933.24688
55              2.17275                 964.70011
56              2.10914                 934.34681
57              2.15011                 950.34685
58              2.17226                 957.96886
59              2.18168                 959.94008
60              2.11695                 929.34149
61              2.46803                 1080.99626
62              2.35530                 1029.26523
63              2.09239                 912.28030
64              2.23483                 972.15236
65              2.36501                 1026.41304
66              2.59780                 1124.84567
67              2.21290                 955.97366
68              2.15923                 930.62942
69              2.15128                 925.04825
70              2.18316                 936.57693
71              2.12676                 910.25114
72              2.11734                 904.10375
73              2.07810                 885.26975
74              2.16641                 920.72382
75              2.11242                 895.66820
76              2.13628                 903.64686
77              2.16101                 911.94749
78              2.15495                 907.23395
79              2.08001                 873.60420
80              2.12340                 889.70669
81              2.09617                 876.19906
82              2.21262                 922.66212
83              2.14954                 894.20989
84              2.15548                 894.52379
85              2.23857                 926.76674
86              2.08401                 860.69613
87              2.08148                 857.57017
88              2.07639                 853.39752
89              2.14784                 880.61604
90              2.13956                 875.08045
91              2.09796                 855.96890
92              2.06908                 842.11353
93              2.10363                 854.07297
94              2.06909                 837.98023
95              2.08402                 841.94327
96              2.13253                 859.40919
97              2.06121                 828.60441
98              2.05265                 823.11265
99              2.07235                 828.93960
100             2.05734                 820.87706

 ac convergence test in progress...
../_images/tutorials_05_fm_planets_81_3.png
101             2.69116                 1071.08128
102             2.05989                 817.77633
103             2.07748                 822.68406
104             2.23612                 883.26701
105             2.10922                 831.03110
106             2.06535                 811.68412
107             2.05835                 806.87438
108             2.06989                 809.32855
109             2.04803                 798.73053
110             2.08677                 811.75159
111             2.05931                 799.01267
112             2.10048                 812.88692
113             2.07984                 802.81785
114             2.06703                 795.80616
115             2.04899                 786.81370
116             2.06094                 789.33964
117             2.04725                 782.04874
118             2.04670                 779.79156
119             2.06307                 783.96546
120             2.06342                 782.03580
121             2.04304                 772.27063
122             2.06675                 779.16287
123             2.05629                 773.16654
124             2.07103                 776.63475
125             2.05197                 767.43678
126             2.05511                 766.55715
127             2.06125                 766.78649
128             2.08546                 773.70492
129             2.16224                 800.03065
130             2.12664                 784.73090
131             2.07121                 762.20381
132             2.09308                 768.15999
133             2.14571                 785.33023
134             2.10268                 767.47674
135             2.14343                 780.20961
136             2.17731                 790.36353
137             2.14328                 775.86808
138             2.06101                 744.02641
139             2.15655                 776.35980
140             2.08707                 749.25921
141             2.18048                 780.61041
142             2.17616                 776.88912
143             2.10955                 750.99838
144             2.08361                 739.68013
145             2.09245                 740.72695
146             2.09421                 739.25507
147             2.14371                 754.58486
148             2.06617                 725.22497
149             2.11631                 740.70850
150             2.10190                 733.56310

 ac convergence test in progress...
../_images/tutorials_05_fm_planets_81_5.png
151             2.61897                 911.39982
152             2.37419                 823.84462
153             2.83107                 979.55160
154             2.17449                 750.19974
155             2.09674                 721.27925
156             2.19557                 753.08154
157             2.21266                 756.72938
158             2.06697                 704.83609
159             2.08820                 709.98970
160             2.23117                 756.36697
161             2.50535                 846.80661
162             2.13338                 718.94872
163             2.08458                 700.41888
164             2.09756                 702.68427
165             2.07810                 694.08607
166             2.06960                 689.17780
167             2.32841                 773.03046
168             2.54923                 843.79678
169             2.80098                 924.32307
170             2.76579                 909.94392
171             2.76688                 907.53697
172             2.83906                 928.37229
173             2.90796                 947.99366
174             2.88044                 936.14332
175             2.81847                 913.18493
176             2.86848                 926.51775
177             3.20779                 1032.90709
178             2.80812                 901.40524
179             2.51933                 806.18624
180             2.78288                 887.73936
181             2.83167                 900.47138
182             2.81588                 892.63491
183             2.94505                 930.63580
184             2.90089                 913.77940
185             2.87036                 901.29430
186             2.81578                 881.34070
187             2.76816                 863.66467
188             3.04537                 947.10883
189             2.89733                 898.17199
190             2.89537                 894.66840
191             2.92775                 901.74700
192             2.89785                 889.63995
193             2.82564                 864.64431
194             2.84148                 866.65140
195             2.81781                 856.61454
196             2.90158                 879.17995
197             2.84025                 857.75671
198             2.84084                 855.09224
199             2.82461                 847.38240
200             2.84134                 849.56096

 ac convergence test in progress...
../_images/tutorials_05_fm_planets_81_7.png
Auto-corr tau/N = [0.07087539 0.06974356 0.06146271]
tau/N <= 0.02 = [False False False]

201             3.52743                 1051.17414
202             2.73947                 813.62318
203             2.63211                 779.10545
204             2.64475                 780.20007
205             2.79837                 822.72196
206             2.70746                 793.28695
207             2.65428                 775.05093
208             2.75838                 802.68829
209             2.84864                 826.10676
210             2.53620                 732.96296
211             2.60016                 748.84522
212             2.48939                 714.45378
213             2.63463                 753.50275
214             2.72670                 777.11007
215             2.70954                 769.50822
216             2.87432                 813.43228
217             2.77169                 781.61545
218             2.58293                 725.80249
219             2.64148                 739.61552
220             2.64630                 738.31686
221             2.74591                 763.36242
222             2.65482                 735.38514
223             2.73050                 753.61690
224             2.91694                 802.15740
225             2.72512                 746.68425
226             2.74792                 750.18298
227             2.63281                 716.12568
228             2.71890                 736.82109
229             2.74596                 741.40866
230             2.62960                 707.36186
231             2.85010                 763.82707
232             2.97704                 794.86861
233             2.95865                 787.00090
234             2.89675                 767.63875
235             2.89458                 764.16833
236             2.98838                 785.94315
237             2.94704                 772.12343
238             2.87541                 750.48331
239             2.89778                 753.42306
240             3.06505                 793.84717
241             2.85928                 737.69450
242             2.89128                 743.05845
243             2.98253                 763.52717
244             2.87274                 732.54895
245             2.80507                 712.48753
246             2.83471                 717.18188
247             2.90946                 733.18468
248             2.98783                 749.94633
249             2.92854                 732.13500
250             2.87079                 714.82696

 ac convergence test in progress...
../_images/tutorials_05_fm_planets_81_9.png
Auto-corr tau/N = [0.06574988 0.06474131 0.06063177]
tau/N <= 0.02 = [False False False]

251             3.47166                 860.97094
252             3.09527                 764.53095
253             3.02699                 744.64003
254             2.99275                 733.22277
255             2.91808                 712.01201
256             2.90441                 705.77114
257             2.81173                 680.43914
258             3.02314                 728.57794
259             2.88697                 692.87232
260             2.92427                 698.90053
261             2.89583                 689.20873
262             2.78900                 660.99395
263             2.82519                 666.74508
264             2.63651                 619.57867
265             2.58678                 605.30722
266             2.53243                 590.05666
267             2.70067                 626.55474
268             2.72029                 628.38745
269             2.73645                 629.38327
270             2.70971                 620.52359
271             2.73447                 623.45848
272             2.52986                 574.27731
273             2.71623                 613.86708
274             2.46491                 554.60587
275             2.50596                 561.33482
276             2.40330                 535.93657
277             2.54190                 564.30136
278             2.54490                 562.42246
279             2.66556                 586.42386
280             2.77727                 608.22213
281             2.75815                 601.27692
282             2.62267                 569.11982
283             2.74276                 592.43638
284             2.62823                 565.06924
285             2.81274                 601.92636
286             2.76452                 588.84233
287             2.76717                 586.63962
288             2.71862                 573.62840
289             2.71747                 570.66954
290             2.64078                 551.92281
291             2.69700                 560.97683
292             2.73664                 566.48365
293             2.74256                 564.96736
294             2.63584                 540.34638
295             2.53768                 517.68672
296             2.67744                 543.51931
297             2.72299                 550.04479
298             2.54489                 511.52289
299             2.59405                 518.81060
300             2.75976                 549.19244

 ac convergence test in progress...
../_images/tutorials_05_fm_planets_81_11.png
Auto-corr tau/N = [0.05992968 0.05819945 0.05582293]
tau/N <= 0.02 = [False False False]

301             3.67342                 727.33756
302             2.76478                 544.66068
303             2.69858                 528.92109
304             2.61604                 510.12799
305             2.52045                 488.96769
306             2.52258                 486.85736
307             2.50870                 481.67040
308             2.53884                 484.91806
309             2.66045                 505.48455
310             2.53255                 478.65157
311             2.76081                 519.03266
312             2.76065                 516.24211
313             2.76599                 514.47340
314             2.84052                 525.49694
315             2.74399                 504.89490
316             2.76381                 505.77650
317             2.84358                 517.53083
318             3.05685                 553.28985
319             2.81842                 507.31488
320             2.83538                 507.53230
321             2.76428                 492.04237
322             2.79499                 494.71270
323             2.68584                 472.70749
324             2.76566                 483.99050
325             2.77391                 482.65947
326             2.72442                 471.32535
327             2.73903                 471.11402
328             2.90955                 497.53322
329             2.76018                 469.23111
330             2.72374                 460.31189
331             2.85692                 479.96239
332             2.71598                 453.56816
333             2.79313                 463.65991
334             2.70836                 446.87874
335             2.75596                 451.97810
336             2.76465                 450.63876
337             2.73433                 442.96114
338             2.69640                 434.12008
339             2.84754                 455.60704
340             2.79521                 444.43887
341             2.59578                 410.13371
342             2.56583                 402.83578
343             2.54488                 397.00050
344             2.49920                 387.37678
345             2.47097                 380.52861
346             2.53352                 387.62887
347             2.46429                 374.57269
348             2.50170                 377.75625
349             2.71581                 407.37120
350             2.67446                 398.49424

 ac convergence test in progress...
../_images/tutorials_05_fm_planets_81_13.png
Auto-corr tau/N = [0.05508451 0.05191419 0.05251266]
tau/N <= 0.02 = [False False False]

351             3.43288                 508.06624
352             2.73257                 401.68764
353             2.58356                 377.20034
354             2.51574                 364.78172
355             2.51717                 362.47219
356             2.62229                 374.98761
357             2.40624                 341.68594
358             2.54885                 359.38757
359             2.54833                 356.76550
360             2.50697                 348.46911
361             2.56279                 353.66530
362             2.72524                 373.35829
363             2.65975                 361.72614
364             2.75432                 371.83333
365             2.58017                 345.74265
366             2.60259                 346.14420
367             2.42601                 320.23319
368             2.45330                 321.38243
369             2.40520                 312.67561
370             2.41590                 311.65149
371             2.41999                 309.75923
372             2.38995                 303.52352
373             2.43677                 307.03315
374             2.40730                 300.91275
375             2.56228                 317.72297
376             2.59147                 318.75118
377             2.57029                 313.57575
378             2.51667                 304.51755
379             2.65217                 318.26016
380             2.54162                 302.45254
381             2.64226                 311.78633
382             2.47525                 289.60390
383             2.42360                 281.13748
384             2.41674                 277.92487
385             2.43464                 277.54885
386             2.68985                 303.95350
387             2.53620                 284.05418
388             2.45539                 272.54829
389             2.48156                 272.97215
390             2.46928                 269.15119
391             2.43104                 262.55178
392             2.41755                 258.67753
393             2.38274                 252.57076
394             2.36740                 248.57689
395             2.43024                 252.74527
396             2.34727                 241.76850
397             2.24049                 228.52978
398             2.23457                 225.69117
399             2.22771                 222.77110
400             2.24155                 221.91305

 ac convergence test in progress...
../_images/tutorials_05_fm_planets_81_15.png
Auto-corr tau/N = [0.05289581 0.04833608 0.04925793]
tau/N <= 0.02 = [False False False]

401             2.82485                 276.83530
402             2.25909                 219.13183
403             2.24930                 215.93309
404             2.24570                 213.34121
405             2.24639                 211.16085
406             2.25156                 209.39480
407             2.26234                 208.13574
408             2.24111                 203.94128
409             2.25257                 202.73112
410             2.24761                 200.03711
411             2.31194                 203.45028
412             2.29979                 200.08156
413             2.35230                 202.29754
414             2.31974                 197.17756
415             2.28121                 191.62181
416             2.21409                 183.76930
417             2.24480                 184.07335
418             2.25174                 182.39053
419             2.24899                 179.91904
420             2.24385                 177.26423
421             2.24402                 175.03333
422             2.25827                 173.88710
423             2.38987                 181.63012
424             2.43192                 182.39400
425             2.55201                 188.84844
426             3.11255                 227.21615
427             2.74537                 197.66671
428             2.42016                 171.83150
429             2.38266                 166.78585
430             2.41526                 166.65315
431             2.52072                 171.40910
432             2.54994                 170.84578
433             2.96115                 195.43583
434             2.85865                 185.81238
435             2.60083                 166.45299
436             2.65571                 167.30948
437             2.60763                 161.67275
438             2.61283                 159.38269
439             2.45856                 147.51378
440             2.30478                 135.98184
441             2.35852                 136.79439
442             2.93609                 167.35719
443             3.35506                 187.88347
444             3.35322                 184.42688
445             3.03787                 164.04514
446             3.13412                 166.10820
447             2.97198                 154.54280
448             2.94016                 149.94841
449             3.12721                 156.36060
450             3.04110                 149.01375

 ac convergence test in progress...
../_images/tutorials_05_fm_planets_81_17.png
Auto-corr tau/N = [0.04968753 0.04632667 0.04714928]
tau/N <= 0.02 = [False False False]

451             3.83669                 184.16131
452             3.04320                 143.03026
453             3.14967                 144.88496
454             2.97998                 134.09887
455             2.98020                 131.12867
456             3.19212                 137.26107
457             3.06396                 128.68636
458             3.01083                 123.44415
459             2.97735                 119.09408
460             2.81378                 109.73734
461             2.86901                 109.02230
462             2.75484                 101.92927
463             2.77912                 100.04814
464             2.78004                 97.30133
465             2.76424                 93.98402
466             2.76505                 91.24668
467             2.73225                 87.43194
468             2.81760                 87.34551
469             2.75159                 82.54776
470             2.75168                 79.79869
471             2.75524                 77.14672
472             2.71894                 73.41135
473             2.61586                 68.01249
474             2.84165                 71.04135
475             2.74670                 65.92082
476             2.68962                 61.86121
477             2.92619                 64.37625
478             2.77344                 58.24224
479             2.64283                 52.85662
480             2.59341                 49.27471
481             2.69371                 48.48682
482             2.53827                 43.15054
483             3.05704                 48.91269
484             2.39190                 35.87852
485             2.75769                 38.60762
486             2.70853                 35.21093
487             2.29752                 27.57024
488             2.55598                 28.11574
489             2.60823                 26.08234
490             2.38936                 21.50420
491             2.37861                 19.02885
492             2.65684                 18.59784
493             2.58740                 15.52441
494             2.32334                 11.61671
495             2.45033                 9.80130
496             2.61941                 7.85822
497             2.66607                 5.33214
498             2.49493                 2.49493
499             2.41644                 0.00000
We have reached the limit # of steps without convergence
Running time:  0:21:06.103810
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

If you ran the previous box and wish to write your results, set write=True in the next box. This will pickle the MCMC chain.

[36]:
write=False

if write:
    import pickle
    output = {'chain':chain}
    with open('../datasets/my_MCMC_results', 'wb') as fileSave:
        pickle.dump(output, fileSave)

Note that an alternative to the box above, is to provide output_dir in the call to mcmc_negfc_sampling and set save to True. This will save the results as a pickle including additional keys: apart from ‘chain’, it will also save ‘input_parameters’, ‘AR’ (acceptance ratio), and ‘lnprobability’.

Pickled results can be loaded from disk like this:

[37]:
import pickle
with open('../datasets/MCMC_results','rb') as fi:
    myPickler = pickle.Unpickler(fi)
    mcmc_result = myPickler.load()

print(mcmc_result.keys())
chain = mcmc_result['chain']
dict_keys(['chain', 'input_parameters', 'AR', 'lnprobability'])

The most accurate approach would involve setting a large enough maximum number of iterations and using FFT-based rotation for PCA. The latter in particular, may change a bit the most likely parameter values given the better flux conservation. However, these changes would involve over ~3 orders of magnitude longer calculation time. It is therefore intractable for a personal laptop and not shown in this notebook. If you have access to a supercomputer feel free to adopt these changes though. The results after 500 iterations are nonetheless good enough for illustrative purpose:

5.3.3.2. Visualizing the MCMC chain: corner plots and walk plots

Let’s first check that the walk plots look ok:

[38]:
show_walk_plot(chain)
../_images/tutorials_05_fm_planets_90_0.png

Then based on the walk plot, let’s burn-in the first 30% of the chain, to calculate the corner plots:

[39]:
burnin = 0.3
burned_chain = chain[:, int(chain.shape[1]//(1/burnin)):, :]
show_corner_plot(burned_chain)
../_images/tutorials_05_fm_planets_92_0.png

For the purpose of this tutorial and to limit computation time we set the maximum number of iterations to 500 for 100 walkers. The look of the corner plots may improve with more samples (i.e. higher number of iterations, for a given burn-in ratio). This can be tested by setting the max. number of iterations arbitrarily high for the autocorrelation-time convergence criterion to be met.

5.3.3.3. Highly probable values and confidence intervals

Now let’s determine the most highly probable value for each model parameter, as well as the 1-sigma confidence interval. For this, let’s first flatten the chains:

[40]:
isamples_flat = chain[:, int(chain.shape[1]//(1/burnin)):, :].reshape((-1,3))

Then use the confidence function:

[41]:
val_max, conf = confidence(isamples_flat, cfd=68, gaussian_fit=False, verbose=True, save=False,
                           gt=gt, ndig=1, title=True)
percentage for r: 69.79202279202278%
percentage for theta: 69.02279202279202%
percentage for f: 68.94017094017094%
/Users/valentin/GitHub/VIP/vip_hci/fm/negfc_mcmc.py:1541: UserWarning: The figure layout has changed to tight
  plt.tight_layout(w_pad=0.1)


Confidence intervals:
r: 30.57489764243517 [-0.20421878422770945,0.08279139901123145]
theta: 239.91771004343246 [-0.14527599962792692,0.18317408648738365]
f: 390.7650582608513 [-20.054061636888093,25.285555976945886]
../_images/tutorials_05_fm_planets_98_3.png

Using the confidence function with the gaussian_fit=True option, it is possible to fit a Gaussian to the posterior distribution of each parameter, and infer associated uncertainty values.

[42]:
mu, sigma = confidence(isamples_flat, cfd=68, gaussian_fit=True, verbose=True, save=False,
                       gt=gt, ndig=1, title=True)
percentage for r: 69.79202279202278%
percentage for theta: 69.02279202279202%
percentage for f: 68.94017094017094%


Confidence intervals:
r: 30.57489764243517 [-0.20421878422770945,0.08279139901123145]
theta: 239.91771004343246 [-0.14527599962792692,0.18317408648738365]
f: 390.7650582608513 [-20.054061636888093,25.285555976945886]

Gaussian fit results:
r: 30.525747462344565 +-0.13676103149725966
theta: 239.94793373875024 +-0.15744912409220196
f: 394.1338634773366 +-21.669634439677626
../_images/tutorials_05_fm_planets_100_1.png

It is recommended to report the results as confidence intervals (i.e. with possibly asymmetric uncertainties) as long as the bin interval is small enough. Here, we also fitted the residual posterior distribution of each parameter to a Gaussian distribution (this shape is the expected one if the noise has been well whitened, but is not necessarily guaranteed at all separations depending on the adopted \(n_{\rm pc}\)). In case, these distributions look Gaussian, the inferred \(\sigma\) value may be a more accurate uncertainty value for the different parameters.

We can see that the confidence intervals inferred by NEGFC for the different parameters encompass the ground truth values used for injection (in particular for the flux).

Note that depending on your choice of mu_sigma, you may have to calculate separately another source of uncertainty. Indeed, with the original expression for \(\chi^2\) (Wertz et al. 2017; mu_sigma=False), only the uncertainty associated to photon noise was reflected in the MCMC results. With the new expression (mu_sigma=True), the \(\chi^2\) expression also takes into account the residual speckle noise at the radial separation of the companion candidate. Our tests suggest that similar final uncertainties can be obtained in either of these 2 ways:

  • the uncertainties obtained with MCMC when setting mu_sigma=True;
  • the uncertainties obtained by combining quadratically the uncertainties obtained with MCMC (setting mu_sigma=False and fm = 'sum'), and the residual speckle uncertainties inferred as in Sec. 5.3.4 (Wertz et al. 2017).

5.3.4. Residual speckle uncertainty

Only needed if using ``mu_sigma=False`` in your call to mcmc_negfc_sampling!

Residual speckle noise can also bias the best parameter estimates found for the companion (if not taken into account in the MCMC). To evaluate the uncertainty associated to this additional source of noise, it is recommended to inject a large number of fake companions at the same radius and flux as estimated for the true companion but different azimuths, and then estimate their parameters using simplex-NEGFC. The distribution of differences with respect to injected parameters can then give us an idea of the residual speckle noise uncertainty. This is done in VIP with the speckle_noise_uncertainty function (see also Sec. 3.3 in Wertz et al. 2017 for more details).

Let’s use the planet parameters inferred by the MCMC-NEGFC algorithm:

[43]:
pl_par = (val_max['r'],val_max['theta'],val_max['f'])
pl_par
[43]:
(30.57489764243517, 239.91771004343246, 390.7650582608513)
[44]:
algo_options={'ncomp':opt_npc_ann, 'annulus_width':4*fwhm_naco, 'imlib':imlib_rot,
              'interpolation':interpolation}
speckle_res = speckle_noise_uncertainty(cubefc, pl_par, np.linspace(0,359,360), angs, pca_annulus,
                                        psfn, fwhm_naco, aperture_radius=2, fmerit='sum',
                                        algo_options=algo_options, transmission=None, mu_sigma=None,
                                        wedge=None, weights=None, force_rPA=False, nproc=None,
                                        simplex_options=None, bins=None, save=False, output=None,
                                        verbose=True, full_output=True, plot=True)

#######################################################
###            SPECKLE NOISE DETERMINATION          ###
#######################################################

Number of steps: 360

Process is running for angle: 0.00
Process is running for angle: 18.00Process is running for angle: 36.00

Process is running for angle: 54.00
Process is running for angle: 72.00
Process is running for angle: 37.00
Process is running for angle: 1.00
Process is running for angle: 19.00
Process is running for angle: 55.00
Process is running for angle: 73.00
Process is running for angle: 38.00
Process is running for angle: 2.00
Process is running for angle: 20.00
Process is running for angle: 3.00
Process is running for angle: 56.00
Process is running for angle: 74.00
Process is running for angle: 39.00
Process is running for angle: 21.00
Process is running for angle: 4.00
Process is running for angle: 57.00
Process is running for angle: 75.00
Process is running for angle: 40.00
Process is running for angle: 22.00
Process is running for angle: 58.00
Process is running for angle: 41.00
Process is running for angle: 76.00
Process is running for angle: 5.00
Process is running for angle: 23.00
Process is running for angle: 42.00
Process is running for angle: 6.00
Process is running for angle: 77.00
Process is running for angle: 24.00
Process is running for angle: 59.00
Process is running for angle: 7.00
Process is running for angle: 43.00
Process is running for angle: 25.00
Process is running for angle: 78.00
Process is running for angle: 60.00
Process is running for angle: 8.00
Process is running for angle: 26.00
Process is running for angle: 44.00
Process is running for angle: 61.00
Process is running for angle: 79.00
Process is running for angle: 27.00
Process is running for angle: 45.00
Process is running for angle: 62.00
Process is running for angle: 9.00
Process is running for angle: 46.00
Process is running for angle: 63.00
Process is running for angle: 28.00
Process is running for angle: 80.00
Process is running for angle: 10.00
Process is running for angle: 47.00
Process is running for angle: 64.00
Process is running for angle: 11.00
Process is running for angle: 29.00
Process is running for angle: 81.00
Process is running for angle: 30.00
Process is running for angle: 65.00
Process is running for angle: 48.00
Process is running for angle: 82.00
Process is running for angle: 12.00
Process is running for angle: 31.00
Process is running for angle: 49.00
Process is running for angle: 66.00
Process is running for angle: 13.00
Process is running for angle: 83.00
Process is running for angle: 32.00
Process is running for angle: 50.00
Process is running for angle: 67.00
Process is running for angle: 84.00
Process is running for angle: 14.00
Process is running for angle: 33.00
Process is running for angle: 68.00
Process is running for angle: 51.00
Process is running for angle: 85.00
Process is running for angle: 15.00
Process is running for angle: 34.00
Process is running for angle: 69.00
Process is running for angle: 86.00
Process is running for angle: 16.00
Process is running for angle: 52.00
Process is running for angle: 35.00
Process is running for angle: 70.00
Process is running for angle: 87.00
Process is running for angle: 53.00
Process is running for angle: 17.00
Process is running for angle: 90.00
Process is running for angle: 71.00
Process is running for angle: 108.00
Process is running for angle: 88.00
Process is running for angle: 126.00
Process is running for angle: 91.00
Process is running for angle: 109.00
Process is running for angle: 127.00
Process is running for angle: 89.00
Process is running for angle: 144.00
Process is running for angle: 92.00
Process is running for angle: 128.00
Process is running for angle: 110.00
Process is running for angle: 162.00
Process is running for angle: 145.00
Process is running for angle: 93.00
Process is running for angle: 111.00
Process is running for angle: 163.00
Process is running for angle: 94.00
Process is running for angle: 129.00
Process is running for angle: 112.00
Process is running for angle: 146.00
Process is running for angle: 164.00
Process is running for angle: 130.00
Process is running for angle: 147.00
Process is running for angle: 113.00
Process is running for angle: 165.00
Process is running for angle: 95.00
Process is running for angle: 131.00
Process is running for angle: 148.00
Process is running for angle: 114.00
Process is running for angle: 166.00
Process is running for angle: 132.00
Process is running for angle: 149.00
Process is running for angle: 96.00
Process is running for angle: 115.00
Process is running for angle: 167.00
Process is running for angle: 97.00
Process is running for angle: 150.00
Process is running for angle: 133.00
Process is running for angle: 116.00
Process is running for angle: 168.00
Process is running for angle: 98.00
Process is running for angle: 134.00
Process is running for angle: 151.00
Process is running for angle: 99.00
Process is running for angle: 169.00
Process is running for angle: 135.00
Process is running for angle: 117.00
Process is running for angle: 152.00
Process is running for angle: 100.00
Process is running for angle: 136.00
Process is running for angle: 170.00
Process is running for angle: 118.00
Process is running for angle: 101.00
Process is running for angle: 137.00
Process is running for angle: 119.00
Process is running for angle: 171.00
Process is running for angle: 153.00
Process is running for angle: 102.00
Process is running for angle: 138.00
Process is running for angle: 172.00
Process is running for angle: 120.00
Process is running for angle: 139.00
Process is running for angle: 154.00
Process is running for angle: 173.00
Process is running for angle: 121.00
Process is running for angle: 140.00
Process is running for angle: 103.00
Process is running for angle: 155.00
Process is running for angle: 122.00
Process is running for angle: 141.00
Process is running for angle: 174.00
Process is running for angle: 104.00
Process is running for angle: 156.00
Process is running for angle: 175.00
Process is running for angle: 123.00
Process is running for angle: 105.00
Process is running for angle: 142.00
Process is running for angle: 157.00
Process is running for angle: 124.00
Process is running for angle: 158.00
Process is running for angle: 143.00
Process is running for angle: 176.00
Process is running for angle: 106.00
Process is running for angle: 125.00
Process is running for angle: 159.00
Process is running for angle: 180.00
Process is running for angle: 177.00
Process is running for angle: 198.00
Process is running for angle: 107.00
Process is running for angle: 181.00
Process is running for angle: 160.00
Process is running for angle: 178.00
Process is running for angle: 199.00
Process is running for angle: 216.00
Process is running for angle: 179.00
Process is running for angle: 161.00
Process is running for angle: 182.00
Process is running for angle: 217.00
Process is running for angle: 200.00
Process is running for angle: 234.00
Process is running for angle: 252.00
Process is running for angle: 183.00
Process is running for angle: 218.00
Process is running for angle: 201.00
Process is running for angle: 235.00
Process is running for angle: 253.00
Process is running for angle: 184.00
Process is running for angle: 219.00
Process is running for angle: 202.00
Process is running for angle: 236.00
Process is running for angle: 254.00
Process is running for angle: 185.00
Process is running for angle: 203.00
Process is running for angle: 220.00
Process is running for angle: 237.00
Process is running for angle: 186.00
Process is running for angle: 255.00
Process is running for angle: 204.00
Process is running for angle: 238.00
Process is running for angle: 221.00
Process is running for angle: 187.00
Process is running for angle: 256.00
Process is running for angle: 205.00
Process is running for angle: 239.00
Process is running for angle: 257.00
Process is running for angle: 222.00
Process is running for angle: 206.00
Process is running for angle: 188.00
Process is running for angle: 240.00
Process is running for angle: 258.00
Process is running for angle: 223.00
Process is running for angle: 207.00
Process is running for angle: 189.00
Process is running for angle: 259.00
Process is running for angle: 224.00
Process is running for angle: 241.00
Process is running for angle: 190.00
Process is running for angle: 208.00
Process is running for angle: 225.00
Process is running for angle: 260.00
Process is running for angle: 242.00
Process is running for angle: 191.00
Process is running for angle: 209.00
Process is running for angle: 243.00
Process is running for angle: 261.00
Process is running for angle: 226.00
Process is running for angle: 210.00
Process is running for angle: 244.00
Process is running for angle: 227.00
Process is running for angle: 192.00
Process is running for angle: 211.00
Process is running for angle: 262.00
Process is running for angle: 245.00
Process is running for angle: 193.00
Process is running for angle: 228.00
Process is running for angle: 263.00
Process is running for angle: 246.00
Process is running for angle: 194.00
Process is running for angle: 212.00
Process is running for angle: 229.00
Process is running for angle: 264.00
Process is running for angle: 247.00
Process is running for angle: 213.00
Process is running for angle: 230.00
Process is running for angle: 195.00
Process is running for angle: 265.00
Process is running for angle: 248.00
Process is running for angle: 214.00
Process is running for angle: 231.00
Process is running for angle: 266.00
Process is running for angle: 196.00
Process is running for angle: 249.00
Process is running for angle: 232.00
Process is running for angle: 267.00
Process is running for angle: 215.00
Process is running for angle: 197.00
Process is running for angle: 250.00
Process is running for angle: 268.00
Process is running for angle: 233.00
Process is running for angle: 270.00
Process is running for angle: 251.00
Process is running for angle: 288.00
Process is running for angle: 271.00
Process is running for angle: 306.00
Process is running for angle: 324.00
Process is running for angle: 289.00
Process is running for angle: 269.00
Process is running for angle: 307.00
Process is running for angle: 325.00
Process is running for angle: 272.00
Process is running for angle: 290.00
Process is running for angle: 308.00
Process is running for angle: 326.00
Process is running for angle: 273.00
Process is running for angle: 342.00
Process is running for angle: 291.00
Process is running for angle: 309.00
Process is running for angle: 327.00
Process is running for angle: 343.00
Process is running for angle: 274.00
Process is running for angle: 310.00
Process is running for angle: 328.00
Process is running for angle: 344.00
Process is running for angle: 292.00
Process is running for angle: 275.00
Process is running for angle: 311.00
Process is running for angle: 345.00
Process is running for angle: 293.00
Process is running for angle: 329.00
Process is running for angle: 312.00
Process is running for angle: 276.00
Process is running for angle: 294.00
Process is running for angle: 346.00
Process is running for angle: 330.00
Process is running for angle: 313.00
Process is running for angle: 277.00
Process is running for angle: 295.00
Process is running for angle: 347.00
Process is running for angle: 331.00
Process is running for angle: 278.00
Process is running for angle: 314.00
Process is running for angle: 296.00
Process is running for angle: 332.00
Process is running for angle: 315.00
Process is running for angle: 348.00
Process is running for angle: 279.00
Process is running for angle: 316.00
Process is running for angle: 349.00
Process is running for angle: 333.00
Process is running for angle: 280.00
Process is running for angle: 297.00
Process is running for angle: 334.00
Process is running for angle: 317.00
Process is running for angle: 281.00
Process is running for angle: 350.00
Process is running for angle: 298.00
Process is running for angle: 335.00
Process is running for angle: 282.00
Process is running for angle: 318.00
Process is running for angle: 351.00
Process is running for angle: 336.00
Process is running for angle: 299.00
Process is running for angle: 319.00
Process is running for angle: 283.00
Process is running for angle: 337.00
Process is running for angle: 352.00
Process is running for angle: 320.00
Process is running for angle: 300.00
Process is running for angle: 338.00
Process is running for angle: 284.00
Process is running for angle: 353.00
Process is running for angle: 321.00
Process is running for angle: 301.00
Process is running for angle: 285.00
Process is running for angle: 339.00
Process is running for angle: 354.00
Process is running for angle: 322.00
Process is running for angle: 340.00
Process is running for angle: 286.00
Process is running for angle: 323.00
Process is running for angle: 355.00
Process is running for angle: 302.00
Process is running for angle: 287.00
Process is running for angle: 341.00
Process is running for angle: 356.00
Process is running for angle: 303.00
Process is running for angle: 357.00
Process is running for angle: 304.00
Process is running for angle: 358.00
Process is running for angle: 305.00
Process is running for angle: 359.00
residuals (offsets):  [ 0.01414171  0.0262174  -0.02752407  0.00633408  0.00851382 -0.01750519
 -0.03294881 -0.05479893 -0.06141104 -0.07184825 -0.06401787 -0.07708748
 -0.09329984 -0.06945247 -0.02782975 -0.01587231  0.02667203  0.05267615
  0.10288292  0.07627987  0.08018347  0.12334401  0.04019966 -0.01770346
 -0.05793479 -0.08912758 -0.0794795  -0.14093781 -0.11793176 -0.119779
 -0.10652351 -0.08996851 -0.08150916 -0.01038372 -0.00545641  0.00270321
  0.04474867  0.01589799  0.02295404  0.0140863   0.01795283 -0.0171472
 -0.01237075  0.00517289 -0.01035374  0.03275646  0.03156872  0.07182242
  0.05508761  0.05658805  0.05980567  0.04488462  0.04459095  0.05390772
  0.04833359  0.0420769   0.0273723  -0.00647562 -0.02855583 -0.04734734
 -0.08889904 -0.08824835 -0.08527666 -0.08017446 -0.09020929 -0.03558856
 -0.00263108  0.01303852  0.00397229  0.0003241   0.031516    0.01254575
  0.01684572  0.02448462 -0.00600674  0.01099637  0.01884888  0.01250341
  0.01915659 -0.05120402 -0.02162102 -0.03636504 -0.01526888 -0.03621322
  0.00222169 -0.0158525   0.02130287  0.07028592  0.05980563 -0.10668692
 -0.00839235 -0.11306102 -0.10291854 -0.08462054 -0.0621191  -0.04049092
 -0.0274534  -0.01483692 -0.00782895  0.02105108  0.01492988  0.06281677
  0.05777138  0.09733645  0.00715255  0.01590097  0.01551381  0.01763579
 -0.00535186 -0.00123858 -0.06567812 -0.0272517  -0.06878691 -0.06145989
 -0.06189969 -0.07929754 -0.04947056 -0.08905203 -0.08176747 -0.04172215
 -0.01835565  0.01168755  0.00662449 -0.00123935  0.0255089  -0.00039274
  0.04184577  0.08905198  0.05344004  0.04358257  0.03197981  0.03093149
  0.03826845  0.01343565  0.03666553  0.0044708   0.00737599  0.01446193
 -0.02055789  0.00765128  0.02713428  0.0266011  -0.00693228 -0.02869506
 -0.06813617 -0.07251775 -0.07042657 -0.05246621 -0.05980442 -0.1352383
 -0.02519737 -0.00744193  0.0051315   0.00648816  0.00827969 -0.01457707
 -0.02572175 -0.04071202 -0.04307454 -0.00939839  0.01882996  0.03074528
  0.02937349  0.0333489   0.01294039  0.00785659 -0.02121328  0.0262367
  0.00969314 -0.01652188 -0.05404947 -0.03568632 -0.03962839 -0.06266726
 -0.08633723 -0.07082849 -0.02199952 -0.03905974 -0.01120959  0.02798086
  0.07354342  0.03140005  0.05278003  0.07006672  0.04922261  0.03948324
  0.00366109  0.00320068 -0.03376512 -0.09167506 -0.05655047 -0.01984591
  0.00015404 -0.01793363  0.0142114   0.02020543  0.02085339  0.05366862
  0.07429004  0.05022277  0.05069189  0.03686727  0.05232911  0.01745188
  0.01070037 -0.03965554 -0.00938616 -0.09157149 -0.08866494 -0.08731417
 -0.04725207 -0.08217073 -0.05372192 -0.08561344 -0.10030488 -0.02973702
 -0.0139864   0.00057709  0.01094436  0.04596034  0.06776114  0.10794664
  0.08595627  0.09753041  0.08342088  0.09656291  0.03307842  0.04704147
  0.00989003  0.03983763 -0.00504385 -0.00950893 -0.02334027  0.00687799
 -0.02206765  0.00058425 -0.02133179 -0.0286293  -0.00461759  0.01834439
 -0.01312724 -0.02674027 -0.02486957 -0.01769951 -0.02191695  0.04678306
 -0.00122649 -0.02041071 -0.02153878  0.05965003 -0.00314831 -0.02949817
 -0.04340635 -0.05143624 -0.0430844  -0.04466024 -0.0541319  -0.08655558
 -0.07282757 -0.09032474 -0.04076583  0.02626367  0.04030204  0.05574441
  0.04111735  0.04741127  0.03231528  0.07808787  0.05137703 -0.00359852
  0.00784718  0.00062876 -0.05086493 -0.02292529 -0.05333868  0.00469407
 -0.02433034 -0.03071813 -0.03684427 -0.02011971 -0.02383375 -0.00073175
 -0.03290598 -0.04328762 -0.05100856 -0.04850863 -0.02487256  0.0065488
  0.02546421  0.01782216  0.03123057  0.03979113 -0.01524292 -0.02185833
 -0.04684812 -0.0566921  -0.07006926 -0.06406619 -0.05324135 -0.03124257
  0.07103417  0.06970876  0.10698092  0.12005147  0.12616922  0.07205042
  0.04007435  0.01631929 -0.01106865 -0.04755221 -0.07247758 -0.04021001
 -0.04576157 -0.0728561  -0.07912741 -0.08320474 -0.04389924 -0.03971008
 -0.01691231  0.01856232  0.03230188  0.08903604  0.10283423  0.10298523
  0.1128555   0.09034228  0.02644974  0.0579124   0.0516035   0.03053519
  0.02384122 -0.00603993 -0.02856481 -0.04721303 -0.04102646 -0.04982591
 -0.10188523 -0.08192184 -0.056746   -0.04563918 -0.04630798 -0.03487262
  0.00066623  0.02864871  0.04359411  0.03937825  0.08614059  0.10158204
  0.13035782  0.12599861  0.08409184  0.04292223  0.01278754 -0.02476346
 -0.01492819 -0.05529099 -0.08109458 -0.06363078 -0.04669634 -0.03738522] [ 2.03643038e-04  6.04056345e-02  7.15148194e-02  2.06498111e-02
  3.43528960e-02  1.92689198e-02 -2.32146359e-02  1.57809159e-02
 -2.95959034e-02  3.13600176e-02  8.78503060e-03  2.91826427e-02
 -3.20643518e-02 -5.19404818e-03 -1.75386692e-02 -8.53406205e-03
 -5.28376401e-02  1.74234807e-04 -3.81397377e-02 -1.78067604e-02
  5.32637447e-03 -4.86433796e-02 -4.50112202e-03  7.31847238e-03
 -2.67333732e-03 -3.04248115e-02  2.89893170e-02  5.69220004e-02
  4.15228342e-02  2.46478833e-02 -4.00100550e-03 -2.65103683e-02
 -4.27952374e-02 -9.14967186e-03 -4.22927946e-02  3.77680215e-04
 -2.04263493e-02  2.24371434e-02  2.92873312e-02  2.30897631e-02
  2.75905926e-02  8.74832710e-03  4.51187914e-02  5.25445719e-02
  6.18002070e-02  6.22349904e-02  1.84110750e-02  5.59055961e-02
  2.68830633e-02 -4.67217622e-02 -5.93801272e-02 -1.05506614e-01
 -1.19418820e-01 -1.07944135e-01 -1.11842052e-01 -1.26989265e-02
  6.70746225e-03  3.35881206e-02  9.84788553e-02  1.65844279e-01
  1.19681498e-01  6.98669022e-02  9.55178802e-02  2.94275709e-02
  1.87921302e-03 -1.48187766e-02  6.33951061e-03 -2.23708257e-02
 -2.28641212e-02 -7.51787794e-02 -9.00448602e-02 -1.04615672e-01
 -1.28973848e-01 -9.55473057e-02 -9.86892296e-02 -6.90071945e-02
 -6.25661587e-02 -3.06038477e-02  8.88208702e-02  5.90856424e-02
  9.73346595e-02  1.18506872e-01  1.76414804e-01  1.22095296e-01
  1.40062766e-01  9.69248646e-02  7.45612802e-02  1.29453814e-02
 -2.36717606e-02  3.78799964e-02 -7.13768397e-03 -6.68700839e-04
  1.61364971e-02  1.55504154e-02 -3.03245132e-02 -2.42066481e-02
 -2.83339521e-02 -5.66273078e-03 -9.24539225e-02 -1.07447433e-01
 -6.50379641e-02 -9.63993511e-02 -1.02594025e-01 -4.78073126e-02
 -5.71143936e-02 -5.73119022e-03  7.90750069e-03  6.79323410e-02
  9.84303508e-02  1.27610462e-01  9.45617555e-02  8.78676755e-02
 -2.23552218e-02 -2.80860892e-02 -4.31649938e-02 -1.26126324e-02
 -5.95048548e-02  3.48089464e-03  7.66239009e-03  3.00911747e-03
  2.29781814e-02  1.20557140e-02  5.77597541e-02  2.13403150e-03
  4.71912585e-03 -2.81731909e-03 -2.25736128e-03  1.92446632e-03
 -6.32995962e-02 -1.96391671e-02  2.07201189e-02  3.06192249e-02
  6.79247802e-02  6.43607136e-02  7.25270510e-02  5.00583327e-02
  8.11297986e-02  6.25848054e-02  5.70107104e-02  4.05597071e-02
 -3.71796688e-03  1.43312286e-02 -5.31489768e-02 -1.13935319e-02
 -2.72397466e-02 -7.23426240e-02 -9.11867795e-02 -7.51726329e-02
 -7.17502062e-02 -9.52528077e-02 -1.36004298e-02  4.64664994e-02
  3.98523710e-02  1.16390570e-01  8.33659694e-02  5.47396837e-02
  7.01762008e-03 -6.35136258e-03 -4.92262411e-02 -4.92755603e-02
 -4.11331478e-02 -6.92442723e-02 -7.29559437e-02 -2.36210229e-02
 -3.58518421e-02  2.89076040e-02  6.54248385e-02  8.65840585e-02
  7.48527526e-02  6.50025366e-02  7.15796805e-02  1.63837628e-02
 -2.35828639e-02 -1.00892468e-01 -9.54976296e-02 -1.88200499e-01
 -1.32562722e-01 -1.25464525e-01 -6.05824622e-02 -9.02201347e-04
  9.55371359e-02  9.72649544e-02  4.45310872e-02  8.33533404e-02
  2.27797028e-02 -4.81986779e-02 -5.13054984e-02 -9.11292976e-02
 -6.26394489e-02 -1.92513093e-02  2.19246813e-03 -2.04420049e-03
  1.43533600e-03  2.30623730e-02 -2.16641053e-02 -2.39032821e-02
  1.40286242e-02  4.07161251e-02  3.10363829e-02  5.24045999e-02
  5.72339362e-02  2.56012727e-02 -8.57718367e-03 -3.54395323e-02
 -2.15445008e-02 -3.96315439e-02 -4.43152526e-02 -2.94764217e-02
  2.62935473e-02  8.67183780e-03  1.08944412e-02  3.83453464e-02
  1.44936334e-02 -2.36900571e-03  1.12531278e-02 -2.22154716e-02
 -5.68184267e-04  1.45833563e-03  2.05545483e-02  3.20206168e-02
  6.22514352e-03 -1.04006381e-02  1.13911618e-02 -4.13978907e-02
  3.90415571e-02 -1.64365528e-02  4.71940729e-02  1.08361781e-01
  4.62650541e-02  2.90854030e-02  1.92660505e-03 -9.95344122e-03
 -5.73910034e-02 -6.45246816e-02 -6.16377100e-02  7.74349213e-03
 -4.73007665e-03  1.04807699e-02  3.64080762e-02  6.16104214e-02
 -2.91159532e-02 -4.80991156e-03 -8.32892470e-02 -3.74854495e-02
 -4.74161814e-02 -1.04048173e-01 -4.05056078e-02  2.03383147e-02
  1.09600905e-02  5.51290351e-02  5.36702832e-02  4.96720974e-02
  3.06858001e-02 -2.61961456e-02 -4.17524011e-02 -8.92086758e-02
 -1.17708498e-01 -9.36077484e-02 -5.55982441e-02 -1.68039861e-02
  5.67615306e-02  6.51231322e-02  1.50934271e-01  1.65575706e-01
  1.66743691e-01  1.29922525e-01  4.55175860e-02 -3.83921922e-02
 -8.67744311e-02 -1.96721186e-01 -2.39042072e-01 -2.97479284e-01
 -2.91613002e-01 -2.01543614e-01 -1.28326955e-01 -3.34181980e-02
 -3.82931663e-04  5.64501106e-02  1.36173331e-01  9.77770716e-02
  3.73373887e-02 -3.65700166e-02 -9.75443494e-03  2.32698898e-03
  1.73804345e-01  1.80830081e-01  1.69188105e-01  1.10104094e-01
  1.63048812e-01  8.09333022e-02 -8.78703899e-03 -7.27250293e-02
 -1.02576236e-01 -1.31234892e-01 -1.87326204e-01 -1.31535266e-01
 -1.49462921e-01 -9.20203277e-02 -1.37145651e-02  7.55277718e-02
  1.29401764e-01  9.28656142e-02  1.22825788e-01  9.69359169e-02
  7.41231922e-02  3.93379673e-02  1.66978050e-02 -3.68099187e-02
 -2.98841191e-02 -4.57764202e-02 -5.62008555e-02 -7.22843638e-02
 -6.11601937e-03 -9.60955183e-03 -3.96788076e-02 -1.58527167e-02
  3.08294867e-02  5.23420046e-02  5.23392734e-02  3.02751470e-02
  1.97240085e-02  3.32237382e-02 -4.02080330e-02 -2.84669587e-02
 -5.72788259e-02 -3.96696175e-02 -4.90873969e-03  4.26421192e-02
  6.93441429e-02  6.21695047e-02  5.46680977e-02  2.52276948e-02
  1.01504763e-01  5.70233558e-02  1.39772048e-02  5.28382846e-03
 -2.79881014e-02 -3.63417119e-02 -2.67575610e-02 -7.19982593e-02
 -6.91638890e-02 -6.36869854e-02 -4.51738349e-02 -3.60339298e-02
 -3.32294804e-03  1.69095583e-02  1.01050718e-01  7.52340504e-02
  5.10626495e-02  7.03337637e-02  4.81707663e-02 -2.13872486e-02
 -4.67918358e-02 -1.53715715e-01 -1.62017655e-01 -1.59213346e-01
 -1.23993149e-01 -5.16842711e-02  3.55498781e-02  1.25588714e-01] [-1.21535789e+01 -8.95141500e+00 -5.37132590e+00 -4.15038104e+00
  2.91852577e+00 -3.46349538e+00 -2.59925234e+00  2.96288383e+00
 -1.92204197e+00  1.49864556e+00  7.87531367e-01  4.73316591e+00
  8.55648417e+00  2.65499273e+00 -1.63454697e+00  2.78853122e-01
 -3.27524280e+00 -3.80845389e+00 -8.68164258e+00 -5.87270401e-01
 -2.08697723e+01 -1.26290329e+01 -7.74417786e+00 -8.15714299e+00
 -9.42300729e+00 -3.29461052e+00 -5.89056725e+00  3.20651451e+00
  4.32157869e+00  3.51621295e+00  3.55377814e+00  4.36934888e+00
  7.83393901e-02 -2.75787711e+00 -5.34874229e+00 -4.40827035e+00
 -4.64256336e+00 -7.38126878e-01  9.55727095e-01  4.59746538e+00
 -2.64774185e+00  2.66852428e+00  1.41776361e+00  2.59246685e+00
  1.14697284e+01  7.70201451e+00  1.14399201e+01  5.79889291e+00
  9.82519353e+00  9.56584372e+00  6.46348751e+00 -3.11046403e+00
 -5.92209340e+00 -1.05811956e+01 -1.89463069e+01 -1.92918321e+01
 -1.94574968e+01 -2.23672961e+01 -1.43862480e+01 -8.40604099e+00
 -6.39582257e+00 -5.58309781e+00  1.21001706e+00  3.17459742e+00
  8.90485680e+00  1.20409746e+01  1.36529507e+01  7.28342767e+00
  6.68703256e+00  7.93712386e+00  1.36187051e+00 -2.51985345e+00
 -3.77838934e+00 -9.92431844e+00 -1.25511802e+01 -1.79810503e+01
 -2.23365827e+01 -2.27569690e+01 -1.93978738e+01 -2.48953506e+01
 -1.83351792e+01 -1.46546535e+01 -1.16235064e+01  3.10166078e+00
  5.94084564e+00  1.09092789e+01  1.82972030e+01  2.32951202e+01
  2.07681023e+01  2.96029557e+01  1.77229963e+01  2.29747076e+01
  1.86304997e+01  1.70426614e+01  1.74812647e+01  1.34997396e+01
  9.52577682e+00  6.26422515e+00  7.09185534e+00 -4.43313287e-01
 -2.32373589e-01 -5.53500566e+00 -1.80476268e+01 -2.12974251e+01
 -6.86390262e+00 -1.58716637e+01 -1.16031886e+01 -1.04061728e+01
 -1.17857454e+01 -3.35033045e+00 -1.73452649e+00  2.54701988e+00
  5.65753310e+00  8.18429542e-01 -1.14445435e+00 -4.55010754e+00
 -6.54442116e+00 -5.47986901e+00 -8.07007293e+00 -5.38303847e+00
 -7.37748114e+00 -3.40506459e+00 -3.78781371e+00  8.76421939e-01
 -1.81432262e-01  9.72665465e-02 -1.26336983e-01 -4.08470630e+00
 -4.96295794e-01 -3.66414351e+00 -6.92731580e+00 -7.32308183e+00
 -3.52803888e+00 -4.13899229e-01  2.02844233e+00 -3.00791795e-03
  7.17519879e+00  1.37434542e+01  1.18346976e+01  1.30479712e+01
  1.61480088e+01  1.36970587e+01  1.42464388e+01  7.87110803e+00
  7.05312003e+00  7.94223428e+00  2.60302814e+00  5.72313904e-01
 -4.51569096e+00 -7.39999636e+00 -6.31805088e+00 -5.69191130e+00
 -6.65914540e-01  8.96304162e+00  9.51334800e+00  9.34879244e+00
  1.54613739e+01  9.70688993e+00  1.21391229e+01  8.12644771e+00
  7.01634234e+00  2.65422027e+00  2.38271950e+00  1.11826880e+00
 -1.72991814e+00  2.88764669e+00  4.19115859e+00  4.07910272e+00
  8.64708036e+00  1.11493720e+01  1.45037696e+01  1.51203397e+01
  1.93096455e+01  1.58130286e+01  7.03478576e+00 -1.09180267e-01
 -1.21049476e+00 -5.04366055e+00 -9.06993398e+00 -7.60378198e+00
 -6.21366003e-01 -8.45347907e-01  4.81028268e+00  2.94101862e+00
  8.04742847e+00  6.66700906e+00  3.63884836e+00  2.04664858e+00
 -9.12822643e-01  1.21207246e+00 -3.68067622e+00 -1.18960687e+00
 -1.71297454e+00  2.01440159e+00 -4.92689313e+00 -3.79022159e+00
 -4.86220150e+00 -3.99545699e+00 -3.51158085e+00  1.53882049e+00
  2.57845372e+00  3.24678498e+00  1.78586037e+00  4.58754590e+00
  1.52184360e+00  6.57111489e-01 -2.42374278e+00  1.70627162e-01
 -4.50206612e-01 -4.57976950e-01  6.05590511e-01 -4.37754564e+00
 -4.14785398e+00 -4.70886567e+00 -1.03825534e+00 -7.08613361e+00
 -5.94175798e+00 -4.40656766e-01 -7.74084164e+00 -6.31381107e+00
 -4.17438682e+00  2.13305471e+00  1.15978724e+00  3.91268163e+00
 -1.54937621e+00 -5.54479464e+00  7.24492034e-01  5.56585352e+00
  6.90164108e+00  1.15053394e+01  9.36015865e+00  1.14735160e+01
  3.97062797e+00  3.72122151e+00  1.91061647e-01  3.58569443e+00
  1.85459515e-01 -4.29199223e+00  3.77561928e+00  2.44289531e-01
  1.52594796e+00  1.37930234e+00 -4.02128431e+00 -6.47788057e-01
 -4.49350583e+00 -7.57904179e+00 -1.10415894e+01 -1.03728427e+01
 -1.25704769e+01 -5.93037282e+00 -7.13157837e+00 -3.63441618e+00
  1.11691044e+00  3.18821947e+00  1.99769110e-01 -1.74653886e+00
 -2.60104941e+00 -5.65973386e+00 -1.27161011e+01 -1.69575090e+01
 -1.54728558e+01 -7.57831344e+00 -4.58024093e+00  3.84218848e+00
  1.47006896e+01  1.91164297e+01  3.02614107e+01  3.27701897e+01
  2.87721631e+01  1.56508424e+01  1.42729051e+01  2.91558077e+00
 -5.38312362e+00 -1.94393121e+01 -2.69399929e+01 -3.59197325e+01
 -3.56152892e+01 -3.48120647e+01 -2.95379858e+01 -2.44661724e+01
 -1.91289487e+01 -1.29724261e+01 -1.32505444e+01 -1.37166238e+01
 -1.58154563e+00  1.11369229e+00  8.55720623e-01  7.99209022e+00
  1.98294386e+01  2.04510928e+01  2.42940408e+01  2.58947744e+01
  1.69767938e+01  6.01389981e+00  8.05268820e+00 -1.39401472e+00
 -1.05155329e+01 -1.81174519e+01 -1.78686813e+01 -2.21444433e+01
 -1.12858850e+01 -1.18053323e+01 -8.58170406e+00  6.51280542e+00
  9.46144059e+00  1.23457962e+01  1.48173162e+01  1.19734882e+01
  9.53050687e+00  9.06338828e+00  4.23352037e+00  3.28171933e+00
  7.41549049e+00 -8.30567965e+00 -4.29381774e+00 -2.11226062e+00
 -2.00593161e+00 -1.09636442e+00  3.87310558e+00  5.41725702e+00
  4.47119642e+00  6.63772859e+00  4.21849862e+00 -3.41724851e-01
 -1.52616657e+00 -5.80214278e+00 -9.38906751e+00 -8.05807558e+00
 -2.10984878e+00  1.96354579e+00  4.90250850e+00  3.79646330e+00
  5.51121580e+00  6.76037233e+00  1.02417189e+01  1.10307347e+01
  1.05984669e+01  9.09027597e+00  6.32608451e+00  7.39601020e+00
  3.66191486e+00  5.46929987e-01 -2.86699071e+00  1.93603139e+00
 -2.13760876e+00 -4.31180537e+00 -3.51062816e-01  5.25554600e+00
  8.40249990e+00  9.53984451e+00  9.91353680e+00  9.66477084e+00
  1.21796192e+01  3.28826282e+00 -2.67936655e+00 -6.47135703e+00
 -1.92797874e+01 -2.28661560e+01 -1.96652376e+01 -1.58607566e+01]
[[ 1.41417065e-02  2.03643038e-04 -1.21535789e+01]
 [ 2.62173985e-02  6.04056345e-02 -8.95141500e+00]
 [-2.75240652e-02  7.15148194e-02 -5.37132590e+00]
 ...
 [-6.36307850e-02 -5.16842711e-02 -2.28661560e+01]
 [-4.66963397e-02  3.55498781e-02 -1.96652376e+01]
 [-3.73852196e-02  1.25588714e-01 -1.58607566e+01]]
(360, 3)
percentage for r: 68.88888888888887%
percentage for theta: 69.44444444444446%
percentage for f: 71.11111111111113%


Confidence intervals:
r: -0.02412997296552928 [-0.056519922705534256,0.07912789178774797]
theta: 0.0014640689854274963 [-0.07307504175148613,0.07307504175148616]
f: -0.620744716053828 [-10.4942936707127,10.4942936707127]

Gaussian fit results:
r: -0.005358008112249547 +-0.05277034604807372
theta: -0.000828655169225096 +-0.07269053838995253
f: -0.05129156730094909 +-10.60503438051984
../_images/tutorials_05_fm_planets_106_1.png

Again, if you wish to write the result (to avoid having to run the previous box again), just set write=True. This will pickle the result:

[45]:
write=False

if write:
    output = {'speckle_res':speckle_res}
    with open('../datasets/my_speckle_residuals_results', 'wb') as fileSave:
        pickle.dump(output, fileSave)

Load pickled results from disk:

[46]:
with open('../datasets/speckle_residuals_results','rb') as fi:
    myPickler = pickle.Unpickler(fi)
    sp_unc_result = myPickler.load()

print(sp_unc_result.keys())

speckle_res = sp_unc_result['speckle_res']
sp_unc, mean_dev, p_simplex, offset, chi2, nit, success = speckle_res
dict_keys(['speckle_res'])

The speckle uncertainty associated to each parameter is contained in sp_unc which corresponds to the 1\(\sigma\) width of a Gaussian distribution fitted to the offset distribution (i.e. the differences with respect to injected ground truths):

[47]:
sp_unc
[47]:
array([ 0.05350939,  0.07376633, 10.64773837])

For comparison, the uncertainties found by the MCMC procedure were:

[48]:
sigma
[48]:
array([ 0.13676103,  0.15744912, 21.66963444])

5.3.5. Final uncertainties

The final uncertainties on the planet parameters should include both statistical and systematic uncertainties. The former include both the photon noise and residual speckle noise uncertainties discussed above. The latter include both the uncertainty on the star location (which may be non-negligible when using a coronagraph) and instrumental calibration errors, including:

  • the uncertainty on the plate scale (for \(r\), when converting to angular separation) - note that it is proportional to the radial separation itself;
  • the uncertainty on the PA of true north (for \(\theta\)).

The uncertainty on the star location is of the order of 0.3px in individual NACO+AGPM images (this is the typical precision by manual recentering during the observation). Given the shift plots in Tutorial 2, it appears the autocorrelation timescale is of the order of ~5 frames. Therefore, considering that there are 61 frames in the datacube, the uncertainty on the star location in the final combined image must be roughly:

[49]:
cen_unc_indiv = 0.3 #px
cen_unc = cen_unc_indiv/np.sqrt(61/5) #px
cen_unc
[49]:
0.08588975014708022

The latter can be translated into an uncertainty on \(\theta\) by division by the radial separation of the companion. The stellar centering uncertainties on each planet parameter can thus be expressed as:

[50]:
star_unc = np.array([cen_unc*pxscale_naco,
                     np.rad2deg(cen_unc/val_max['r']),
                     0]) # uncertainty on each of the 3 parameters due to stellar centering

where the multiplication by pxscale_naco has converted the radial separation in arcsec.

For the instrumental calibration errors, we adopt the values quoted in Absil et al. (2013). Note that the uncertainty related to the plate scale is directly proportional to the radial separation of the companion.

[51]:
dr_unc = 0.00004 # plate scale uncertainty in arcsec per px
tn_unc = 0.09    # deg
syst_unc = np.array([val_max['r']*dr_unc,
                     tn_unc,
                     0])

The final uncertainties are then the different sources of uncertainty added quadratically - after conversion of the radial separation to arcsec:

[52]:
sigma[0] *= pxscale_naco
sp_unc[0] *= pxscale_naco

if mu_sigma:
    final_unc = np.sqrt(np.power(sigma,2) + np.power(star_unc,2) + np.power(syst_unc,2))
else:
    final_unc = np.sqrt(np.power(sigma,2) + np.power(sp_unc,2) + np.power(star_unc,2) + np.power(syst_unc,2))
[53]:
msg = "The final estimates on the radial separation, PA and flux of the companion, including uncertainties, are: \n"
msg+= "r = {:.2f}+-{:.2f} mas (GT: {:.2f} mas), \n"
msg+= "PA = {:.2f}+-{:.2f} deg (GT: {:.2f} deg) \n"
msg+= "f = {:.2f}+-{:.2f} ADUs (GT: {:.2f} ADUs)"
print(msg.format(val_max['r']*pxscale_naco*1000, final_unc[0]*1000, rad_fc*pxscale_naco*1000,
                 val_max['theta'], final_unc[1], theta_fc,
                 val_max['f'], final_unc[2], flux_fc))
The final estimates on the radial separation, PA and flux of the companion, including uncertainties, are:
r = 831.33+-4.56 mas (GT: 829.29 mas),
PA = 239.92+-0.24 deg (GT: 240.00 deg)
f = 390.77+-21.67 ADUs (GT: 400.00 ADUs)

Let’s consider the Gaussian fit instead:

[54]:
msg = "Considering a Gaussian fit to the posterior distributions, the final estimates on the radial separation, PA and flux of the companion, including uncertainties, are: \n"
msg+= "r = {:.2f}+-{:.2f} mas (GT: {:.2f} mas), \n"
msg+= "PA = {:.2f}+-{:.2f} deg (GT: {:.2f} deg) \n"
msg+= "f = {:.2f}+-{:.2f} ADUs (GT: {:.2f} ADUs)"
print(msg.format(mu[0]*pxscale_naco*1000, final_unc[0]*1000, rad_fc*pxscale_naco*1000,
                 mu[1], final_unc[1], theta_fc,
                 mu[2], final_unc[2], flux_fc))
Considering a Gaussian fit to the posterior distributions, the final estimates on the radial separation, PA and flux of the companion, including uncertainties, are:
r = 830.00+-4.56 mas (GT: 829.29 mas),
PA = 239.95+-0.24 deg (GT: 240.00 deg)
f = 394.13+-21.67 ADUs (GT: 400.00 ADUs)