5356

Research Article

Vol. 57, No. 19 / 1 July 2018 / Applied Optics

Statistics of normalized Stokes
polarization parameters
NATHAN HAGEN
Utsunomiya University, Department of Optical Engineering, Center for Optical Research and Engineering (CORE), 7-1-2 Yoto, Utsunomiya,Tochigi 321-0954, Japan (nh@hagenlab.org)
Received 27 April 2018; revised 31 May 2018; accepted 31 May 2018; posted 4 June 2018 (Doc. ID 330398); published 27 June 2018

The normalized Stokes parameters are formed from the ratio of the polarization components to the intensity
component of light. Such ratio distributions are known to have an undefined mean and variance, and yet
researchers in the polarization community work with these normalized parameters all the time. We show thatwhile in theory the normalized parameters have a pathological probability density, in practice they are quite wellbehaved. We provide expressions for their approximate densities and confirm the results with laboratory
measurements. © 2018 Optical Society of America
OCIS codes: (260.5430) Polarization; (280.0280) Remote sensing and sensors; (230.5440) Polarization-selective devices;(030.4280) Noise in imaging systems.
https://doi.org/10.1364/AO.57.005356

1. INTRODUCTION
Stokes polarimeters aim to estimate the polarization state oflight incident on the instrument’s entrance pupil via the four
Stokes parameters s  s0 , s 1 , s2 , s 3 T . Rather than simply using
them in this general form, users often convert the latter threeparameters into a normalized form s˜i  si ∕s 0 to make themeasier to interpret: linear polarized light at 0°, for example, givess˜  1, 0, 0, so that it is easy to see the general polarization
state at a glance. To distinguish between s0 and the remaining
three Stokes parameters s 1 , s2 , s 3 , we will use only the index ito refer to the latter.
The statistical properties of Stokes polarimeter measurements has been the subject of a number of research studies,in which the probability density for the four Stokes parameter
estimators is typically approximated with a four-dimensional(4D) multivariate Gaussian probability density [1,2]. From this
4D density we can predict the precision and accuracy for agiven measurement setup. While the statistics of the raw
Stokes parameters have received a good deal of attention,
the statistics of the normalized parameters have not. One likely
reason that analysis of the normalized Stokes parameters has
been neglected is that they are represented by the Gaussian ratio
distribution—a pathological probability density that is wellknown to have no first or higher moments. For such a distribution, one cannot define a mean or variance, such that data of
this sort may seem of little utility. Nevertheless, we use them all
the time.
In the discussion below we will see that, by using our priorknowledge of polarization measurements and of the physicallimits of the normalized Stokes parameters, we will be able
1559-128X/18/195356-08 Journal © 2018 Optical Society of America

to define and calculate appropriate mean and variance values.
Indeed, for almost all practical situations for polarimetric measurements, the densities of the normalized parameters can beclosely approximated with a Gaussian density. Thus, we
provide a theoretical justification for use of the normalized
Stokes parameters and show how they can be characterized.
2. GAUSSIAN RATIO DISTRIBUTION
Since each normalized Stokes parameter involves taking theratio si ∕s0 , the statistical properties of s˜i can be derived fromtaking the ratio of one Gaussian distribution with another.
To facilitate the discussion, we will use the variable y in placeof the denominator s 0, and the variable x in place of thenumerator—any one of the other Stokes parameters si .
Thus, we consider the probability density prz produced fromthe quotient z  x∕y of two Gaussian-distributed variables x
and y.
The density prz can be expressed in terms of the means μand standard deviations σ of the variables x and y, and theircorrelation coefficient ρ, as [3]pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1 − ρ2−c
exp pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
prz πσ x σ y a2 z
2 1 − ρ21bzd zgz
erf pﬃﬃﬃ ,(1)
 pﬃﬃﬃﬃﬃ
2π σ x σ y a3 z
2
where


 2
z
2ρz
1 1∕2
−

,
az 
σ 2x σ x σ y σ 2y(2)Research Articlebz 

Vol. 57, No. 19 / 1 July 2018 / Applied Optics

μx z ρμx  μy z μy
−
 2,
σy
σx σy
σ 2x

(3)

μ2x 2ρμx μy μ2y
−
 2,
σ 2x
σx σy
σy

(4)

c


 2
