L Montier

S Plaszczynski

F Levrier

M Tristram

D Alina

I Ristorcelli

J.-P Bernard

V Guillet

Polarization measurements analysis II. Best estimators of polarization fraction and angle

Keywords: Polarization -Methods, data analysis -Methods, statistical

With the forthcoming release of high precision polarization measurements, such as from the Planck satellite, it becomes critical to evaluate the performance of estimators for the polarization fraction and angle. These two physical quantities suffer from a well-known bias in the presence of measurement noise, as has been described in part I of this series. In this paper, part II of the series, we explore the extent to which various estimators may correct the bias. Traditional frequentist estimators of the polarization fraction are compared with two recent estimators: one inspired by a Bayesian analysis and a second following an asymptotic method. We investigate the sensitivity of these estimators to the asymmetry of the covariance matrix which may vary over large datasets. We present for the first time a comparison among polarization angle estimators, and evaluate the statistical bias on the angle that appears when the covariance matrix exhibits effective ellipticity. We also address the question of the accuracy of the polarization fraction and angle uncertainty estimators. The methods linked to the credible intervals and to the variance estimates are tested against the robust confidence interval method. From this pool of polarization fraction and angle estimators, we build recipes adapted to different use-cases: we provide the best estimators to build a mask, to compute large maps of the polarization fraction and angle, and to deal with low signal-to-noise data. More generally, we show that the traditional estimators suffer from discontinuous distributions at low signal-to-noise ratio, while the asymptotic and Bayesian methods do not. Attention is given to the shape of the output distribution of the estimators, and is compared with a Gaussian distribution. In this regard, the new asymptotic method presents the best performance, while the Bayesian output distribution is shown to be strongly asymmetric with a sharp cut at low signal-to-noise ratio. Finally, we present an optimization of the estimator derived from the Bayesian analysis using adapted priors.

Introduction

The complexity of polarization measurement analysis has been described by Serkowski (1958) when discussing the presence of a systematic bias in optical measurements of linear polarization from stars, and then by Wardle & Kronberg (1974) addressing the same issue in the field of radio astronomy. The bias of polarization measurements happens when one is interested in the polarization intensity P ≡ (Q 2 + U 2 ) or the polarization fraction p ≡ P/I and the polarization angle ψ = 1/2 atan(U/Q) (where I, Q and U are the Stokes parameters), quantities which become systematically biased in the presence of noise. Working with the Stokes parameters Q and U as far as possible avoids this kind of bias. Once a physical modeling of p and ψ is available, and can be translated into Q and U, a likelihood analysis can be performed directly on the Stokes parameters. For the other cases, where no modeling is available, Simmons & Stewart (1985) proposed the first compilation and comparison of methods to deal with the problem of getting unbiased polarization estimates of the polarization fraction and angle, with their associated uncertainties. Then Naghizadeh-Khouei & Clarke (1993) extended the work of Simmons & Stewart (1985) to the characterisation of the polarization angle uncertainties, and Vaillancourt (2006) proposed a method to build confidence limits on polarization fraction measurements. More recently, Quinn (2012) suggested us-ing a Bayesian approach to get polarization estimates with low bias. In all these studies the authors made strong assumptions: no noise on the intensity I and no correlation between the Q and U components, which were also assumed to have equal noise properties. Montier et al. (hereafter PMA I, 2013 in preparation) have quantified the impact of the asymmetry and the correlation between the Q and U noise components on the bias of the polarization fraction and angle measurements. They have shown that the asymmetry of the noise properties can not be systematically neglected as is usually done, and that the uncertainty of the intensity may significantly affect the polarization measurements in the low signal-to-noise (SNR) regime.

In the context of the new generation of polarization data, such as Planck1 (Planck Collaboration I 2011), Blast-Pol (The Balloon-borne Large Aperture Submillimeter Telescope for Polarimetry, [START_REF] Fissel | Society of Photo-Optical Instrumentation Engineers[END_REF], PILOT [START_REF] Bernard | EAS Publications Series[END_REF] or ALMA (Pérez-Sánchez & Vlemmings 2013), which benefit from a much better control of the noise properties, it is essential to take into account the full covariance matrix when deriving the Montier et al.: Polarization measurements analysis II

f (I, p, ψ | I 0 , p 0 , ψ 0 , Σ) = 2|p| I 2 (2π) 3 σ 3 exp            - 1 2         
I -I 0 p I cos 2ψp 0 I 0 cos 2ψ 0 p I sin 2ψp 0 I 0 sin 2ψ 0

         T Σ -1          I -I 0 p I cos 2ψ -p 0 I 0 cos 2ψ 0 p I sin 2ψ -p 0 I 0 sin 2ψ 0                     , (1)
f 2D (p, ψ | I 0 p 0 , ψ 0 , Σ p ) = p πσ 2 p exp         - 1 2       p 2 cos
              , (2)
polarization measurement estimates. In recent works no correction for the bias of the polarization fraction were applied (e.g. Dotson et al. 2010), or only high SNR data were used for analysis (>3) to avoid these issues (e.g. Vaillancourt & Matthews 2012). Two issues are immediately apparent. First, this choice of the SNR threshold may not be relevant for all measurements, and the asymmetry between the orthogonal Stokes noise components could affect the threshold choice. Secondly, the question remains of how to deal with low signal-to-noise data. Using simply the measurements of the polarization parameters (we will call them the "naïve" ones) as estimators of the true values leads to very poor performance, as they lack any information on the noise power. Instead, we would like to perform some transformation on the polarization parameters, in order to remove bias and improve the variance.

This work is the second of a series on the 'Analysis of polarization measurements'. Its aim is to describe how to recover from a measurement (p, ψ) the true polarization fraction p 0 and polarization angle ψ 0 with their associated uncertainties, taking into account the full covariance matrix Σ. We will compare the performance of the various estimators available, and study the impact of the correlation and ellipticity of the covariance matrix on these estimates. We stress that we adopt a frequentist approach to investigate the properties of these estimators, even when dealing with the method inspired by the Bayesian analysis. This means that the estimators are defined as single value estimates, instead of considering the probability density function (pdf) as the proper estimate, as it is usually done in Bayesian methods. The performance of these estimators will be evaluated using three main criteria: the minimum bias, the smallest risk function, and the shape of the distribution of the output estimates. The choice of the most appropriate estimator may vary with the application at hand, and a compromise among them may be chosen to achieve good overall performance. Throughout this work we will make the following two assumptions: i) circular polarization is assumed to be negligible, and ii) the noise on Stokes parameters is assumed to be Gaussian. We also define four regimes of the covariance matrix to quantify its asymmetry, in terms of effective ellipticity (ε eff ) as described in PMA I: the extreme (1<ε eff <2), the low (1<ε eff <1.1), the tiny (1<ε eff <1.01) and the canonical (ε eff =1) regimes.

The paper is organized as follows: we first review in Sect. 2 the expression and the limitations of the polarization estimators, which are extended to take into account the full covariance matrix. We discuss in Sect. 3 the meaning of the polarization uncertainties and we present the different uncertainty estimators. We then compare the performance of the estimators of the polarization fraction in Sect. 4, and of the polarization angle in Sect. 5. In Sect. 6, we discuss some aspects of the problem when the total intensity I is not perfectly known. We conclude with general recipes in Sect. 7.

Polarization estimators

Early work on polarization estimators was based on the Rice (1945) distribution which provides the probability to find a measurement p, for a given true value p 0 and the noise estimate σ p of the Q and U Stokes parameters. The noise values of the Stokes parameters were assumed to be equal (σ p =σ Q /I 0 =σ U /I 0 ), and the total intensity was assumed to be perfectly known, I = I 0 . As we would like to include the full covariance matrix, we use the generalized expression of the pdf from PMA I, which provides the probability to get the measurements (I,p, ψ), given the true values (I 0 , p 0 ,ψ 0 ) and the covariance matrix Σ. Following the notations of PMA I, the expression of the pdf in 3D, including the intensity terms, denoted f (I, p, ψ|I 0 , p 0 , ψ 0 , Σ), is given by Eq. 1, and the pdf in 2D, f 2D (p, ψ|I 0 , p 0 , ψ 0 , Σ p ), by Eq. 2 when the intensity I 0 is assumed to be perfectly known. We also note the introduction of the covariance matrix reduced in 2D,

Σ p = 1 I 2 0 σ 2 Q σ QU σ QU σ 2 U = σ 2 p,G 1 -ρ 2      ε ρ ρ 1/ε      , (3)
where ε = σ Q /σ U is the ellipticity and ρ = σ QU /σ Q σ U is the correlation between the Q and U noise components, leading to an effective ellipticity given by:

ε eff = 1 + ε 2 + (ε 2 -1) 2 + 4ρ 2 ε 2 1 + ε 2 -(ε 2 -1) 2 + 4ρ 2 ε 2 . ( 4
)
With these notations we have Det(Σ p ) = σ 4 p,G and

