No need to add "x" (the average of aperture fluxes) for contrast calculation

[VChristiaens]

Hi Markus,

I recently re-inspected the definition of the contrast curve in VIP and the definition in Mawet et al. (2014). After some head scratching, I came to the conclusion that there was a confusion between what the steps in the paper describe and how the actual contrast achieved in an HCI image should be calculated. In my interpretation, the article describes the steps needed to calculate the minimum flux a test aperture should have in a processed image to consider a significant detection there (for the chosen FPF), which involves adding "x" (the average of the fluxes measured in the apertures) in the last step.

However to get the contrast, "x" should not be added. Instead, I'd argue it should just be: contrast = tau s sqrt(1 + 1/n_ap) / throughput / star_photometry

In general x ~ 0, so it usually doesn't matter much, except in some symptomatic cases where I noted this can make a significant difference, in particular for classical full-frame PCA RDI. When an imperfectly correlated PSF is provided, oversubtraction can lead to negative annular regions, where the estimated contrast then gets way too optimistic, as a negative x gets added. In the worst cases, x could even get larger in absolute value than s * sqrt(1 + 1/n_ap), and lead to negative contrasts.

Another way to see this is that the actual (throughput-affected) companion flux in processed images is its aperture flux minus the average of aperture fluxes at the same radius (bar{x1}-bar{x2} in the notations of the paper) - be it in symptomatic cases or not. Therefore adding "back" x (x2) is not necessary, when calculating the achieved contrast.

I remember you mentioned the missing "x" in a VIP issue, which prompted a change on our side. I have now reverted the change in the latest release of VIP. I thought I'd let you know here as well since this also concerns how Applefy calculates contrasts.

Cheers, Valentin

[markusbonse]

Hi Valentin,

thanks for checking this again! An imperfect correction of the PSF in RDI is definitely a problem for the t-test. If you see negative structures in the images after PSF subtraction your residual noise is probably highly non-Gaussian and highly non-independent in which case the test is not really applicable. If there is not way to remove the systematic noise not adding back x might be a work around. But from the statistics point of view it is important.

The main reason for this is the choice of the test: A two sample t-test. The test takes two noisy samples X and Y. In our case the sample which contains the signal of the planet only contains a single value: Y_1 = signal after PCA + noise at the position of the planet.

The second part "noise at the position of the planet" is important as the test accounts the uncertainty in both samples X and Y. If we would know the true noise-free signal of the planet we should rather use a one sample t-test and test if our noise sample X is significantly different from the true planet signal.

Given only the residual image we don't have the noise free planet signal, so a one sample t-test is not applicable. For contrast curve calculations it would be possible. But the test statistic of a one sample test is very different changing the whole statistic. If we want our contrast curves to be consistent with test we apply to our residual images we have to add back x to make sure both samples are noisy (as assumed by the test).

But there are two additional aspects which might be important for your RDI scenario:

1. The placement of the apertures: How you place the apertures in your image has a large effect on the final contrast. Applefy accounts for this by placing apertures at various different positions.
2. We have to be careful how we compute the throughput: The throughput is the faction of the planet light that remains after PCA. This means that for the calculation of the throughput x has to be subtracted to get the attenuation of the signal in a noise free environment. If we don't subtract x in this step the curves are probably biased in cases like RDI with a bad PSF subtraction.

I hope these comments help. I'll keep the implementation of applefy unchanged. But I'm happy to discuss further.

Cheers, Markus

[VChristiaens]

Hi Markus,

Thanks for your detailed reply - and sorry for the delay in mine. I agree with an important point you raise: the assumptions to apply the t-test may not be met in some symptomatic cases, and this is definitely something to be cautious about. In some cases, the results may just be meaningless because of that.

That being said, one could in principle have the case of an algorithm leaving non-zero mean residual in concentric annuli, that still correspond to aperture fluxes that are i.i.d and normally distributed - just with a non-zero mu. I am not saying this is necessarily the case of PCA-RDI, but for the sake of better interpreting the role of "x", let us consider that case. The only conditions for the test to be valid are i.i.d and normally distributed samples, but there is no requirement for mu to be zero.

I agree that the addition of x, the estimator of mu, is required to calculate the flux an aperture should have to correspond to a significant detection for a given FPF threshold, considering the 2-sample t-test statistics. But my point is that this aperture flux is not the companion flux - the latter is rather the aperture flux minus mu (or well its estimator x for not knowing the exact value of mu), due precisely to the fact that the mean of the noise distribution could be non-zero.

In essence, I’m suggesting that for the same reason you said x should not be added back when estimating the throughput, I don't think it should be added back when it comes to the actual companion flux to be considered in the contrast definition - i.e. what needs to be considered for the ratio with the stellar flux. Don’t you think?

Cheers, Valentin