b z − ca2 z,d z  exp
21 − ρ2 a2 zgz 

bz
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ,az 1 − ρ2

(5)

(6)

and “erf ” is the standard error function. Figures 1–4 show examples of prz for four different sets of input variables. Fromthe examples, we can see that the density can be bimodal,
highly asymmetric, or nearly Gaussian in shape dependingon the form of prx and pry.
Figure 1 shows prz for uncorrelated x and y when μy  0.
The distance between the two peaks increases as jμx − μy j increases, and when μx  μy  0, prz becomes single modeand symmetric. The gray bars in the figure indicate a histogramof 106 samples obtained by computer simulation from generating

Fig. 1. Ratio distribution for parameter set μx , μy , σ x , σ y , ρ 
−2, 0, 1, 1, 0. The density becomes symmetric when μy  0 andρ  0. The gray rectangles show the histogram data obtained fromsimulating the ratio of Gaussian-distributed x and y variables. The blue
curve gives the exact ratio density prz from Eq. (1), the red curve theapproximate form [Eq. (8)], and the green curve the Gaussian approximation [Eq. (11)]. The Gaussian approximation curve here has zeroamplitude and infinite variance.

Fig. 2. Ratio density becomes asymmetric when μy > 0.

5357

x and y separately and taking their ratio. The blue curve in the
figure shows the exact probability density [Eq. (1)], while the red
curve shows the approximate form [Eq. (8)], and the green curve
the Gaussian approximation to the density [Eq. (11)].
Figure 2 shows prz for μy > 0, where we can see that thedensity is now asymmetric. Figure 3 shows a similar setup inwhich the data is correlated with ρ  0.25. The normal
approximation is still not useful for either of these cases, but
as μy ∕σ y increases, we find that the normal approximation fits
more and more closely to the exact distribution, as in Fig. 4.
While the normal approximation becomes a closer fit to theexact distribution as μy ∕σ y increases, it remains true that theexact distribution has an undefined mean and infinite variance.
While the complex five-parameter μx , μy , σ x , σ y , ρ dependence of prz makes it difficult to analyze its behavior,
Geary and co-workers [3–5] have each shown that this fiveparameter set can be transformed into an equivalent twoparameter form. We can define these two parameters ξ, η as
ξ

μx
,
σx

η

μy ∕σ y  − ρμx ∕σ x 
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ:
1 − ρ2(7)

If inputs x and y are uncorrelated (i.e., ρ  0), then ξ and η aresimply the reciprocal coefficients of variation: ξ  μx ∕σ x andη  μy ∕σ y . In addition, for analyzing the behavior of the ratio
distribution, we need only consider x ≥ 0, since flipping thesign of x simply flips the density about the origin.

Fig. 3. Example ratio density when ρ ≠ 0. If we use the same
parameter inputs but assume ρ  0 (dashed blue curve), we cansee that the correlation has significant effects on the density.

Fig. 4. Ratio density becomes more Gaussian-like as η  μy ∕σ ybecomes large.

5358

Research Article

Vol. 57, No. 19 / 1 July 2018 / Applied Optics

Although this transformation from five to two parameters iscompletely general, some reports in the research literature havepropagated an error that the generality does not hold [3,6,7].
Another situation in which the five-parameter set is reduced to
two is for shot-noise-limited data. In this case, the distributions
for x and y will have σ 2x  μx , σ 2y  μy , and ρ  0.
Marsaglia has also suggested that if y can be assumed to be
non-negative, then the density can be approximated by one
having a simpler analytic structure [7]:prM z  0.39894228 exp−0.5w2 ju∕v 3 j,

(8)9u  σ x μy σ x − μx ρσ y   zσ y μx σ y − σ x ρμy , =v  σ 2x − 2zσ x σ y ρ  z 2 σ 2y 1∕2 ,;w  zμy − μx ∕v:

(9)

where

This approximation is shown as a red curve in the figures.
The exact form of prz in Eq. (1) shows that it can be considered as the sum of two terms, one of which is a weighted
Cauchy density pr1 z and the other an asymmetric bimodal
Gaussian-like density pr2 z [8]. The Cauchy density pr1 z
has an amplitude A that decreases rapidly with increasing ξand η:
 2