σ 2 p,G = σ 2 Q I 2 0 1 -ρ 2 ε , (5)
which represents the equivalent radius of a circular Gaussian distribution with the same integrated area as the elliptical one. We also define σ p =σ Q /I 0 =σ U /I 0 when ε eff =1. Finally the pdfs of p and ψ, f p and f ψ , are obtained by marginalization of f 2D over ψ and p, respectively. The expressions for the 1D pdfs f p and f ψ depend on the full set of initial parameters (I 0 , p 0 , ψ 0 ) in the general case, unlike the case under the canonical simplifications (see appendix C of PMA I for fully developed analytical expressions).

We describe below the various estimators of the polarization fraction and angle listed in Table 1. We stress that most of the expressions derived in this work have been obtained when restricting the analysis in the 2D case, assuming furthermore that the true intensity I 0 is perfectly known, except for the Bayesian estimator where we present a 3D development (see Sect. 6).

Maximum Likelihood estimators

The Maximum Likelihood (ML) estimators are defined as the values of p 0 and ψ 0 which maximize the pdf calculated at the polarization measurements p and ψ. When computed using the 2D pdf f 2D to fit p 0 and ψ 0 simultaneously, this estimator gives back the measurements, whatever the bias and the covariance matrix are, and is inefficient at correcting the bias of the data.

After marginalization of the pdf f 2D over ψ, the 1D ML estimator of p 0 , pML , is now defined by

0 = ∂ f p ∂p 0 p | p 0 , ψ 0 , Σ p p 0 = pML . ( 6
)
Note that the expression of f p is independent of the measurement ψ, but still theoretically depends on the true value ψ 0 which is unknown. In the canonical case (ε eff =1) ψ 0 disappears from the expression, but it must be considered as a nuisance parameter in the general case. One way to proceed in such a case is to compute the mean of the solutions pML for ψ 0 varying in the range -π/2 to π/2. As already stressed by Simmons & Stewart (1985), this estimator yields a zero estimate below a certain threshold of the measurement p, which implies a strong discontinuity in the resulting distribution of this p 0 estimator. Nevertheless, contrary to the 2D ML estimators, the p ML estimator does not give back the initial measurements, and is often used to build polarization estimates.

Similarly, the 1D ML estimator of ψ 0 , ψML , is given after marginalization of f 2D over p by

0 = ∂ f ψ ∂ψ 0 ψ | p 0 , ψ 0 , Σ p ψ 0 = ψML . (7)
As mentioned for the ML estimator pML , the unknown parameter p 0 in the above expression has to be considered as a nuisance parameter when solving Eq. 7. We stress that because the canonical simplifications have always been assumed in the literature, bias on the ψ measurements has not been previously considered and the ψML estimator has not yet been used and qualified to correct this kind of bias. This analysis is done in Sect. 5.

Most Probable estimators

The Most Probable estimators of p 0 and ψ 0 are the values for which the pdf f 2D reaches its maximum at the measurements values (p,ψ). Notice that the Most Probable estimators ensure that the measurement values (p,ψ) are the most probable values of the pdf computed for this choice of p 0 and ψ 0 , i.e.they take the maximum probability among all possible measurements with this set of p 0 and ψ 0 . As a comparison the ML estimators ensure that the measurement values (p,ψ) take the maximum probability for this choice of p 0 and ψ 0 , compared to the probability of the same measurement values (p,ψ) for all other possible sets of p 0 and ψ 0 .

The 2D Most Probable estimators (MP2), pMP2 and ψMP2 , are defined as the values of p 0 and ψ 0 simultaneously satisfying the two following relations:

0 = ∂ f 2D ∂p p, ψ | p 0 , ψ 0 , Σ p p 0 = pMP2 ψ 0 = ψMP2 (8) and 0 = ∂ f 2D ∂ψ p, ψ | p 0 , ψ 0 , Σ p p 0 = pMP2 ψ 0 = ψMP2 . ( 9
)
These relations can be solved, using the fully developed expression of f 2D including the terms of the inverse matrix Σ -1 p , as provided in Appendix A. When canonical simplifications are assumed, this yields

ψMP2 = ψ , pMP2 = (p -σ 2 p /p) for p > σ p 0 for p ≤ σ p , (10)
as found in Quinn (2012). We observe that the MP2 estimate of the polarization fraction is systematically lower than the measurements, so that this estimator tends to over-correct p, as it will be shown in Sect. 4. After marginalization over p or ψ, the 1D Most Probable (MP) estimators, pMP and ψMP , are defined independently by:

0 = ∂ f p ∂p p | p 0 , ψ 0 , Σ p p 0 = pMP (11) and 0 = ∂ f ψ ∂ψ ψ |p 0 ψ 0 , Σ p ψ 0 = ψMP . (12)
The 1D and 2D estimators are not expected to provide the same estimates. Under the canonical assumptions, the MP estimator of p is commonly known as the Wardle and Kronberg's (Wardle & Kronberg 1974) estimator.

As mentioned earlier, the MP estimator yields a zero estimate below a certain threshold of p (Simmons & Stewart 1985), which implies a strong discontinuity in the resulting distribution of these estimators for low SNR measurements.

Asymptotic estimator

The Asymptotic estimator (AS) of the polarization fraction p is usually defined in the canonical case by

pAS = p 2 -σ 2 p p > σ p 0 p ≤ σ p . ( 13
)
The output distribution of the AS estimator appears as the asymptotic limit of the Rice (1945) distribution when p/σ p tends to ∞, just as the ML and MP estimators, and given by pdf

p σ p → N          p 0 σ p 2 + 1, 1          , (14)
where N(µ, σ) denotes the Gaussian distribution of mean µ and variance σ 2 . As with the previously presented estimators, this one suffers from a strong discontinuity at pAS =0.

In the general case, when the canonical simplification is not assumed, it has been shown by Plaszczynski et al. (2014, hereafter P14) that the expression of the asymptotic estimator can (green) and AS (red). We assume the covariance matrix to be canonical, and a SNR of p 0 /σ p =1.

be extended to a general expression by changing the term σ 2 p in Eq. 13 into a 'noise-bias' parameter b 2 defined by

b 2 = σ ′2 U cos 2 (2ψ 0 -θ) + σ ′2 Q sin 2 (2ψ 0 -θ) I 2 0 , ( 15
)
where θ represents the position angle of the iso-probability bivariate distribution, and

σ ′2 U , σ ′2 Q the rotated variances θ = 1 2 atan        2ρσ Q σ U σ 2 Q -σ 2 U        , (16)
σ ′2 Q = σ 2 Q cos 2 θ + σ 2 U sin 2 θ + ρσ Q σ U sin 2θ , (17)
σ ′2 U = σ 2 Q sin 2 θ + σ 2 U cos 2 θ -ρσ Q σ U sin 2θ . ( 18
)
and ψ 0 is the true polarization angle, which can be approximated asymptotically by the naïve measurement ψ or, even better, by the estimate ψML of Sect. 2.1. It has been shown that this equivalent 'noise-bias' b 2 ensures the minimal bias of pAS .

Discontinuous estimators

The estimators of p introduced above (ML, MP and AS) exhibit a common feature: below some cutoff value the estimator yields exactly zero. This means that the estimator distribution is discontinuous and is a mixture of a discrete one (at p=0) and a continuous one (for p > 0), This type of distribution is illustrated in Fig. 1 for a SNR p 0 /σ p =1 and a canonical covariance matrix. The distribution of the naïve measurements is built using a Monte-Carlo simulation, starting from true polarization parameters p 0 and ψ 0 . The other three distributions of p are obtained after applying the ML, MP and AS estimators. A non negligible fraction of the measurements provide null estimates of p. As shown in Fig. 2, this fraction of null estimates reaches 40% at low SNR with the MP and AS estimators, and more than 50% with the ML estimator for SNR<1. It converges to 0% for SNR >4.

If taken into account as reliable estimate of p, null estimates will somewhat artificially lower the statistical bias of the p estimates compared to the true value p 0 , as detailed in Sect. 4. A null value of these estimators should be understood as an indicator of the low SNR of this measurement, which has in fact to be included into any further analysis as an upper limit value. In practice, the user seldom has various realizations at hand. Using these estimators then leads to a result with upper limits mixed with non-zero estimates in the analysis. Such complications may be especially hard to handle when studying polarized maps of the interstellar medium. On the other hand, it would be disastrous to omit those estimates in any statistical analysis, since weaklypolarized points would be systematically rejected. To avoid such complications, we explore below other estimators which avoid this issue and lead to continuous distributions. This is especially important in the range of SNR between 2 and 3, where the discontinuous estimators still yield up to 20% of null estimates.

Modified Asymptotic estimator

A novel powerful asymptotic estimator has been introduced by Plaszczynski et al. (2014) to correct the discontinuous distribution of the standard estimators while still keeping the asymptotic properties. It has been derived from the first order development of the Asymptotic estimator, which has been modified to ensure positivity, smoothness and asymptotical convergence at high SNR. The Modified Asymptotic (MAS) estimator is defined as follows:

pMAS = p -b 2 • 1 -e -p 2 /b 2 2p , (19)
where the 'noise-bias' b 2 is given by Eq. 15 and computed using a polarization angle assessed from each sample using the asymptotic estimator ψ.

P14 also provides a sample estimate of the variance of the estimator that is shown to represent asymptotically the absolute risk function (defined in Sec. 3.1) of the estimator:

σ 2 p,MAS = σ ′2 Q cos 2 (2ψ -θ) + σ ′2 U sin 2 (2ψ -θ) I 2 0 . ( 20
)
This estimator focuses on getting a "good" distribution, that transforms smoothly from a Rayleigh-like to a Gaussian one, the latter being reached in the canonical case for an SNR of about 2.

Bayesian estimators

The pdfs introduced in Sect. 2 provide the probability to observe a set of polarization measurements (I, p, ψ) given the true po-larization parameters (I 0 , p 0 , ψ 0 ) and the covariance matrix Σ. Because we are interested in the opposite, i.e. getting an estimate of the true polarization parameters given a measurement and the knowledge of the noise properties, we use Bayes Theorem to build the posterior distribution. The posterior pdf B is given in the 3D case by

B(I 0 , p 0 , ψ 0 | I, p, ψ, Σ) = (21) f (I, p, ψ | I 0 , p 0 , ψ 0 , Σ) • κ(I 0 , p 0 , ψ 0 ) +∞ 0 1 0 π/2 -π/2 f (I, p, ψ | I ′ 0 , p ′ 0 , ψ ′ 0 , Σ) κ(I ′ 0 , p ′ 0 , ψ ′ 0 ) dI ′ 0 d p ′ 0 dψ ′ 0 ,
where κ(I 0 , p 0 , ψ 0 ) is the prior distribution, which represents the a priori knowledge of the true polarization parameters and has to be positive everywhere and normalized to 1 over the definition ranges of I 0 , p 0 and ψ 0 . When no a priori knowledge is provided, we have to properly define a 'flat' prior, or noninformative prior, which encodes the ignorance of the prior. A class of non-informative priors is given by the Jeffreys' prior [START_REF] Jeffrey | Theory of probability[END_REF] where the ignorance is defined under symmetry transformations that leave the prior invariant. As discussed by Quinn (2012) for the two dimensional case, this kind of prior can be built as a uniform prior in cartesian space (Q 0 ,U 0 ) or in polar space (p 0 , ψ 0 ), both expressing the ignorance of location.

We will prefer the latter, uniform in polar space, which ensures uniform sampling even for small values of p 0 . While p 0 and ψ 0 are only defined on a finite range ([0, 1] and [-π/2, π/2], respectively), the intensity I 0 may be infinite in theory, which leads to an issue when defining the ignorance prior. In practice, an approximation of the ignorance prior for I 0 will be chosen as a tophat centered on the measurement I and chosen to be sufficiently wide to cover the wings of the distribution until it becomes negligible. Such uniform priors lead to the expression of B given in Eq. 23, where the normalization factor has been omitted. We emphasize that the definition of the ignorance prior introduced above becomes data-dependent, which is not strictly following the Bayesian approach. Furthermore, the question of the definition range of the prior and the introduction of non-flat priors will be discussed in Sect. 4.3, in the context of comparing the performance of the estimators inspired by the Bayesian approach.

Similarly, the posterior pdf in 2D (i.e., when the total intensity is perfectly know, I = I 0 ) is defined by

B 2D (p 0 , ψ 0 | p, ψ, Σ p ) = f 2D (p, ψ | p 0 , ψ 0 , Σ p ) • κ(p 0 , ψ 0 ) 1 0 +π/2 -π/2 f 2D (p, ψ|p ′ 0 , ψ ′ 0 , Σ p ) κ(p ′ 0 , ψ ′ 0 ) d p ′ 0 dψ ′ 0 .
(24) The analytical expressions of the posterior pdf B 2D with a flat prior is given in Eq. 25, where the normalization factors have been omitted and the intensity has been assumed perfectly known (I = I 0 ). Illustrations of this posterior pdf are presented in Appendix B. We also introduce B p and B ψ the Bayesian posterior pdfs of p and ψ in 1D, respectively, and defined as the marginalization of the B 2D over ψ and p, respectively.

We use the Bayesian posterior pdf in 2D B 2D to build two frequentist estimators: the MAP and the MB.

The MAP2 and MAP estimators in 2D and 1D, respectively, are simply defined as the (p 0 , ψ 0 ) values corresponding to the maximum of the posterior pdf, B 2D , and B p and B ψ , respectively. We recall that these estimators match exactly the ML estimators of Sect. 2 in one and two dimensions, respectively, when a flat prior is assumed. Hence the MAP2 estimators yield back the polarization measurements, whereas the MAP estimators provide a simple way to compute the ML estimates.

The Mean Bayesian Posterior (MB) estimators are defined as the first order moments of the posterior pdf:

pMB ≡ +π/2 -π/2 1 0 p 0 B 2D (p 0 , ψ 0 | p, ψ, Σ p )d p 0 dψ 0 (26) and ψMB ≡ ψ+π/2 ψ-π/2 1 0 ψ 0 B 2D (p 0 , ψ 0 | p, ψ, Σ p )d p 0 dψ 0 . ( 27
)
Notice that in the definition of ψMB the integral over ψ 0 is performed over a range centered on the measurement ψ. This has to be done to take into account the circularity of the posterior pdf over the

ψ 0 dimension. Note that B 2D (p 0 , ψ 0 | p, ψ, Σ p ) = B 2D (p 0 , ψ 0 + π | p, ψ, Σ p ).
We stress that the frequentist estimators inspired by a Bayesian approach, pMB and ψMB , introduced above in the 2D case can be easily extended to the 3D case by integrating the pdf B(I 0 , p 0 , ψ 0 | I, p, ψ, Σ) of Eq. 21 over the I, p and ψ dimensions. This is extremely powerful when the uncertainty of the intensity I has to be taken into account in the estimate of the polarization parameters, which is highly recommended in some circumstances, such as a low SNR on I (<5) or the presence of an unpolarized component on the line of sight (see Sect. 6 and PMA I for more details).

Uncertainties

Variance and risk function

It is important not to confuse the variance (noted V) of an estimator with its absolute risk function (noted R). For any distribution of the random p variable the definitions are :

V ≡ E (X -E[X]) 2 , (28) R ≡ E (X -X 0 ) 2 , ( 29
)
where E[X] is the expectation of the random variable X and X 0 is the true value. Introducing the absolute bias in E[X] = X 0 + B and expanding both relations, the link between the variance and the absolute risk function is simply:

V = R -B 2 . ( 30
)
Therefore, for a constant absolute risk function, the variance decreases with the absolute bias and both are equal when the estimator is unbiased. The variance does not require knowing

B(I 0 , p 0 , ψ 0 | I, p, ψ, Σ) ∝ Det(Σ -1 ) 2(π) 3 exp            - 1 2          I -I 0 p I cos(2ψ) -p 0 I 0 cos(2ψ 0 ) p I sin(2ψ) -p 0 I 0 sin(2ψ 0 )          T Σ -1          I -I 0 p I cos(2ψ) -p 0 I 0 cos(2ψ 0 ) p I sin(2ψ) -p 0 I 0 sin(2ψ 0 )                     , ( 23
)
B 2D (p 0 , ψ 0 | p, ψ, Σ p ) ∝ 1 πσ 2 p,G exp         - 1 2 p cos(2ψ) -p 0 cos(2ψ 0 ) p sin(2ψ) -p 0 sin(2ψ 0 ) T Σ -1 p p cos(2ψ) -p 0 cos(2ψ 0 ) p sin(2ψ) -p 0 sin(2ψ 0 )         , (25)
the true value of the random variable, which makes it useful to provide an uncertainty estimate, but it has to be used extremely carefully in the presence of bias. In such cases, the variance will always underestimate the uncertainty. Furthermore, it is known that the variance is not appropriate for providing uncertainties with non-Gaussian distributions, which is the case for the polarization fraction and angle. In such circumstances, the confidence intervals (see Sect. 3.3) are the preferred method for obtaining robust uncertainties. The variance, however, is often used as a proxy of the uncertainty in the high regime of the SNR. We will detail in Sect. 4.5 and 5.3 in which conditions this can still be applied.

Posterior uncertainties

One of the main benefits of the Bayesian approach is to provide simple estimates of the uncertainties associated with the polarization estimates. One option is to build credible intervals around the MAP estimates as first proposed by Vaillancourt (2006), and the other option is to use the variance of the pdf.

Given a polarization measurement (p, ψ) and the posterior pdf B 2D (p 0 , ψ 0 |p, ψ, Σ p ), the lower and upper limits of the λ% credible intervals are defined as the lower and upper limits of p 0 and ψ 0 ranging the iso-probability region Ω(λ, p, ψ) over which the integral of B equals λ%, so that

Ω(λ,p,ψ) B 2D (p 0 , ψ 0 | p, ψ, Σ p ) d p 0 dψ 0 = λ 100 . ( 31
) These intervals, [p low MAP2 , p up MAP2 ] and [ψ low MAP2 , ψ up MAP2 ]
, estimated from the 2D expression of B 2D are defined around the MAP2 estimates pMAP2 and ψMAP2 , which are equal to the measurements (p, ψ).

A similar definition can be given in the one-dimensional case, which leads to different results. The lower and upper limits, p low MAP and p up MAP , around pMAP are defined as follows