1
ξ  η2:(10)A  exp −π2Thus,pcomponent
pr1 z becomes negligible when ρ is small
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
and ξ2  η2 > 2.5. Figure 5 shows an example of the fulle(a)density function prz and its two components pr1 and pr2 .
In linear scale, one cannot see the second mode of pr2 , but
in the logarithmic plot one can easily see the separate modes.
While the second component pr2 z is inherently bimodal,
we can describe it as “Gaussian-like” because the density
asymptotically approaches a Gaussian density as η becomeslarge (see Appendix A). While it is always bimodal, one modeof pr2 shrinks rapidly in amplitude as ξ and η increase from 0,
while the remaining mode does not shrink in amplitude andapproaches a Gaussian shape as η increases. Thus, for large η,a Gaussian approximation becomes a good one, and we canestimate the mean and variance parameters for the approximately Gaussian-distributed z by taking a third-order
Taylor expansion of z  x∕y in the neighborhood of
x, y  μx , μy , giving [9–12]μz G 

μx σ 2y μx ρσ x σ y 3 − 2 ,μyμyμyσ 2 σ 2y μ2x 2ρσ x σ y μx:
σ 2z G  2x  4 −
μy
μ3y
μy

(11)

This is a close match to the exact density whenη > 3, jξj < η∕3.
The Gaussian approximation is shown as a green curve inFigs 1–4. While Fig. 4 shows that the Gaussian approximation
closely matches the exact ratio density for large η values, theratio density’s mean and variance are nevertheless undefined,a fact that is easily confirmed in simulation. This is the resultof the Cauchy component pr1 . For large η it has a very low
amplitude and a peak that is often a large distance from themain mode, making it difficult to detect visually. However,its extended tails produce an infinite variance, despite its
low amplitude.
3. RATIO DISTRIBUTION FOR STOKESPARAMETERS(b)

Fig. 5. Example density showing the Cauchy component pr1 and thebimodal component pr2 for density prz [Eq. (1)], given inputs ξ  1
and η  1.5. While the low-amplitude second mode of pr2 is not visiblein (a) the linear plot, it is easily visible in (b) the logarithmic plot.

For Stokes parameters, s 0 is always non-negative, and the case ofs0  0 is uninteresting, so that the denominator y of the ratio
distribution can be considered strictly positive. In addition, by
the physical constraints of how polarimetric measurements canbe performed, we know that σ x is always larger than σ y . Thus,
the poor fit of the Gaussian approximation to the exact distribution (e.g., Figs 1–3) apply primarily to nonphysical situations. Finally, if we consider shot-noise-limited data where
the average signal si is 25 photoelectrons or more, thenη ≳ 5, and prz closely matches a Gaussian density in shape.
Thus, while the ratio distribution retains its undefined mean
and variance in theory, in practice any of the rare data points
lying far outside the physical region will be considered outliers
rather than valid measurements.
The setup for Fig. 4 is a realistic situation that can occur for
fully polarized light. For shot-noise-limited measurement, thissetup corresponds to a weak light level of only 100 photons
collected over the duration of the measurement. For a more
typical value of 104 (modern detector array pixels typicallycollect >104 photons before saturation), the Gaussian approximation becomes so good that the two curves are almost
indistinguishable. The Gaussian approximation [Eq. (11)] isResearch Article

Vol. 57, No. 19 / 1 July 2018 / Applied Optics

therefore a good representation for the vast majority of polarimetric measurements.
For optimized full Stokes polarimeters, we find that thepvariﬃﬃﬃ
ance in s1 , s2 , and s3 is 3 times that in s0 , giving σ x  3σ y[2,13,14].
pﬃﬃﬃ (For linear Stokes polarimeters, this becomesσ x  2σ y .) As a result, Eq. (11) becomeshs i hs i 1 pﬃﬃﬃ 1
h˜si iG  i  i 2 − 3ρ 2 ,hs0 i hs0 i ηηpﬃﬃﬃ
2
hs i 13 2 3ρ 1
σ 2s˜i G  i 2 2  2 −,(12)
hs0 i ηη
hs 0 i η2
where η  hs0 i∕σ s0 ≈ SNR, subscript “G” indicates theGaussian approximation, and angle brackets h·i indicate taking
the mean value. For shot-noise-limited data, we can also write
hs0 i ≈ SNR 2 , so that in the limit η ≫ 1 these relations become
3  h˜s i i2
:
(13)
SNR 2
For a linear Stokes polarimeter, the 3  hσ s˜i i2  here
becomes 2  hσ s˜i i2 .
Thus, in order to estimate the statistics of the normalized
Stokes parameters, we need only check that hs0 i > 25 photoelectrons to apply either the Gaussian approximation of
Eq. (11) or the simplified form [Eq. (13)].
h˜si iG ≈ h˜si i,

σ 2s˜i ≈4. PROBABILITY DENSITY FOR THE DEGREE
OF POLARIZATION
From the individual probability densities of the normalized
Stokes parameters, we can also form the density of the degreeof polarization (DOP) p, estimated as
qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ(14)p  s˜21  s˜22 ,for a linear Stokes polarimeter. Taking the Taylor expansion of
Eq. (14), the variance of p can be written in terms of the meansand covariances of the s˜i as 2    ∂p  2  ∂p 2 2
 ∂p  ∂p 
2




σ p ≈   σ s˜1    σ s˜2  2  cov˜s1 , s˜2 
∂˜s1
∂˜s2
∂˜s1 ∂˜s2
1≈ 2 hs˜1 i2 σ 2s˜1  h˜s2 i2 σ 2s˜2  2h˜s1 ih˜s 2 icov˜s1 , s˜2 :(15)pNote that the s˜1 and s˜2 estimates generally have significant correlation, so that use of the correlation term in σ 2p is important
for accurate estimates.
To find the dependence of p statistics on the original (unnormalized) Stokes vector, we substitute the definitions of σ 2s˜1and σ 2s˜2 from Eq. (12) into Eq. (15), obtainingh˜s i2h˜s i2σ 2p  A1 σ 2s0 1 2  σ 2s1  A2 σ 2s0 2 2  σ 2s2 , (16)
hs0 i
hs0 i
where
Ai 

hs i i2 ∕hs 0 i2:
hs1 i2  hs2 i2(17)

Section 6 shows an experimental distribution for p and comparison with Eq. (15).

5359

For a full Stokes polarimeter, the DOP estimator becomes
qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
p  s˜21  s˜22  s˜23 ,(18)
so that the corresponding three-variable Taylor expansion is  
 2
 2
 2
 ∂p  ∂p  ∂p  ∂p  ∂p 
σ 2p ≈ s˜21    s˜22    s˜23    2  cov˜s1 , s˜2 
∂˜s 1
∂˜s 2
∂˜s3
∂˜s 1 ∂˜s 2
     ∂p  ∂p  ∂p  ∂p 
 2  cov˜s 1 , s˜3   2  cov˜s2 , s˜3 ∂˜s1 ∂˜s3∂˜s2 ∂˜s3
1
 2 hs˜1 i2 σ 2s˜1  h˜s2 i2 σ 2s˜2  h˜s 3 i2 σ 2s˜3  2h˜s 1 ih˜s2 icov˜s 1 , s˜2 p 2h˜s1 ih˜s3 icov˜s1 , s˜3   2h˜s 2 ihs˜3 icov˜s2 , s˜3 :(19)

Goudail and Bénière analyzed the variance of the DOP
estimator directly from N measurements of a rotating linear
polarizer, without first forming the normalized Stokes parameters [15,16]. In the presence of independent identically
distributed Gaussian (IG) noise, or in the presence of pure
Poisson (PP) noise, they obtained the expressions
σ2σ 2p IG  16 20 2  μ2p μ0andσ 2p PP 22 − μ2p : (20)μ0Comparing these to Eqs. (15) and (16), we find that the
Poisson-noise variances estimated by Goudail and Bénière
are larger by a factor of 4 due to the use of N (the numberof measurements) in their definition of the measured intensities. Once this factor is removed, the expressions agree
closely.
5. RECTIFICATION OF PARAMETER ESTIMATES
Beyond simply taking the ratio of s i and s0 to get the normalized parameters, it is also a common practice among the polarization community to rectify the calculated s˜i such that they lie
within the physical range p ≤ 1, enforcing that the DOP cannot exceed 1. Any values obtained outside the physical range areforced back within range, so that, for example, if we obtain an
estimate s˜1  1.1, then this result is reset to s˜1  1. In terms of
the probability density of the resulting normalized parameters,
this creates a mixed continuous-discrete density with δ functions at the edges of the physical range [17]:prr z  A−1 δz  1  rectz∕2prz  A1 δz − 1,(21)
where δz is the Dirac delta function, rectz is the standard
rectangle function, used here to limit the range between −1 and1, andZ −1
Z ∞przdz,A1 
przdz:
(22)
A−1 
−∞