p up MAP p low MAP B p (p 0 | p, Σ p ) d p 0 = λ 100 , (32)
with the constraint that the posterior probability function is identical for p low MAP and p up MAP . Similarly, the lower and upper limits, ψ low MAP and ψ up MAP , around ψMAP are given by

ψ up MAP ψ low MAP B ψ (ψ 0 | ψ, Σ p ) dψ 0 = λ 100 . ( 33
)
We recall that this integral has to be computed around the measurement value ψMAP to take into account the circularity of the posterior pdf with the polarization angle. Notice that the credible intervals built in 1D or 2D are not supposed to be identical, as ( pMAP2 , ψMAP2 ) and ( pMAP , ψMAP ) are not equal in the general case.

The second definition of the uncertainty comes from the second moment of the 1D posterior probability density functions B p and B ψ , as follows:

σ 2 p,MB ≡ 1 0 (p 0 -pMB ) 2 B p (p 0 | p, Σ p ) d p 0 , (34)
and

σ 2 ψ,MB ≡ ψ+π/2 ψ-π/2 (ψ 0 -ψMB ) 2 B ψ (ψ 0 | ψ, Σ p ) dψ 0 . ( 35
)
The operation of subtraction between the two polarization angles must be done with care, restricting the the maximum distance to π/2. At very low SNR, i.e. an almost flat uniform pdf, the uncertainty reaches the upper limit σ ψ,MB ≤ π/ √ 12 rad = 51. • 96. We stress that these 1-σ estimates may not be associated with the usual 68% confidence intervals of the normal distribution, because of the asymmetrical shape of the posterior distribution and because of the circularity of the angular variable.

Confidence intervals

So far we have considered point estimation of the true p 0 value which is somewhat tricky in the low SNR regime because of the non-Gaussian nature of the estimator distribution.

A different approach, that takes into account the entire shape of the distribution is to build confidence regions (or intervals), which allows at some given significance level, to obtain bounds on the true value given some estimator value. Simmons & Stewart (1985) have built the so-called Neyman "confidence belt" for the naïve estimator in the canonical case. PMA I proposed the construction of two-dimensional (p 0 , ψ 0 ) intervals, for the general covariance matrix case. The classical construction suffers from a standard issue: at very low SNR the confidence interval lies entirely in the unphysical p < 0 region, and both previous studies provide over-conservative regions.

P14 has implemented the Feldman-Cousins prescription (Feldman & Cousins 1998) which is based on using a likelihood ratio criterium in the Neyman construction. This allows building intervals that always lie in the physical region without ever being conservative. They provided these intervals for the MAS estimator including analytical approximations to the upper and lower limits for 68, 95 and 99.5% significance levels.

p estimator performance

Methodology

We investigate in this section the capability at providing polarization fraction estimates with low bias of the seven p estimators introduced in the previous sections: the naïve measurements p, the Maximum Likelihood (ML), the Most Probable (MP and MP2), the Asymptotic (AS), the Modified Asymptotic (MAS) and the Mean Bayesian Posterior (MB) estimators. Their performance is first quantified in terms of relative bias and risk function of the resulting estimates. Given true polarization parameters (p 0 , ψ 0 ) and a covariance matrix Σ p , we build a sample of one million simulated measurements (p,ψ) by adding noise on the true Stokes parameters using the covariance matrix. We define the relative bias and risk function on p as follows:

Bias p ≡ p -p 0 σ p,G and Risk p ≡ ( p -p 0 ) 2 σ 2 p,G , ( 36
)
where p is the polarization fraction estimate computed on the simulated measurements p, p 0 is the true polarization fraction, <> denotes the average computed over the simulated sample, and σ p,G is the estimate of the noise of the polarization fraction. The choice of σ p,G to scale the absolute bias and risk function, as a proxy of the p uncertainty, is motivated by the fact that it depends only on the effective ellipticity and not on ψ 0 . Notice that this choice can lead to a relative risk function falling below 1 at low SNR, due to the fact that σ 2 p,G >Var in this regime. The accuracy of the p estimators is also quantified regarding the , risk function (middle) and Jarque-Bera test (bottom) of the pure measurements (naïve, black), ML (dashed blue), MP (dashed light green), MP2 (dashed green), AS (dashed red), MAS (orange) and MB (pink) p estimators in the canonical case, as a function of the the SNR p 0 /σ p . The dashed lines stand for the discontinuous estimators presenting a peak of their output distribution at p=0.

shape of their output distributions. We use the Jarque and the MB (pink) p estimators in the canonical case (ε eff =1), for three levels of the SNR p 0 /σ p =1,2 and 5 (from top to bottom).

In the case of the MB estimator, we show two setups of p 0 =1% and 50% to illustrate the dependence of the output distribution on the p 0 value, due to the prior used in the Bayesian approach ( pMB ∈ [0, 1] so that pMB /p 0 ∈ [0, 1/p 0 ]). The other estimators are not sensitive to the true value p 0 .

distribution, and defined by

JB = n 6        µ 2 3 µ 3 2 +       µ 4 µ 2 2 -3       2 /4        , ( 37
)
where n is the number of samples and µ i is the naïve estimate of the ith central moment of the distribution. This test is based on the joint hypothesis of the skewness and the excess kurtosis being zero simultaneously. A value JB=0 means a perfect agreement with the normality to the 4th order, but does not prevent departure from the normality at higher orders. This JB estimator tends to a χ 2 test with two degrees of freedom when n becomes large enough. Hence the JB has to satisfy the following condition JB < χ 2 α , once chosen a significance level α. For a significance level α=5% and 1%, we get the conditions JB < 5.99 and JB < 9.21, respectively.

Canonical case

We first assume the canonical simplification of the covariance matrix (ε eff =1). The relative Bias p and Risk p quantities are shown on Fig. 3 for the seven p estimators. We recall that the discontinuous estimators, shown in dashed line (ML (blue), MP (light green), MP2 (green) and AS (red)), have an output distribution presenting a strong peak at zero, which leads to artificially lower the statistical relative Bias p when simply including null values instead of using upper limits, as discussed in Sect. 2.4. Effectively these estimators show the lowest relative biases (top panel of Fig. 3) compared to the MAS (orange) and MB (pink) estimators. Hence the ML and MP2 estimators seem to statistically over-correct the data, below SNR=3. Consequently, the ML, MP and AS p estimators have to be used with an extreme care to deal with null estimates. We suggest here to focus on the two continuous estimators: MAS and MB.

MAS provides the better performances in terms of relative bias over the whole range of SNR, while MB appears less and less efficient at correcting the bias when the SNR tends to zero. At larger SNR (>2), MB tends to slightly over-correct with a small negative relative bias (2% of σ p ) up to SNR ∼5, while MAS converges quickly to a null relative bias for SNR > 3.

The MB estimator clearly minimizes the risk function (in the range 0.7<SNR<3.2), as expected for this kind of posterior estimator. At larger SNR (>3.2) both MAS and MB have roughly the same behavior, even if the risk function associated to MAS appears slightly lower.

The resulting pMB distribution is highly asymmetric at low SNR (see upper panels of Fig. 4), with a sharp cutoff at 0.8σ p . Moreover, we note that the output pMB distribution depends not only on the SNR p 0 /σ p , but also on the value of the true polarization fraction p 0 . We report two cases, p 0 =1% (pink) and 50% (dotted pink) in Fig. 4. This comes from the prior of the Bayesian method, which bounds the estimate pMB between 0 and 1. As a consequence, the normality of the Bayesian distribution is extremely poor, as pointed in the bottom panel of Fig. 3, where we show that the JB test of the MB estimator is larger than 9.21 (consistent with a χ 2 0.01 test) over the whole range of SNR explored here (up to SNR∼5). On the contrary, the resulting pMAS distribution of Fig. 4 looks much better, mimicking the Rayleigh distribution for low SNR and going neatly to the Gaussian regime, as pointed out by P14. The JB of the MAS estimator is the lowest for SNR >3 (see bottom panel of Fig. 3), illustrating the consistency between the MAS distribution and the normal distribution. Notice that all distributions, naïve, MAS and MB, converge to a Gaussian distribution at higher SNR.

κ ∈ [0, 100p 0 ] κ ∈ [0, 10p 0 ] κ ∈ [0, 5p 0 ] κ ∈ [0, 3p 0 ] κ ∈ [0, 2p 0 ]

Impact of the Bayesian prior

The choice of the prior is crucial in the Bayesian approach, and we have seen how it is hard to define a non-informative prior in Sect. 2.6. The MB estimator studied up to now assumes a flat prior in p 0 between 0 and 1, equivalent to no a priori knowledge. In practice when dealing with astrophysical data, we can bound the expected true values of the polarization fraction between much tighter limits. We know, for example, that the polarization fraction of the synchrotron signal peaks at ∼75%, but never reaches this maximum due to line-of-sight averaging. The maximum polarization fraction of the dust thermal emission is still a debated issue, but is unlikely to be larger than 20 to 30% (Benoît et al. 2004). Appropriate priors can then be introduced to take into account this a priori physical knowledge into the MB estimator.