1

The “r” subscript indicates the post-rectification version of theprobability density. Figure 6 shows an example rectified density[18], in which the curve represents a probability density, and a
dot a probability value (for the δ-function).
If we consider all three normalized Stokes parameters
together, the distribution takes a more complicated form than
the 1D density of Eq. (21), since the rectification now takes placeat the surface of a sphere, outside of which p > 1. We define the(longitude, latitude) representation of the polarization state on

5360

Research Article

Vol. 57, No. 19 / 1 July 2018 / Applied Optics

6. EXPERIMENTAL RESULTS

Fig. 6. Example Gaussian probability density (μ  0.5, σ  0.6)
rectified to the range −1 ≤ s˜i ≤ 1. The curve gives a probability density (as indicated on the y axis), but the dots give probability values(and thus have arbitrary scale on this axis). The dashed curves indicate
the pre-rectified density.

the Poincaré sphere with α, ϵ. The 3D rectified density cantherefore be written as
˜  Aα, ϵδp − 1  Hpprs,˜
prr s(23)

where p is the radial coordinate in 3D polarization space, δ· isthe Dirac delta function, H · is the Heaviside step function,used here to limit the range of the probability density to be
between 0 and 1, andZ ∞prsdp˜(24)Aα, ϵ 1is the 2D amplitude weighing the δ-functions in sphericalangular coordinates:α  arctan˜s2 ∕˜s1 ,ϵ  arccos˜s 3 ∕p:(25)

The subscript “r” here marks the density as “rectified.”
This post-rectification density prr s is obtained by first substituting in the pre-rectification density prs and solving theresulting radial integral to get the δ-function amplitude at thesphere’s surface (p  1).
The nonlinear rectification operation eliminates the pathology of the ratio distribution from the probability density for thenormalized Stokes parameters, both in practical terms and intheory. With rectification, all moments of the posterior density
for the parameters become defined.
When the pre-rectification density for the DOP prp
is negligible outside the physical range (i.e., prp ≈ 0 for
p > 1), rectification has little effect on the apparent shape ofthe distribution, since the long tail that is removed is typicallyof such low amplitude that it can be neglected in practice.
When the DOP is near 1, however, the density prp can have
half of its curve outside the physical limit, and when this data is
rectified by the user, the post-rectification density is no longer
approximately Gaussian (e.g., Fig. 8).
Prior to rectification, the estimates for the s˜i are approximately unbiased, so that rectification will introduce bias.
The width of the distribution is reduced from its original,
so that the variance shrinks (i.e., improves). Finally, the prerectified data are all necessarily farther from the true value thanare the post-rectification values, so that the overall error is
reduced.

To confirm the expressions obtained in Sections 3 and 4, weconstructed a simple linear Stokes polarimeter, as shown in
Fig. 7. Polarizer P 1 is used to define the incident polarization
azimuth θ, so that it together with the light source can be considered the polarization state generator. Polarizer P 2 is rotatedto azimuth angles 0°, 45°, 90°, and 135°, with image data collected at each position to estimate the incident polarization
state. The diffuser (a plastic disk abraded on both sides) is used
to ensure that uniform unpolarized light is incident on P 1 ,while the optical bandpass filter (651.2–661.2 nm, CVI
Laser LLC p/n 50523270) restricts the measurement to a spectral region where the polarizers have high extinction. The detector array is a 6.8 mm × 5.4 mm (1280 × 1024 pixel) camera
core (Basler ace acA1280-60gm). This lensless layout is used toprevent any unwanted (and difficult to predict) polarization
effects from lenses, and the light source was placed at sufficientdistance (55 mm) to ensure that the angles of incidence on thepolarizers are <5° to minimize angle-dependent effects [19].
To model Poisson noise correctly, we calibrate the detector
output to photon number using the slope of the photon transfer curve [20].
One worry about this experimental design is that having arotating polarizer in front of the detector can introduce artifacts
due to polarization dependence in the detector itself. In such acase, it is possible to remove this issue by mounting polarizer P 2
and the camera core together and rotating the two in tandem.
This step proved to be unnecessary here.
If we write the measured intensity I at each of the four orientation angles of P 2 as I 0 , I 45 , I 90 , and I 135 , then the linear
Stokes parameters can be estimated from91
I 0  I 45  I 90  I 135  >ŝ 0  2ϵ=ŝ 1  1ϵ I 0 − I 90 
:
(26)
>
;
1
ŝ 2  ϵ I 45 − I 135 
By collecting a set of 5000 images and analyzing the noise
behavior over the time axis, and repeating the collection process
for a progressive series of light intensities, we can verify that thenoise is Poisson dominated, with a relatively small Gaussiancomponent derived from detector noise. Figure 8 shows thehistograms of the Stokes parameter estimates calculated for
an incident linear polarization state of θ  22.5°, demonstrating that the resulting si are also accurately approximated as
Gaussian.

Fig. 7. Experimental layout for measuring the Stokes parameters forlinearly polarized light oriented at θ, where the incident polarization
state is defined by P 1 θ.

Research Article

Vol. 57, No. 19 / 1 July 2018 / Applied Optics

5361

Table 1 shows numerical measurements of parameter standard
deviations for incident polarization azimuth angles of 0° and
22.5°, showing that measurements agree with the predicted values to within 3%. It is important to note that, when making thes˜i variance estimates from Eq. (11), the third term in the equation—incorporating the effect that correlation among the input
Stokes parameters has on the s˜i variances—can be significant
and should not be neglected. While the 0-45-90-145 measurement configuration of Eq. (26) obtains Stokes parameter estimates with negligible covariances in the case of independent
identically distributed Gaussian noise, the covariances are significant for the more common cases of pure Poisson noise andmixed Poisson–Gaussian noise, as shown in Table 1 [2].
7. IMPACT OF PARAMETER VARIANCES ON
SYSTEM OPTIMIZATION

Fig. 8. Experimental data showing the normalized histograms of
measured Stokes parameters s0 , s1 , and s 2 ; the normalized Stokes
parameters s˜1 and s˜2 ; and the DOP p, for an input state of linear polarization at 22.5°. Each histogram is shown together with its corresponding Gaussian fit [solid curve, Eq. (11)]. The rectified DOP distributionis also shown at the bottom right (note that the vertical axis is in logarithmic scale).

From the input Stokes parameters, we calculate the normalized Stokes parameters s˜i and DOP p (Fig. 8). The corresponding histograms indicate that the normalized parameter and
DOP distributions are closely approximated as Gaussian.
Table 1. Measured and Predicted Means and Standard
Deviations of the Measurement Data Shown in Fig. 8 for
Input Linear Polarization Orientations 0° and 22.5°a
Meas.
Mean
s0 :
s1 :
s2 :
s˜1 :
s˜2 :
p:

6713
6553
−106
0.9762−0.0157
0.9764

s0 :
s1 :
s2 :
s˜1 :
s˜2 :
p:

6583
4753
4474
0.70090.66490.9662

a

Pred.
Mean

Meas.
Stdev

Pol. State: 0°
—
61.32
—
85.20
—
85.83
0.97620.0091−0.0157
0.01280.97630.0090Pol. State: 22.5°
—
61.32
—
87.22
—
83.34
0.70090.01130.66490.01080.96610.0093

Pred.
Stdev

ρs i , s 0 ———
0.00930.01280.00931.000
0.679−0.044——————0.01130.01100.00941.000
0.508−0.007———The inputs s 0 , s 1 , and s 2 are in units of photoelectrons.