We have already observed in Sect. 4.2 how the output distribution of the pMB estimates is impacted by the value of the true p 0 (1% or 50%) due to the upper limit (p 0 <1) of the prior, see Fig. 4. We explore here a family of simple priors defined by κ(p ′ 0 ) = 1/(kp 0 ) for p ′ 0 ∈ [0, kp 0 ] and 0 otherwise, where we ad- just the upper limit of the prior as a function of the expected true value. We performed Monte Carlo simulations in the canonical case by setting the true value at p 0 =1% and varying the upper limit of the prior (k = 2, 3, 5, 10, and 100). The statistical relative Bias p of the MB estimators associated with each version of the priors are shown on Fig. 5. The smaller the upper limit, the lower the relative Bias p , as expected. However the upper limit of the prior has to be very constraining (k ≤ 3) to observe a decrease of the relative bias in the range of SNR between 1.5 and 3. This requires very good a priori knowledge. Using more relaxed priors (k ≥ 5) will not improve significantly the performances of the MB estimator at SNR>1. When dealing with maps of polarized data, an interesting approach would be to start by estimating the histogram of p values in the map and use it as a prior into our MB estimators, even if this moves away from a strictly Bayesian approach again by introducing a data-dependent prior. As a first guess, the prior can be set to the histogram of the naïve estimates of p, but a more sophisticated prior would be an histogram of p deconvolved from the errors, using a Maximum Entropy method for example.

We illustrate the performance of the MB estimator with this kind of prior on Figs. 6 and7. We start with a sample of 10 000 independent true values (p 0,i ) ranging between 0 and 20% of polarization fraction, with a distribution shown as the grey shaded histogram in Fig. 7 on which a random realization of the noise is added with the same noise level over the whole sample, leading to varying SNRs through the sample. We explore two extreme cases of the Bayesian prior, corresponding to i) an idealistic perfect knowledge of the input distribution and ii) its first guess provided by the naïve estimates. Hence the prior is chosen as the input distribution of the true p 0,i values (dashed pink) and the output distribution of the naïve estimates (dotted pink). We compare the performance of these two new versions of the MB estimators with the naïve (black), MAS (orange) and flat prior MB (solid pink) estimators, in terms of relative bias in Fig. 6.

We stress that the relative bias values are not defined as previously done in Sect. 4.1, but refer now to the mean of the difference between each sample of true value p 0,i and its associated estimate pi . The pink shaded region provides the domain of the possible improvement of the MB estimators, by setting an appropriate prior as close as possible to the true distribution. The improvements may seem spectacular, leading to a statistical relative bias close to zero at all SNRs in the best configuration (dashed line). Caution is warranted, however, when looking at the output distributions associated with these new MB estimators on Fig. 7, shown for three levels of the noise chosen so that the mean SNR is p 0 /σ p,G =1, 2 and 3. At low SNR (≃1), the output distribution of the MB estimator with a perfect prior (dashed line) is extremely peaked around the mean value of the sample p 0 , but does not match the input distribution at all. Even at higher SNR (2-3), the three MB output distributions suffer from the same feature already mentioned in Sect. 4.2, a sharp cutoff at low values of p. Using a prior that is too constraining will yield dramatic cuts of extremes values of the input distribution. By contrast, the naïve prior is quite effective in that it allows the MB estimator to recover the upper limit of the input distribution reasonably well at a SNR 2, while the other estimators fail to do so at such low SNR.

The performance of the MB estimator with an evolved prior will also strongly depend on the initial true distribution of the polarization fraction. For example we duplicated the analysis made above with a different initial distribution (p 0,i ) centered on 20% of polarization fraction instead of 10% (see Fig. 8). In this configuration, the output distributions of the Bayesian estimators are not as much affected by the cut-off at low p as observed in Fig 7 .  The MB estimator with the naïve prior appears extremely effective, even at low SNR (∼2).

Robustness to the covariance matrix

In PMA I we have discussed extensively the impact of the asymmetry of the covariance matrix on the measurements of the polarization fraction. In particular, we have stressed that once the effective ellipticity departs from the canonical case, the bias on the polarization fraction now depends on the true polarization angle ψ 0 , which remains unknown. We would like to explore in this section how the performance of the various p estimators are sensitive to the effective ellipticity of the covariance matrix.

We illustrate the dependence of the p estimators on the true polarization angle ψ 0 in Fig. 9. Given true polarization parameters (p 0 =0.1 and ψ 0 ranging between -π/2 and π/2) and a covariance matrix characterized by ε eff =2 and θ=0 (left panel), and a SNR p 0 /σ p,G =1, we first set the polarization measurements (p, ψ) to the maximum of the pdf f 2D (left panel). We apply then the six estimators on these measurements to get the p estimates for each ψ 0 between -π/2 and π/2. With this particular setting, the MP2 (green) estimator gives back the true polarization fraction p 0 whatever the polarization angle ψ 0 , by definition of this estimator and the choice of the measurement in this example. On the contrary, the MP (light green) and the ML (blue) estimators are extremely sensitive to the true polarization angle ψ 0 , yielding estimates spanning between 0 and 1.4p 0 , while the AS (red) and MAS (orange) estimators yield results spanning between 1 to 1.8p 0 when ψ 0 varies. The MB (pink) estimator provides stable estimates in the range 1.4 to 1.5 p 0 , which is consistent with the fact that the posterior estimators minimize the risk function. This of course has a cost, and the MB estimator provides here the largest averaged relative bias compared to the other methods, with the exception of the naïve (black) one.

More generally, for each value of the true polarization angle ψ 0 between -π/2 and π/2, we build a sample of 10 000 simulated measurements using the same setup of the covariance matrix as above. Then we compute the statistical average of the naïve, MAS and MB estimates (black, orange and pink lines, respectively) obtained on this simulated sample, with their associated 1-σ dispersion (black, orange and pink dot-dash lines, respectively), as shown in the right panel of Fig. 9. The averaged MB estimates present the same characteristic as shown on the left panel. By contrast, the averaged MAS estimates are independent from the unknown ψ 0 true polarization angle. The MAS 1-σ dispersion is, however, slightly larger than the MB 1-σ dispersion.

The impact of the effective ellipticity of the covariance matrix is then analysed statistically for the MAS and MB estimators only in Fig. 10. Instead of looking at the accuracy of the p estimators around one particular measurement (the most probable one) as done in Fig. 9, for each set of true polarization parameters (p 0 =0.1, ψ 0 ), with ψ 0 ranging between -π/2 and π/2, we perform Monte Carlo simulations. For each set of true polarization parameters, we build a sample of 100 000 simulated measurements on which we apply the MAS and MB estimators to finally compute the statistical relative Bias p and Risk p , as defined in Sect. 4.1. This is done for various setups of the covariance matrix chosen to cover the whole range of the extreme and low regimes. The minimum and maximum relative Bias p and Risk p are then computed over the whole range of ψ 0 and effective ellipticity ε eff in each regime of the covariance matrix to build the shaded regions of Fig. 10

ψ 0 [rad] < MB > < MAS > < Naive > 1-σ MB 1-σ MAS 1-σ Naive p 0 true p 0 /σ p,G =1 ε eff =2 θ=0
p Fig. 9. Illustration of the robustness of the p estimators against the unknown ψ 0 parameter when the covariance matrix departs from the canonical value. The covariance matrix is setup with ε eff =2 and a SNR p 0 /σ p,G = 1, and a true polarization fraction p 0 =0.1. For each value of ψ 0 , we first illustrate (on the left panel) the performance of the 7 estimators on one particular measurement set to the maximum of the pdf. We focus then on the statistical average estimates p computed over 10 000 Monte-Carlo realizations for the naïve, MAS and MB estimators (right panel). The full lines stand for the mean, and the dot-dash lines for the 1-σ dispersion. of the MAS estimator is less impacted by a change of ellipticity for SNR>2 than the MB estimator, even in the extreme regime of the covariance matrix. The dependance of the risk function on the ellipticity is almost identical for the two estimators around their respective canonical curve. The thickness of the risk function region is slightly smaller for the MB estimator than for the MAS estimator at low SNR (<3), while it is the opposite for larger SNR (>3), as already observed in the canonical case.

ε eff =1 1 < ε eff < 1.1 1 < ε eff < 2 MB (ε eff =1) 0 1 2 3 4 5 p 0 /σ p,G 0 1 2 3 Risk p MAS ε eff =1 1 < ε eff < 1.1 1 < ε eff < 2 Naive (ε eff =1) MB (ε eff =1) 0 1 2 3 4 5 p 0 /σ p,G -0.5 0.0 0.5 1.0 1.5 Bias p Naive MB ε eff =1 1 < ε eff < 1.1 1 < ε eff < 2 MAS (ε eff =1) 0 1 2 3 4 5 p 0 /σ p,G 0 1 2 3 Risk p MB ε eff =1 1 < ε eff < 1.1 1 < ε eff < 2 Naive (ε eff =1) MAS (ε eff =1)

Polarization fraction uncertainty estimates

The questions of estimating the polarization uncertainties and how uncertainties are propagated are essential in reliable polarization analysis. The best approach consists of building the confidence intervals to retrieve robust estimates of the lower and upper limits of the 68, 95 or 99.5% intervals, which is valid even when the distribution is not Gaussian. As already mentioned in sect. 3.3, building optimized confidence intervals including the full knowledge of the covariance matrix may represent a challenge for large samples of data. Hence P14 provides analytic approximations of such confidence intervals for the MAS estimator, which can be extremely useful.

A commonly used approach, however, is to provide the 1-σ dispersion, assuming the Gaussian distribution of the p estimates as a first approximation. We have already stressed the difference between the risk function and the variance, and the limitations of the latter to derive robust uncertainties in the presence of bias. We compare below the performance of the usual uncertainty estimates introduced in Sect. 3 to provide robust 68% tolerance intervals: MAS variance, credible intervals MAP and 1-σ a posteriori dispersion MB.

Starting with a true p 0 value, we have performed Monte-Carlo simulations in the low regime of the covariance matrix, by exploring the whole range of the true polarization angle ψ 0 , with an SNR spanning from 0 to 30. For each simulated measurement (p,ψ), we compute the p estimates with their uncertainty estimators σ p. We then compute the a posteriori probability to find the true p 0 inside the interval [ pσ low p , p + σ up p ]. In the case of the MAP estimator, the lower and upper limits of the interval, pMAPσ low pMAP and pMAP + σ up pMAP , are set to p low MAP and p up MAP , respectively, (with λ=68 as defined in Sect. 3.2), which can be asymmetric. We report the results compared to the expected 68% level in Fig. 11. We recall that this comparison approach is frequentist again, while anything derived from the Bayesian pdf is used to build single estimates and to be compared with the confidence intervals.

As pointed out in Sect. 3.1, the theoretical variance associated with the MAS estimator still tends to provide slightly lower probabilities than the expected 68% at low SNR, mainly due to the asymmetry of the distribution. The variance associated to the MB estimator, which is more biased at low SNR, gives extremely low probability to recover the true p 0 value at low SNR (< 0.5). By contrast, it provides probabilities greater than 68% (as high as 90%) for SNR between 0.5 and 2. This comes from the fact that the MB variance statistically over-estimates by a factor of 2 the exact variance of the a posteriori pMB distribution at low SNR (<2). Thus the MB uncertainty estimator yields conservative estimates of the uncertainty for SNR >0.5. At high SNR (>3) all these uncertainty estimators provide compatible estimates of the probability close to 68%.

Because the true SNR is always unknown (see Sect. 4.6), the probability to find the true p 0 value in the confidence interval is also shown as a function of the measured SNR in Fig. 12. This much more realistic picture shows that the variance estimates provide reliable probability for measured SNR larger than ∼6.

Polarization signal-to-noise ratio

In any real measurement, the true SNR p 0 /σ p,G remains unknown. From observations, we only have access to the measured SNR, which can be obtained by the ratio p/σ p associated with each estimator, or by a confidence interval approach (see P14), which is much more robust at a low true SNR. We show in Fig. 13 the accuracy of the measured SNR compared to the true SNR for the four following methods: the naïve estimate plus Classical estimate of the uncertainty, the MAS estimate with the associated variance, the MB estimate and its variance, and the ML estimate with the MAP credible intervals. We observe that all methods agree only for a true SNR larger than 3, giving back the true SNR in this regime. Below this true SNR, the measured SNR becomes extremely biased whatever the method used, due to the bias of the measurement p itself, but also due to the bias introduced by the variance as an estimate of the uncertainty when the output distribution departs from the Gaussian regime.

ψ estimator performance

Methodology

As pointed out by PMA I, once the covariance matrix is not canonical (ε eft > 1), a bias of the polarization angle measurements ψ appears with respect to the true polarization angle ψ 0 . This bias may be positive or negative. We propose to compare the accuracy at correcting the bias of the polarization angle of the four following ψ estimators: naïve measurements ψ, the ML ψML (which is equivalent to the MAP ψMAP ), the MP2 ψMP2 and the MB ψMB .

Similarly to the p estimators, we define the relative bias and risk function on ψ as follows:

Bias ψ ≡ ψ -ψ 0 σ ψ,0 and Risk ψ ≡ ( ψ -ψ 0 ) 2 σ 2 ψ,0 , (38)
where ψ is the polarization angle estimate computed on the simulated measurements ψ, ψ 0 is the true polarization fraction and angle, <> denotes the average computed over the simulated sample, and σ ψ,0 is the standard deviation of the simulated measurements.

Performance Comparison

We explore the performance of the four ψ estimators at four SNR=0.5, 1, 2 and 5 (from top to bottom) and a covariance matrix with an effective ellipticity ε eff =2, on Fig. 14. The relative Bias ψ (left panels) and Risk ψ (right panels) are plotted as a function of the true polarization angle ψ 0 . While the MB (pink) estimator seems to provide the least biased estimates with the lowest risk function at low SNR (<1), it becomes the least efficient at higher SNR. On the contrary, the ML (or MAP too) presents poor performances at low SNR, but provides impressive results at high SNR, reducing the relative bias close to zero at a SNR of 5. The MP2 estimator does not present any satisfactory properties: strong relative bias and risk function in almost all cases. Hence this ψMP2 estimator can be ruled out. An overview of the performance of the four ψ estimators as a function of the SNR is shown on Fig. 15, after marginalization over all the possible values of the ψ 0 parameter. As the relative Bias ψ can be positive or negative depending on ψ 0 , we compute the average of the absolute value of the relative bias, < |Bias ψ | > as an indicator of the statistical performance of the estimators whatever the true polarization angle is. We observe again on the left panel of Fig. 15 that the MB (pink) estimator provides the lowest relative bias for SNR<1.2, while the ML is especially powerful for SNR>2. All estimators provide almost the same results for the average Risk ψ (left panel), even if MB appears slightly better than the others, including the naïve measurements.

The examples provided above have been computed with an extreme effective ellipticity (ε eff =2) to emphasize the observations, but the same conclusions can be reached for lower values of the ellipticity. See, for example, the case with ε eff =1.1 shown in dotted line in Fig. 15. In the low regime of the covariance matrix, however, the statistical relative bias on ψ is very small, typically smaller than 5% of the dispersion, so that the need to correct the bias on ψ remains extremely limited.

Polarization angle uncertainty estimates

Once a reliable estimate of ψ based on the MB and ML (MAP) estimators has been obtained, we would like to build a robust estimate of the associated uncertainties σ ψ, which should be done by building confidence intervals. Because this last step could represent important efforts in some cases, for example when dealing with the full covariance matrix, we detail other methods below.

One option is to use the uncertainty associated with the MB estimator, σ ψ,MB (see Eq. 35). Another is to use the credible intervals built around the MAP estimates on the posterior pdf. We can keep the lower and upper limits, ψ low MAP and ψ up MAP computed for a 68% credible interval, or build a symmetrized uncertainty:

σ ψ,MAP = 1 2 ψ up MAP -ψ low MAP . (39)
A third option consists in taking the classical uncertainty given in PMA I, derived from the derivatives of the polarization parameters. PMA I has already shown that this ψ uncertainty estimator, associated with the naïve measurements, tends to systematically underestimate the true dispersion of the ψ distribution. We first assume the canonical simplification of the covariance matrix, which implies that the ψ measurements are not statistically biased. We also recall that under such assumptions the ML (MAP) and MB ψ estimators will give back the measurements ψ. We study, however, how the uncertainties associated with these two estimators can be used to get a reliable estimate of the uncertainty σ ψ. Starting from a true (p 0 , ψ 0 ), we simulate a sample of 50 000 simulated measurements p, ψ at a given SNR p 0 /σ p , on which we apply the two ML (MAP) and MB ψ estimators and their associated uncertainty σ ψ,MAP and σ ψ,MB , respectively. From this simulated set we can derive the averaged σ ψ for both methods. Because all estimators give back the measurements in the canonical case, we compare the MAP (blue) and MB (pink) polarization angle uncertainties estimators directly to the true dispersion (black) of the ψ measurements in Fig. 16. We also repeat the average of the classical estimates (dashed line) of the polarization uncertainty estimate, which has been shown by PMA I (see their Fig. 7) to underestimate by a factor of two the true uncertainty at low SNR (<2). We observe that the MAP estimator σ ψ,MAP provides an extremely good estimate of the polarization angle uncertainty compared to the true one over the whole range of SNR, even if slightly conservative up to a SNR of 5. The MB estimator σ ψ,MB provides consistent estimates of the uncertainty from intermediate SNR∼1, but still underestimates at lower SNR (<1).