One question that we can address with the above results deals
with system optimization. If the aim of a polarimetry systemdesigner is to measure the normalized Stokes parameters s˜accurately, does the ideal system design differ from one thataims to measure the unnormalized parameters s? That is, since
s0 occurs in the denominator of all the normalized parameters,
perhaps it will improve the design if we weight the s 0 noise
more than the si noise rather than choosing an equally weighted
variance (EWV).
Reducing the variance on s0 at the expense of the other
Stokes parameters may be possible using imperfect diattenuating elements instead of ideal polarizer elements—improving
system transmission at the cost of reducing the modulation
efficiency. Equation (16) gives the variance of p in terms of
the input Stokes parameters. From this expression we can seethat the two scale factors A1 and A2 depend on the state ofpolarization, but not on the estimator variances, so they do
not play a role in determining how much we should weightσ s0 versus the σ si . Only the parenthetical terms in Eq. (16) tell
us how to weight them. These terms indicate that, when the
DOP of a measurement is small, the s0 variance terms in
Eq. (16) will be dominated by the respective s1 and s 2 variances,which are generally several times larger than the s0 variance inany case. That is,σ 2p ≈ A1 σ 2s1  A2 σ 2s2 :(27)Thus, if we wish to reduce the variance of the DOP measurement, then the best option would be to optimize the instrument for reducing the σ si at the expense of σ s0 rather thanthe other way around. This essentially supports the existing
methodology used in the polarization community, and justifies
choosing ideal polarizer elements in instruments rather than
non-ideal diattenuators.
8. CONCLUSION
While in theory the posterior probability distributions of thenormalized Stokes parameters have pathological properties,
in that the mean and variance are undefined, we have foundthat in practice the distributions are generally well-behaved,
producing probability distributions that are closely approximated as Gaussian. It is only when the light level isextremely low—a mean number of photons per measurement

5362

Research Article

Vol. 57, No. 19 / 1 July 2018 / Applied Optics

less than 10—that the posterior distributions differ greatly fromtheir approximately Gaussian shape.
While the results of Section 3 indicate that the pathology of
the ratio distribution is usually not a problem for polarimetry,
one place where the ratio distribution can cause trouble is whencalculating the orientation angles of physical parameters. In
Stokes polarimetry this occurs when calculating the orientationof linear polarization, given by θ  12 arctan˜s 2 ∕˜s 1  [21]. In thiscase the denominator can take on values close to zero, so thatthe problems of an undefined mean and unbounded variancebecome a real concern. The arctangent function mitigates mostof this trouble by smoothly mapping large positive and negativevalues from the ratio to angles near π and −π. However, when
the numerator s2 is also close to zero, then the arctangent function will have a linear mapping behavior, so that the resulting
orientation angle will have the bimodal distribution shown in
Fig. 1. That is, when the degree of linear polarization is small(s1 and s2 both close to zero), θ will tend to oscillate between
the two modes. A similar situation occurs in Mueller polarimetry, when deriving orientation angles from normalized Mueller
matrix elements. In this case, the linear polarizance azimuth
angle is given from 12 arctanm̃20 ∕m̃10 , linear diattenuation azimuth angle from 12 arctanm̃02 ∕m̃01 , and the linear retardance
azimuth angle from 12 arctanm̃23 − m̃32 ∕m̃22  m̃33 .

APPENDIX B: CODE FOR CALCULATIONS
Code written in Scientific Python is given below for generating
the exact and approximate forms of the ratio distribution probability densities. The variables mx, my, sx, and sy represent
μx , μy , σ x , and σ y respectively.
The exact ratio density [Eq. (1)] (blue curve in Figs. 1–4) [3]:Marsaglia’s approximate ratio density [Eq. (8)] (red curve in
Figs. 1–4) [7]:The Gaussian approximation to the ratio density [Eq. (11)](green curve in Figs 1–4) [10]:APPENDIX A: ASYMPTOTIC GAUSSIANBEHAVIOR OF THE BIMODAL COMPONENTMarsaglia shows that the asymmetric bimodal density pr2 z
can be written in the form [5]
Z q
q
1 2
2
e −2u q  du, (A1)pr2 z 2 η2 ∕22
ξ− 1 0π1  z epﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
where q  ξz  η∕ 1  z 2 . For parameter inputs wherez ≪ 1 and jξj < η, we can approximate q as η1 − 12 z 2 , so thatZ
η
3 2 −q2 ∕2 q −u2 ∕2
pr2 z ≈edu (A2)
1− z e22πe η ∕2 − 1
0


 Z
2
η
ηe −η ∕2
1 2
1 2 2
2
exp − η 1 − z
e −u ∕2 du: (A3)≈ η2 ∕2
2
2
πe
− 1
0
Next, whentheintegral in the expression
ﬃﬃﬃﬃﬃﬃﬃﬃ
R we assume η ≫ 1,p
becomes 0∞ exp−u2 ∕2du  π∕2, so that pr2 further simplifies to2ηe η ∕2
1 2
1 4
2
(A4)
pr2 z ≈ pﬃﬃﬃﬃﬃ exp − η 1 − z  z 
2
4
2π