In the non-canonical case a statistical bias on ψ appears, which can be partially corrected using the appropriate ψ estimators (see Sect. 5.2), leading to an output distribution of the ψ estimates. We quantify the performance of the ψ uncertainty estimators via Monte-Carlo simulations, as done for the p uncertainties. Starting from a set of polarization parameters (p 0 =0.1, -π/2<ψ 0 <π/2), we build a sample of simulated measurements (p, ψ) using various setups of the covariance matrix in the low regime, and various SNRs ranging from 0 to 30. We then compute the a posteriori probability to find the true polarization angle ψ 0 in the interval [ ψ-σ low ψ , ψ+σ metrized. The results are shown as a function of the true SNR p 0 /σ p,G in Fig. 17 and of the measured SNR p/σ p in Fig. 18. We observe that the MAP estimator provides slightly conservative probabilities over the whole range of SNR. The MB estimator gives low probabilities to recover the true polarization angle ψ 0 for a true SNR <1, and a measured SNR<2.

Three-dimensional case

In all of the preceding sections, the total intensity I was assumed to be perfectly known, I = I 0 . in some cases, however, this assumption is not valid as discussed by PMA1. For instance, one needs to subtract from the observed intensity signal any unpolarized component, leading to three main issues: i) the derived polarization fraction may be grossly underestimated if this is not done properly, ii) this subtraction may be subject to a relatively large uncertainty, larger than the noise on the total intensity, and could lead to diverging estimates of the polarization fraction when intensity crosses null values ; iii) this uncertainty on this unpolarized component intensity level should be included in the 3D noise covariance matrix, and propagated to the uncertainty estimates of the polarization fraction. This happens for instance when dealing with the polarization fraction of the Galactic dust component at high latitude, where the total intensity of the signal is strongly contaminated by the unpolarized signal of the Cosmic Infrared Background (CIB).

The Bayesian approach has the definite advantage over other estimators discussed here in that it can deal fairly easily with three-dimensional (I, Q, U) noise. However, an uncertain total intensity still poses problems, which are most acute in low brightness regions, since the noisy I may become null or negative, leading to infinite or negative polarization fractions. With this in mind, it is possible that the choice of the prior in p 0 and I 0 may have a strong impact on the pMB estimate. One may for instance choose to allow for negative I 0 in low-brightness regions, which implies extending the definition range of the polarization fraction to the negative part, leading to a prior defined on [-1,1]. Another possibility in this case, and possible development of the present paper, is to extend the dimensionality of the problem to include the unpolarized intensity component I offset , e.g., with a flat prior between I offset,min and I offset,max , and still imposing I 0 > 0.

Let us stress that the Bayesian approach is also currently the only one that can deal with correlation between total intensity I to Stokes Q and U. We note, however, (i) new and forthcoming polarization data sets have a much better control of these systematics, and (ii) the impact of these correlations between noise components on the polarization fraction and angle bias is quite limited, as shown by PMA1.

Conclusion

We have presented in this work an extensive comparison of the performance of the polarization fraction and angle estimators. While Simmons & Stewart (1985) focused on the common estimators of the polarization fraction, such as the Maximum Likelihood (ML), the Most Probable (MP) and the Asymptotic (AS), and Quinn (2012) suggested to use a Bayesian approach to estimate the polarization fraction, we have generalized all these methods to take into account the full covariance matrix of the Stokes parameters. We have also included in this comparison a novel estimator of the polarization fraction, the Modified Asymptotic (MAS, Plaszczynski et al. 2014). In addition, we have performed for the first time a comparison of the performance of the polarization angle estimators, since a statistical bias of ψ is expected when the covariance matrix departs from its canonical form. We have followed a frequentist methodology to investigate the properties of the polarization estimators, even when dealing with the frequentist estimators inspired by the Bayesian approach.

The question of the performance of a p or ψ estimator depends intrinsically on the analysis we would like to carry out with these quantities. Including or not the full covariance matrix is one of the first questions that must be handled, but the more important aspect relies on the properties of the output distribution of each estimator. In practice, a compromise between three frequentist criteria has to be found: a minimum bias, a minimum risk function, and the shape of the output distribution, in terms of non-Gaussianity. We present below a few recipes associated to typical use cases.

-Build a mask. It is usually recommended to build a mask on the intensity map, instead of using the SNR of the polarization fraction, so that no values of the polarization fraction (especially low values of p) are discarded in the further analysis. It can be useful, however, to build a mask based on the SNR of a polarization fraction map when we are interested in strong values of the polarization fraction only, and we try to reject p estimates artificially boosted by the noise. This is the case when we look for the maximum value of p, for example. In this context we suggest following the prescription of P14, using a combination of the MAS estimator with confidence intervals. This method allows building conservative domains where the SNR is ensured to be greater than a given threshold. P14 provide numerical approximations in the canonical case. If one wants to take into account the specificity of the noise properties in each pixel, confidence intervals can be built for any covariance matrix (including ellipticity and correlation), but it could require intensive computing. Another alternative in that case is to build credible intervals using the posterior distribution (MAP).

-Large maps of the polarization fraction with high SNR on the intensity. Another typical use case is to provide large maps of the polarization fraction with the associated uncertainty, when the intensity is assumed to be perfectly known. Because of their discontinuous distributions presenting a peak at p=0 and their strong dependance to the unknown true polarization angle ψ 0 , the common estimators of p ML, MP and AS are not well designed for this purpose. These estimators could produce highly discontinuous patterns with zero values over the output p map when the SNR goes below 4, which may imply complicated analysis including upper limits values. In order to avoid such issues, we first suggest using the MAS estimator which has been shown to produce the lowest relative bias, with a continuous output distribution which becomes close to a Gaussian for SNR larger than 2. Moreover, the relative risk function associated with the MAS estimator becomes competitive for SNR>3, while the MB estimator minimizes the relative risk function for an intermediate SNR, between 1 and 3. The uncertainties can then be derived again from the confidence or credible intervals, depending on the ellipticity of the covariance matrix. A second option, especially suited for intermediate SNR (2-3), consists in performing a preliminary analysis on the data to build a prior from the p distribution, which can then be injected into the MB estimator. The performance of this method strongly rely of the properties of the initial true distribution. It is particularly efficient for true po-larization fractions largely greater than zero, to avoid the major drawback of the MB estimator presenting a lower-limit proportional to the noise level. Hence the MB (with flat prior) estimator presents a cut-off at 0.8σ p , so that it can never provide null estimates of p. We stress that above a SNR of 4, all methods (except MP2) fall in agreement.

-Combined polarization fraction and angle analysis. The Bayesian estimators of pMB and ψMB may be used to build estimates of the polarization fraction and angle simultaneously, by taking into account the full covariance matrix, including the ellipticity and correlation between Q and U, and the correlation between total and polarized intensity. This could be useful when performing an analysis over large areas with inhomogeneous noise properties, when the SNR on the intensity becomes problematic, or when an important correlation between I and (Q, U) exists. Nevertheless we stress that the output distributions of the MB estimates are strongly asymmetric at low SNR (<3), and that the Bayesian uncertainty estimates can not be used as typical Gaussian 68% tolerance intervals.

-Low SNR on the intensity. We recommend in this case to use the Bayesian estimators which allow simultaneous estimates of the intensity and the polarization parameters, taking into account the full covariance matrix, and include the impact of the uncertainty of the intensity on the polarization fraction estimate.

-Very low SNR studies. Very low SNRs studies may require different approaches. We have seen that at low SNR, all estimators provide biased estimates of the polarization fraction, with highly asymmetrical distributions. The more conservative option in this case is to use the confidence or credible intervals. Similarly the question of assessing the unpolarized level of a set of data (i.e. SNR∼0) has been first raised by [START_REF] Clarke | [END_REF]. They suggested to used Kolmogorov test to compare the measurement distributions with the expectation derived from the Rice distribution with p 0 =0. Another option is to build the likelihood in two dimensions (Q,U) to perform a χ 2 test with Q 0 =U 0 =0. A last method is to use the Bayesian posterior probability B(p 0 |p, σ p ) to assess the probability to have p 0 =0 for a given measurement or a series of measurements by convolving all individual pdfs.

-Polarization angle. Concerning the polarization angle estimates ψ, we have shown that the ML provides the best performance in terms of relative bias and risk function for SNR>1. It corrects a potential bias of ψ when the covariance matrix is not under its canonical form. Because the ML and MAP estimators give equivalent results, the MAP can be used to efficiently build credible intervals and symmetric uncertainties, which have been shown to be in a very good agreement with the output distributions. Nevertheless we stress that the level of the absolute bias of ψ remains extremely limited compared to the dispersion of the polarization angle in most cases (i.e. in the low and tiny regime of the covariance matrix), so that it can be usually neglected.

I 0       p 2 4σ 2 p       1 F 1       3 2 , 1, p 2 2σ 2 p       -1 .
(C.10) As shown in Fig. C.1, this analytical approximation gives less than 0.15% of relative error at low SNR compared to the exact pMB estimate and less than 0.05% for the associated uncertainty. This small departure quickly tends to 0 for a SNR>4. Thus these expressions may be used to speed up the computing time when the canonical simplification may be assumed. We explore in Fig. C.2 to extent to which the canonical simplification may be done in the presence of an effective ellipticity of the covariance matrix. In this more general case, we suggest changing σ p into σ p,G in the Eqs. C.9 and C.10. The relative error between the approximate estimate and the exact bayesian estimate has been explored in two regimes of the covariance matrix, the low (1<ε eff <1.1) and tiny (1<ε eff <1.01) regimes. Three domains are observed in the top panel of Fig. C.2 dealing with the accuracy of the pMB estimate: i) at low SNR (<1) the bias on p is so large that the presence of an effective ellipticity does not affect significantly the estimate in comparison; ii) for an intermediate range of the SNR (1<SNR<4), the effective ellipticity of the Σ p significantly affects the Bayesian estimate so that the departure of the analytical approximation from the exact estimate becomes important; iii) at high SNR(>4) the noise is so low that the Bayesian estimate is not sensitive to the asymmetry of the covariance matrix anymore. Consequently, the approximate analytical expression provides very good estimates of pMB for SNR<1 and SNR>4, and 5% to 0.5% of relative error for intermediate 1<SNR<4 in the low and tiny regimes of the covariance matrix, respectively. Notice that in the extreme regime of the covariance matrix the relative error increases up to 20%.