η
1 2 2
(A5)
ηz :
≈ pﬃﬃﬃﬃﬃ exp
22πThis has the form of a Gaussian density with a width parameter
proportional to η−1 .
This derivation assumes z ≪ 1, which is valid whenμx ≪ μy , i.e., when the measured degree of polarization islow. However, we have been able to observe empirically, andthrough simulation, that this Gaussian approximation stillholds even for DOP near 1 when η ≫ 1 and jξj < η.

Note that this PDF is undefined for μy  0.

REFERENCES AND NOTES
1. F. Goudail, “Noise minimization and equalization for Stokes polarimeters in the presence of signal-dependent Poisson shot noise,” Opt.
Lett. 34, 647–649 (2009).
2. N. Hagen and Y. Otani, “Stokes polarimeter performance: generalnoise model and analysis,” Appl. Opt. 57, 4283–4296 (2018).
3. D. V. Hinkley, “On the ratio of two correlated normal random variables,” Biometrika 56, 635–639 (1969).
4. R. C. Geary, “The frequency distribution of the quotient of two normal
variates,” J. R. Stat. Soc. 93, 442–446 (1930).
5. G. Marsaglia, “Ratios of normal variables and ratios of sums of uniform
variables,” J. Am. Stat. Assoc. 60, 193–204 (1965).
6. D. V. Hinkley, “Correction: ‘On the ratio of two correlated normal
random variables’,” Biometrika 57, 683 (1970).
7. G. Marsaglia, “Ratios of normal variables,” J. Stat. Softw. 16, 4–13
(2006).
8. J. H. Curtiss, “On the distribution of the quotient of two chance variables,” Ann. Math. Stat. 12, 409–421 (1941).
9. H. H. Ku, “Notes on the use of propagation of error formulas,” J. Res.
Natl. Bur. Stand. 70, 263–273 (1966).
10. J. Hayya, D. Armstrong, and N. Gressis, “A note on the ratio of
two normally distributed variables,” Manage. Sci. 21, 1338–1341(1975).
11. M. P. Silverman, W. Strange, and T. C. Lipscombe, “The distribution ofcomposite measurements: how to be certain of the uncertainties inwhat we measure,” Am. J. Phys. 72, 1068–1081 (2004).
12. D. T. Holmes and K. A. Buhr, “Error propagation in calculated ratios,”
Clin. Biochem. 40, 728–734 (2007).
13. F. Goudail, “Optimization of the contrast in active Stokes images,”
Opt. Lett. 34, 121–123 (2009).

Research Article
14. R. Perkins and V. Gruev, “Signal-to-noise analysis of Stokes
parameters in division of focal plane polarimeters,” Opt. Express 18,
25815–25824 (2010).
15. A. Bénière, F. Goudail, M. Alouini, and D. Dolfi, “Precision of degree ofpolarization estimation in the presence of additive Gaussian detector
noise,” Opt. Commun. 278, 264–269 (2007).
16. F. Goudail and A. Bénière, “Estimation precision of the degreeof linear polarization and of the angle of polarization in thepresence of different sources of noise,” Appl. Opt. 49, 683–694(2010).
17. Note that rectified distributions differ from the better-known truncateddistributions, in that, for Eq. (21), the latter define A−1  A1  0 andrescale the remainder of the density to achieve unit area. Thus,

Vol. 57, No. 19 / 1 July 2018 / Applied Optics

5363

the latter are continuous rather than mixed continuous–discrete
distributions.
18. N. D. Socci, D. D. Lee, and H. S. Seung, “The rectified Gaussian distribution,” Adv. Neural Inf. Process. Syst. 10, 350–356 (1998).
19. X. Lu, W. Jin, L. Li, X. Wang, S. Qiu, and J. Liu, “Theory and analysis
of a large field polarization imaging system with obliquely incident
light,” Opt. Express 26, 2495–2508 (2018).
20. S. B. Howell, Handbook of CCD Astronomy, 2nd ed. (Cambridge
University, 2006).
21. Note that orientation angle θ differs by a factor of 2 from α defined in
Eq. (25) because θ represents an angle in physical space, while α
represents the longitudinal angle in the polarization space of the
Poincaré sphere.