Concerning the accuracy of the Bayesian approximate estimate σp,MB of the polarization fraction uncertainty (bottom panel), the agreement is better than 0.1% for SNR<1, and about 8% SNR>1 in the low regime, and 1% in the tiny regime.

Fig. 1 .

 1 Fig.1. Distributions of p estimates obtained with the standard estimators: naïve (black), ML (blue), MP (light green), MP2 (green) and AS (red). We assume the covariance matrix to be canonical, and a SNR of p 0 /σ p =1.

Fig. 2 .

 2 Fig. 2. Statistical fraction of null estimates of p provided by the ML, MP, MP2 and AS estimators applied on Monte-Carlo measurements, as a function of the SNR, in the canonical case.

Fig. 3 .

 3 Fig.3. Comparison of the average relative bias (top), risk function (middle) and Jarque-Bera test (bottom) of the pure measurements (naïve, black), ML (dashed blue), MP (dashed light green), MP2 (dashed green), AS (dashed red), MAS (orange) and MB (pink) p estimators in the canonical case, as a function of the the SNR p 0 /σ p . The dashed lines stand for the discontinuous estimators presenting a peak of their output distribution at p=0.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Impact of the flat prior interval upper limit on the relative Bias p performance of the MB estimator.

Fig. 7 .Fig. 8 .

 78 Fig. 7.Output distributions of the p estimates starting from a distribution of independent true values (p 0,i centered around 10% of polarization fraction (grey shaded region) shown at three levels of noise characterized by the mean SNR p 0,i /σ p,G =1, 2 and 3 (top, middle and bottom, respectively). The naïve (black) and MAS (orange) output distributions are compared to the MB output distributions obtained with three different priors: flat prior between 0 and 1 (solid pink), set to the naïve output distribution (dotted pink) and set to the true input distribution (dashed pink).

  panels) p estimators. The domain of the naïve measurements in each regime is repeated in grey shaded regions, while we show the result in orange shaded regions for the MAS and pink shaded regions for the MB estimators. It appears that the relative Bias p

Fig. 10 .Fig. 11 .Fig. 12 .

 101112 Fig.10. Impact of the effective ellipticity of the covariance matrix on the statistical relative Bias p (left column) and Risk p (right column) quantities in the extreme (light shaded region) and low (dark shaded regions) regimes, for both MAS (orange, top) and MB (pink, bottom) p estimators. The domain of the naïve measurements is repeated in grey shaded regions on both plots. The canonical case of the MAS (and MB) is also repeated on each panel in dashed orange (and pink) lines.

Fig. 13 .

 13 Fig. 13. Average Measured SNR computed over 10 000 Monte-Carlo simualtions as a function of the true SNR for four methods : Naive p/σ p,C (dark), MAP confidence intervals pML /σ p,MAP (blue), MB pMB /σ p,MB (pink) and MAS variance pMAS /σ p,MAS (orange). The covariance matrix is taken in its low regime.

Fig. 14 .

 14 Fig. 14.Comparison of the relative Bias ψ (left) and Risk ψ (right) quantities of the four ψ estimators: Naive (black), ML (blue), MP2 (green) and MB (pink) plotted as a function of the true polarization angle ψ 0 and computed at four SNR p 0 /σ p,G =0.5, 1, 2 and 5. The covariance matrix is set to ε = 2 and ρ = 0 (ε eff = 2).

Fig. 15 .Fig. 16 .

 1516 Fig.15. Statistical relative Bias ψ (left panel) and Risk ψ (right panel) averaged over ψ 0 between -π/2 and π/2, as a function of the SNR on p 0 /σ p,G , for the four ψ estimators: naïve (black), ML / MAP (blue), MP2 (green) and MB (pink). We consider two setups of the covariance matrix here: ε eff =2 (solid line) and and ε eff =1.1 (dotted line).

Fig. 17 .

 17 Fig. 17. Probability to find the true polarization angle ψ 0 inside the interval [ ψσ low ψ , ψ + σ up ψ ], where σ low ψ and σ up ψ are the lower and upper uncertainties for each estimator, ML/MAP (blue) and MB (pink), and plotted as a function of the SNR p 0 /σ p,G . Monte-Carlo simulations have been carried out in the low regime of the covariance matrix. The expected level at 68% is shown as a dashed line.

Fig. 18 .

 18 Fig. 18. Same as Fig. 17, but plotted as a function of the measured SNR p/σ p.

  Fig. C.1. Accuracy of the approximate analytical expression of the Bayesian estimates of the polarization fraction pMB (solid line) and its associated uncertainty σp,MB (dashed line), as a function of the SNR of the measurement p/σ p , where σ p =σ Q /I 0 =σ U /I 0 .

  Fig. C.2. Accuracy of the generalized approximate analytical expression of the Bayesian estimates pMB (top) and σp,MB (bottom), taking into account the full covariance matrix components, in the low (light grey ) and tiny (dark grey) regimes.

Table 1 .

 1 List of the acronyms of the estimators used in this work. The parameters to which each estimator applies, independently (/) or simultaneously (&), are given in the last column.

	Acronym Description	Parameters
	ML	Maximum Likelihood	p / ψ
	MP	Most Probable in 1D	p / ψ
	MP2	Most Probable in 2D	p & ψ
	AS	Asymptotic	p
	MAS	Modified Asymptotic	p
	MAP	Maximum A Posteriori	p / ψ
	MAP2	Maximum A Posteriori in 2D	p & ψ
	MB	Mean posterior Bayesian	I & p & ψ

  for the MAS (top panels) and MB (bottom

		0.30 0.35	MP ML Naive	AS MB MAS	p 0 true	θ=0 ε eff =2 p 0 /σ p,G =1	0.30 0.35
		0.25	MP2				0.25
		0.20					0.20
	p	0.15					0.15
		0.10					0.10
		0.05					0.05
		0.00					0.00
		-π/2 -3π/8 -π/4 -π/8	0.00 ψ 0 [rad]	π/8 π/4 3π/8 π/2	-π/2 -3π/8 -π/4 -π/8	0.00	π/8 π/4 3π/8 π/2

Planck (http://www.esa.int/Planck) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states (in particular the lead countries France and Italy), with contributions from NASA (USA) and telescope reflectors provided by a collaboration between ESA and a scientific consortium led and funded by Denmark.

Acknowledgements. This paper was developed to support the analysis of data from the Planck satellite. The development of Planck has been supported by: ESA; CNES and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN, JA, and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); and PRACE (EU). A description of the Planck Collaboration and a list of its members, including the technical or scientific activities in which they have been involved, can be found at http://www.sciops.esa.int/index.php?project=planck&page= Planck_Collaboration. We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA), part of the High Energy Astrophysics Science Archive Center (HEASARC). HEASARC/LAMBDA is a service of the Astrophysics Science Division at the NASA Goddard Space Flight Center. Some of the results in this paper have been derived using the HEALPix package. We would also like to thank P. Leahy, S. Prunet and M. Seiffert for their very useful comments.

Appendix A: Most Probable in general case

The MP2 estimators, pMP2 and ψMP2 , have to satisfy the Eqs. 8 and 9 simultaneously. These relations can be solved, using the fully developed expression of f 2D including the terms of the inverse matrix Σ -1 p :

with

This analytical solution only depends on the input measurements (p, ψ) and the covariance matrix Σ p . Because the polarization fraction must be positive, there exists a lower limit of the SNR so that pMP2 = 0. In that case ψMP2 is not constrained any more and can be chosen to any possible value, we will set it equal to the measurement ψ. Moreover, this expression can be simplified when ρ = 0, which implies that v 12 = 0, leading to:

.

In the canonical case (v 12 =0, v 11 =v 22 =1/σ 2 p ), we recover the expression derived by Quinn (2012)

Appendix B: Bayesian Posterior pdf

We illustrate the shape of the posterior pdf in Fig. B.1, where B 2D (p 0 , ψ 0 | p, ψ, Σ p ) is shown at four levels of the SNR and five couples of (ε, ρ). It is interesting to notice that the posterior pdf allows the polarization fraction to be null at low SNR, when these values were rejected by the pdf (see Appendix B of PMA I). Moreover the posterior pdf peaks at the location of the measurements used to compute it. As largely emphasized in PMA I, we also recall that once the effective ellipticity of the covariance matrix departs from the canonical simplification, the pdfs are sensitive to the initial true polarization angle ψ 0 .

Appendix C: Mean Bayesian Posterior analytical expression

In the canonical case, the MB estimator of the polarization fraction p takes a simple analytical expression. The Bayesian posterior on p is given in this case by:

Because the uncertainty becomes small compared to the polarization fraction at high SNR, up to 8% of error in σp,MB is still acceptable for this approximation.