- Email: [email protected]
Parametric modeling of photometric signals
Signal Processing 82 (2002) 649 – 661
www.elsevier.com/locate/sigpro
Parametric modeling of photometric signals A. Ferraria; ∗ , J.Y. Tourneretb , G. Alengrina a UMR
6525 Astrophysique, Universite de Nice Sophia-Antipolis, Parc Valrose, 06108 Nice cedex 2, France b ENSEEIHT=TESA, 2 rue Camichel, BP 7122, 31071 Toulouse cedex 7, France
Received 25 February 2001; received in revised form 2 August 2001; accepted 12 September 2001
Abstract This paper studies a new model for photometric signals under high .ux assumption. Photometric signals are modeled by Gaussian autoregressive processes having the same mean and variance denoted Constraint Gaussian Autoregressive Processes (CGARP’s). The estimation of the CGARP parameters is discussed. The Cram2er Rao lower bounds for these parameters are studied and compared to the estimator mean square errors. The CGARP is intended to model the signal received by a satellite designed for extrasolar planets detection. A transit of a planet in front of a star results in an abrupt change in the mean and variance of the CGARP. The Neyman–Pearson detector for this changepoint detection problem is derived when the abrupt change parameters are known. Closed form expressions for the Receiver Operating Characteristics (ROC) are provided. The Neyman–Pearson detector combined with the maximum likelihood estimator for CGARP parameters allows to study the generalized likelihood ratio detector. ROC curves are then determined using computer simulations. ? 2002 Elsevier Science B.V. All rights reserved. Keywords: Abrupt change detection; Photometric signal; Extra solar planet detection
1. Introduction The detection of extrasolar planets is a challenging problem in astronomy (see the extrasolar planets encyclopedia [18]). Among the methods currently pursued to detect extrasolar planets, the transit method (also referred as photometric method or occultation method) may be the only one to >nd earth class planets in the near future. The transit method is based on the detection of photometric .ux variations which results from the transit of a planet in front of a star. Several space-based projects have been proposed to achieve this goal including the US Kepler project, the Eddington european project and the COROT (COnvection ∗ Corresponding author. Tel.: +33-9207-6349; fax: +33-92076321. E-mail address: [email protected] (A. Ferrari).
and ROTation) french project. For a star like the Sun, the typical relative .ux variation for a Jupiter-like giant planet is 10−2 (during 25 h) whereas for a Earth-like planet it is 10−4 (during 13 h). Consequently, the photometric .ux variations caused by terrestrial planets are diGcult to detect using the conventional transit method and new detectors have to be investigated. The >rst problem addressed in this paper is the modeling of the photometric signals received by the satellite. Based on two realistic assumptions, the photometric signals are modeled by Gaussian autoregressive processes having the same mean and variance, denoted Constraint Gaussian Autoregressive Processes (CGARP’s). The second problem addressed in this paper is the detection of extrasolar planets by using standard estimation and detection tools. The theoretical distribution of photometric signals being
0165-1684/02/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 1 ) 0 0 2 1 2 - 2
650
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
generally intractable, most extrasolar planet detectors reduce to measure the brightness drop of a star which results from the transit of one of its planets across its disk [6]. This paper proposes to detect a variation in the received photometric .ux from the CGARP modeling and conventional likelihood ratio based detectors. The paper is organized as follows: • Section 2 studies the theoretical model for photometric signals. The concept of CGARP is introduced from the asymptotic distribution, of the data. • Section 3 studies the maximum likelihood estimation (MLE) of the CGARP parameters and the corresponding Cram2er Rao lower bounds (CRLB’s) for large sample size. These bounds are compared to the mean square errors (MSE’s) of the estimates obtained through Monte Carlo simulations. • Finally, the Neyman–Pearson detector (NPD) for the abrupt change (AC) detection is developed in Section 4. The exact distribution of the test statistic is obtained, allowing computation of Receiver Operating Characteristics (ROC). The practical application where the AC parameters are unknown is then investigated. The generalized likelihood ratio detector (GLRD) is derived and its performance is studied from Monte Carlo simulations.
2. Signal model derivation We assume that the signal is dominated by the photon noise i.e. the read-out noise and the thermal noise for the electronic are negligible. The model derivation relies on the semiclassical theory of photodetection: the light propagation to the detector is described by diNraction theory and the photocount only occurs during the sensor photodetection. When light with a >xed intensity over time is incident on a photodetector, the joint probability of registering the successive photocounts x n ; n = 1; : : : ; N on a single pixel is distributed according to an i.i.d. Poisson distribution: Pideal (X ) =
N e− x n n=1
x n!
;
(1)
where X = (x1 ; : : : ; xN )t and is the light intensity (square of the wave amplitude) integrated between two successive measurements [8, p. 466].
The problem is more complicated when the light wave incident on the photosurface has stochastic attributes. In this case, the distribution (1) is regarded as a conditional probability distribution and the deterministic parameter is replaced by a random variable n with mean E[n ] = [8, p. 467], [16, p. 419]. The distribution of X =(x1 ; : : : ; xN )t conditioned on =(1 ; : : : ; N )t is then an independent sequence with Poisson distribution of parameter . Consequently, the unconditioned distribution of X is ∞ N e−n nx n (2) P(X ) = p() d1 : : : dN ; x n! 0 n=1
where p() is the joint distribution of (1 ; : : : ; N ). This transform relating the photocount probability to the integrated intensity probability is often referred to as the Poisson Mandel transform of p() [20]. This transform was >rst derived in 1958 using classical arguments and rederived in 1964 using a semiclassical method [16]. It is worthy to note that Eq. (2) can also handle the case where the detector is not perfect. If the quantum eGciency of the detector varies with time [16,8], Eq. (1) can be regarded as the ideal sensor signal. As previously, the statistical properties of the observed signal in the real case are described by Eq. (2). In the monovariate case, x1 is known in the actuarial literature as a mixed Poisson process [9]. The ideal case de>ned in Eq. (1) can be obtained from (2) by choosing for p() a product of Dirac delta functions. This case corresponds to a well-stabilized single-mode laser radiation. Theoretical expressions of P(X ) can also be derived in the case of polarized thermal light [8]. Unfortunately, an a priori distribution for is generally very diGcult to choose. Moreover, an analytic expression of (2) cannot be computed for most probability density functions (pdf’s) p(). For these reasons, few studies use the stochastic model described by Eq. (2) for the detection of photometric .ux variation. This paper proposes a simple model for X using two realistic assumptions: A1 : E[n ] = 1;
A2 : var[n ]:
A1 expresses a high .ux assumption and A2 conveys the fact that the variations of the integrated light intensity and the variations of the sensor are small. Based on these assumptions, this paper proposes to model x n
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
as a Gaussian AR process subjected to the constraint E[x n ] = var[x n ] = , denoted Constraint Gaussian AR Process (CGARP). The Gaussian distribution for X is justi>ed by the following proposition. Proposition 1. If the probability distribution of X veri
CXX = C
+ IN ;
(3)
where u is a N ×1 vector of ones, IN is the N ×N identity matrix; the mean of n and C
the covariance matrix of . (2) De
= VD2 V t the eigen decomposition of C
; k () the kth order cumulant of (1 ; : : : ; k ) where = (1 ; : : : ; k ) and the standardized vector: √ Y = (D + IN )−1 V t (X − u); (4) k=2
(E[Y ] = 0; CY Y = IN ). If k () = o( ) for k ¿ 2; Y converges in distribution to the multivariate Gaussian distribution N(0; IN ) when → +∞. Proof. See Appendix A. The vector Y converges in distribution (when → +∞) to the multivariate Gaussian distribution N(0; IN ). Consequently, the distribution of Y can be well approximated for large by its asymptotic distribution (see for instance [4, p. 204]). Hence, Eq. (4) and assumption A1, A2 imply that the distribution of X can be approximated by the multivariate Gaussian distribution N(u; C
+ IN ). Eq. (3) shows that second order stationarity for implies second order stationarity for X . Moreover, assumption A2 implies that the variance of x n , which equals + var[n ], can be approximated by . Finally, parametric AR modeling for X is motivated in the stationary case by the fact that for any continuous spectral density S(f), an AR process can be found with a spectral density arbitrary close to S(f) [4, p. 132]. Classical justi>cations for this model in the stationary Gaussian context such as the Wold decomposition can be also found for example in [17].
651
Based on the previous comments, X is modeled by a pth order CGARP de>ned by xn = −
p
ak x n−k +
k=1
p
ak + e n ;
(5)
k=0
where a0 = 1 and en is an i.i.d. zero mean Gaussian sequence. The variance of en in model (5) is such that var[x n ] = . It is denoted e2 (a; ) in order to take into account its dependence toward the ak and . An analytic expression of e2 (a; ) is diGcult to obtain. However, a formal expression can be obtained by rewriting the Yule Walker equations [19] as a linear system where the unknowns are the signal covariances, c = (; c(1); : : : ; c(p))t . This leads to (A1 + A2 )c = (e2 (a; ); 0; : : : ; 0)t ; where
1 a1 A1 = a2 .. .
0 1 a1 .. .
0 ··· 0 ··· 1 ··· .. .
0 0 0 ; .. .
ap ap−1 · · · a1 1
0 0 A2 = ... 0 0
(6)
ap 0 .. : . ap 0 · · · 0 0 ··· 0 0 a1 · · · ap−1 a2 · · · ap .. .
(7)
De>ning the vector e1 = (1; 0; : : : ; 0)t , we then obtain: e2 (a; ) =
= e2 (a; 1): e1t (A1 + A2 )−1 e1
(8)
Note that e2 (a; ) can be computed using the recursive algorithm studied in [1, p. 117] for the computation of the power of a linear process. It is important to note that in the ideal case, the high .ux assumption implies that the Poisson distribution can be approximated by a Gaussian distribution with same mean and variance [5]. This particular case corresponds to ak = 0; ∀k ¿ 0 in (5). Consequently, a major eNect of the detector imperfections is to correlate the signal measurements.
652
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
3. CGARP parameter estimation 3.1. Maximum likelihood estimation This section is devoted to the maximum likelihood estimation of a and . Denote L(X |a; ) the log-likelihood function of {xp+1 ; : : : ; xN } conditioned on {x1 ; : : : ; xp }: L(X |a; ) N −p log(2 e2 (a; 1)) 2
p 2 p N 1 − 2 ak x n−k − ak : 2e (a; 1)
=−
n=p+1
k=0
k=0
(9) After dropping some constants, the maximization of (9) with respect to for a known AR parameter vector a leads to
p 2 ak 2 + e2 (a; 1)
x n i.e. E[x2n ] = + 2 . This estimator is the positive root of the following second order equation: N 1 x2n = 0: (11) 2 + − N −r n=r+1
It can be easily checked that in the uncorrelated case (ak =0; ∀k ¿ 1), Eq. (10) reduces to Eq. (11). Consequently, the second order moment based estimator does not take into account the correlations between the observed samples x n . Next section derives the CRLB’s for and a which are compared to the MLE mean square errors. 3.2. Cramer Rao lower bounds CRLB’s are convenient tools for determining the achievable accuracy of estimators. Unfortunately, the CRLB’s cannot always be obtained in a simple closed form expression. In these cases, asymptotic expressions can be used to approximate the CRLB’s for large values of the number of samples. This section derives the asymptotic CRLB’s for the parameters a and of the model (5).
k=0 N 1 − N −r
n=p+1
p
2 ak x n−k
= 0:
(10)
k=0
The two roots of this second order polynomial being obviously of opposite sign, an analytic expression ˆ of (a) is given by the positive root of (10). ˆ In order to obtain the MLE’s of (; a), (a) is replaced in (9) and the resulting criterion is maximized over a using a classical optimization algorithm. A critical point of this step is the optimization initialization. We propose to take as initial condition the estimates obtained with classical AR identi>cation algorithms (e.g. [19]) after removing the estimated mean. These results suggest the following remarks: • In the Poisson i.i.d. case, the MLE of equals the sample mean. Here, could obviously also be estimated by the sample mean (method of moments) or the sample variance. However, the MLE of has to be preferred because of its asymptotical good properties [11]. • In order to take into account the constraint E[x n ] = var[x n ] = , an estimator of could also be constructed by estimating the second order moment of
Proposition 2. The asymptotic Cramer Rao lower bounds for the parameters of a CGARP de
N →+∞
=
22 + ∇a ()t Cp−1 ∇a () e2 (a; )
−
22 + ∇a ()t Cp−1 ∇a () e2 (a; )
2
e2 (a; ); (12)
lim N:CRLB(a)
N →+∞
=[Cp−1 − Cp−1 ∇a ()∇a ()t Cp−1 ]e2 (a; );
(13)
where is given by (50); the components of ∇a () by (52) and Cp−1 is computed using the Gohberg– Semencul formula [19]. Proof. See Appendix B. The theoretical expressions of the CRLB’s derived in Eqs. (12) and (13) have been compared to the
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
653
0.7
64
62
a = 0.8 1 a = 0.6 1 a1=0.6 a =0.8
0.6
60
1
0.5
56
1
ρ(a )
−10Log(MSE), −10Log(CRLB)
58
54
MLE(a1)
52
0.4
CRLB(a ) 1
50
0.3
48
46
44 20
0.2 40
60
80
100
120
140
160
0
10
20
30
40
Fig. 1. CRLB of a1 and estimated MSE of the MLE of a1 = 0:8 for p = 1 and = 1000.
−4 −6 −8
−10Log(MSE), −10Log(CRLB)
−10 −12 −14 −16
MLE(λ) CRLB(λ)
−18 −20 −22 −24 20
40
60
80
100
120
140
160
N
Fig. 2. CRLB of and estimated MSE of the MLE of = 1000 for p = 1 and a1 = 0:8.
MSE’s of the parameter MLE’s. For this purpose, 1000 independent realizations of the signal (5) with p = 1, a1 = 0:8, = 1000 have been generated for different values of N . The parameter MLE’s have been determined for each realization and the corresponding MSE’s have been computed. A comparison between the estimated MSE’s and the CRLB’s is depicted in Figs. 1 and 2. These >gures suggest the following comments: • for large values of N , a perfect adequacy between (12), (13) and the MSE’s is observed: expressions
50
60
70
80
90
100
2
σe
N
Fig. 3. Comparison between the CRLB of a1 for the CGARP and the GARP.
(12), (13) are a good approximation for both the parameter CRLB’s and the estimate MSE’s, • for small values of N , the “loss of eGciency” of the MLE and the unveri>ed asymptotic assumption for CRLB’s result in a increasing diNerence between the two curves, • the validity domain for expressions (12), (13) is approximately N ¿ 40, which is in agreement with the results obtained in [7]. Next simulations compare the CRLB’s for the CGARP parameters corresponding to model (5) with the CRLB’s for a standard unconstraint Gaussian AR Process (GARP) with mean , AR parameter vector a and driving noise variance e2 . The later are denoted CRLBu ( ), where u stands for unconstraint and = (; e2 ; a). The unconstraint Fisher Information Matrix for is block diagonal and does not depend of . Consequently, the asymptotic CRLB’s of (e2 ; a) are given by the >rst terms of (45) and the asymptotic CRLB of is CRLBu () =
N(
2 p e
k=0
ak ) 2
:
(14)
The CRLB’s are compared for a >rst order model by means of the following ratios: CRLB(a1 ) CRLB() "(a1 ) = ; "() = : (15) CRLBu (a1 ) CRLBu () Figs. 3 and 4 represent "(a1 ) and "() for diNerent values of a1 . The bounds corresponding to the
654
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661 5
1
10
λ =22 10 5 λ =35 10 4 λ =5.6 104 λ =8.8 103
4
10 a = 0.8 1 a = 0.6 1 a =0.6 1 a =0.8 1
3
10
ρ(λ)
CRLB(λ)
0.95
0.9
2
10
1
10
0
10
0.85 0
10
20
30
40
50
60
70
80
90
100
2
σe
Fig. 4. Comparison between the CRLB of for the CGARP and the GARP.
CGARP parameters are functions of , contrarily to GARP’s parameters. For this reason, the results are plotted as functions of e2 . This choice guarantees that "(a1 ) and "() are computed for two processes having the same variance. The results suggest the following remarks: • The asymptotic CRLB’s are lower for the CGARP’s than for the GARP’s, since "(a1 ) and "() are ¡ 1. This result can be easily generalized for the regression coeGcients noticing that the >rst term of (13) is the asymptotic CRLB of the GARP coeGcients and the second a positive de>nite matrix. • When a21 tends to 1, the variances of the CGARP’s and GARP’s tend to +∞. Consequently, CRLB(a1 ) and CRLBu (a1 ) decrease to 0. Fig. 3 reveals that in this case CRLB(a1 ) decreases faster than CRLBu (a1 ). • Fig. 4 shows that the limit of "() when e2 tends to +∞ is not a function of a1 , contrarily to "(a1 ). These results prove the better identi>ability of CGARP’s with respect to unconstraint GARP’s (specially for the regression coeGcients). Indeed, it is well known that the MLE is asymptotically eGcient. In other words, the variances of ML estimates are close to the corresponding CRLB’s for large number of samples. Consequently, the qualitative behavior of CRLB’s for CGARP’s and GARP’s is similar to the qualitative behavior of CGARP and GARP parameter
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
a1
Fig. 5. CRLB of for the >rst order model (N = 1000).
estimates. Regarding the detection algorithm, the performance of the GLRD is closely related to the variance of the estimated value of 1 , which is smaller for the CGARP model than for the GARP model. Consequently, the GLRD based on CGARP modeling has to be preferred to the GLRD based on GARP modeling. Finally, Fig. 5 depicts the behavior of CRLB() as a function of a1 for diNerent values of . The objective is to evaluate the achievable precision on the star .ux estimation. The parameters of the simulation have been computed from Eddington speci>cations [6]: a sampling time of 30 s, a collecting area U an optical transof 55 cm2 , a bandwidth of 3250 A, mission of 0:923 and a detector quantum eGciency of 0:84. The diNerent star magnitudes (degrees of brightness) used in the simulation are Mv = 8; 10; 12 and 14 and the observation time has been chosen equal to 8.3 hours (N = 1000). These values correspond to a .ux of = 22 × 105 , 35 × 104 , 5:6 × 104 and 8:8 × 103 . Fig. 5 shows that CRLB() increases when increases and when a1 decreases. Eq. (14) proves that this last property is also veri>ed by the GARP’s.
4. Abrupt change detection This section is devoted to the major application of the model under scope: the detection of a decrease in the photometric .ux. In the search for terrestrial planets by occultation this photometric .ux variation
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
results from the decrease of the star diNracted intensity during the transit of a planet. This can be modeled using the notations of Section 2 as AC (at instant r) in the parameter : ∀n ∈ S0 = {1; : : : ; r}
= 0 ;
∀n ∈ S1 = {r + 1; : : : ; N } = 1 (¡ 0 ):
under hypothesis H1 can be written: L(X |H1 ) = − (N − p) log e2 (a; 1) − (r − p) log 0 − (N − r) log 1 r N 1 1 1 − 2 e2n; 0 + e2n; 1 ; e (a; 1) 0 1
(16)
The planet detection problem can then be formulated as the following AC detection problem:
H1 : S1 = ∅ (jump):
(17)
In the case of i.i.d. Poisson distributed data x n , the NPD for problem (17) can be easily derived and yields: H0 rejected if
N 1 x n ¡ ': N −r
n=p+1
(18)
en; i =
p
ak x n−k − i
k=0
4.1. The Neyman–Pearson detector After dropping the constant terms, the log-likelihood function of {xp+1 ; : : : ; xN } conditioned on {x1 ; : : : ; xp }
p
ak
:
(20)
k=0
The log likelihood function under hypothesis H0 is readily obtained from (19). After dropping the constant terms, we obtain: L(X |H0 ) = −(N − p) log e2 (a; 1) − (N − p) log 0
n=r+1
The problem is obviously more complicated when the observations x n are correlated. The next section derives the NPD for problem (17), when x n is the CGARP de>ned in (5). AC detection and estimation for AR processes have been studied for long time (see [3,14] and references therein for an overview). The new contribution here is the development of a detection scheme in the particular case where the jumps occurs on the mean and variance of a CGARP. The study is restricted to oN-line change point detection algorithm [3]. A similar problem was studied in [21], for nonzero mean AR processes multiplied by a sigmoidal function modeling a jump. However, even if this model represents a jump in the mean and the power of an AR process, it can be easily checked that it cannot handle the case where mean and variance are equal. This constraint, as it will be shown below, simpli>es substantially the test statistic.
n=r+1
(19) where
H0 : S1 = ∅ (no jump);
655
−
N 1 e2n; 0 : e2 (a; 1)0
(21)
n=p+1
Using the hypothesis 1 ¡ 0 , the NPD reduces to:
p 2 N H0 1 T= ak x n−k ?' (22) N −r H1 n=r+1
k=0
which suggests the following comments: • the Neyman Pearson test is a uniformly most powerful test, since the test statistic T does not depend on 0 and 1 , • in the constraint i.i.d. case (ak = 0, ∀k ¿ 0), T reduces to the estimated power of x n , whereas in the unconstraint i.i.d. Poisson case, T is the estimated mean (18). When a jump occurs independently on the mean and the variance of the process as in [21], the test statistic T is the diNerence between two positive de>nite quadratic forms which is generally inde>nite. The exact distribution of T is then very diGcult to study and has been approximated by a Gaussian distribution in [21]. In this paper, because of the constraint (8), T reduces to a single positive de>nite quadratic form whose distribution can be determined under hypothesis H0 and H1 : under hypothesis Hi , T is the sum
656
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661 1
1
∆ λ0/λ0=10 −4
0.9
−5
∆ λ0/λ0=6 10
0.7
0.7
0.6
0.6
0.5
0.4
0.3
0.3
0.2
0.2
0.1
0.1 0
0 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
λ =5.6 104 0
0.5
0.4
0
4
0
0.8
∆ λ0/λ0=4 10−5
Pd
Pd
0.8
λ =35 10
0.9
0
0.1
0.2
0.3
0.4
of the square of N − r independent Gaussian random variables with: √ p • mean i ( k=0 ak )= N − r, • variance e2 (a; i )=(N − r). Consequently, (N − r)T=e2 (a; i ) is distributed as a noncentral *2 distribution with N − r degrees of freedom and noncentrality parameter: p i2 ( k=0 ak )2 = (N − r) i e2 (a; i ) p i ( k=0 ak )2 = (N − r) : (23) e2 (a; 1)
t ¿ 0;
0.7
0.8
0.9
1
1 0.9
a1=0.6 0.8
a =0.3 1
a1=0.05
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Pfa
(24)
where fi (t) is a mixture of central *2 pdf’s: √ 1 fi (t) = (t= i )(N −r−2)=4 I(N −r−2)=2 ( i t)e−( 2
0.6
Fig. 7. ROC curves for the NPD as a function of 0 for a >xed value of V0 =0 (V0 =0 = 10−4 , N − r = 1000). Continuous line: p = 1, a1 = 0:2. Dashed line: p = 2, a1 = 0:2, a2 = 0:05.
Pd
Fig. 6. ROC curves for the NPD as a function of V0 =0 for a >xed value of 0 (0 = 35 × 104 ; N − r = 1000). Continuous line: p = 1, a1 = 0:2. Dashed line: p = 2, a1 = 0:2, a2 = 0:05.
The pdf of T can be expressed as
N −r N −r fi −t ; pT (t|Hi ) = 2 e (a; i ) e2 (a; i )
0.5
Pfa
Pfa
i +t)=2
;
(25)
and I, (x) is the modi>ed Bessel function of the >rst kind of order , [10]. Eq. (24) allows us to plot ROC curves, in order to evaluate the NPD performance and the in.uence of the various parameters. Figs. 6 and 7 show the ROC
Fig. 8. ROC curves for the NPD as a function of a1 (p = 1, a1 = 0:2, N − r = 1000, 0 = 5:6 × 104 , V0 =0 = 10−4 ).
curves computed for a transit of an earth-like planet in front of a magnitude 10 star, as a function of the ratio V0 =0 where V0 = 0 − 1 and 0 . As can be expected, the NPD performance increases when V0 =0 and 0 increase. These >gures also show that the qualitative behavior of the NPD is very similar for a >rst order CGARP (p = 1, continuous line) and an higher-order CGARP (p = 2, dashed line). Based on these comments, next simulations have been carried out for >rst order CGARP’s for simplicity. Fig. 8
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
studies the eNect of the signal correlation on the detection performance. The detector performance increases with the signal correlation, as it could be predicted.
657
1
a1=0.8
0.9
a1=0.6
0.8
a1=0.3
0.7
4.2. Generalized likelihood ratio detector The optimal NPD provides a reference to which suboptimal detectors can be compared. However, it requires a priori knowledge of the abrupt change parameters 1 and r (the parameters ak , k = 1; : : : ; p and 0 are assumed to be known). The GLRD has received much attention in practical applications, where the abrupt change parameters are unknown [2,12]. The GLRD is the ratio of the supremum of the likelihood function with respect to the unknown parameters under both hypothesis. The resulting likelihood ratio is compared to a suitable threshold which depends on the Probability of False Alarm (PFA). The GLRD requires the computation of the MLE’s of 1 and r. When r is known, the MLE of 1 , denoted ˆ1 (r), is de>ned as in Section 3.1 by the positive root of
p 2 ak 12 + e2 (a; 1)1 k=0 N 1 − N −r
n=r+1
p
2 ak x n−k
= 0:
(26)
k=0
In order to obtain the MLE of r, the expression of ˆ1 (r) obtained from (26) is replaced in (19) and the resulting criterion is evaluated for r ∈ {p + 1; : : : ; N − 1}. The global maximizer of this criterion de>nes the MLE of r denoted rˆ and the MLE of 1 is ˆ1 (r). ˆ The GLRD statistics is then obtained by replacing r and 1 in L(X |H1 ) − L(X |H0 ) by their MLE’s. A closed form expression of the GLRD statistic distribution is clearly diGcult to derive. Consequently, ROC curves have been computed using Monte Carlo simulations. Fig. 9 shows the ROC curve of the GLRD for p=1, a1 = 0:2, N = 6000, r = 3000, 0 = 14 × 106 (Mv = 6) and V0 =0 = 1:5 × 10−2 . Note that this >gure has been obtained from 100 independent signal realizations. A comparison between Figs. 6 and 9 shows the well-known loss of performance of the GLRD compared to the NPD. However, the GLRD shows
Pd
0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Pfa
Fig. 9. ROC curves for the GLRD. p = 1, N − r = 6000, r = 3000, 0 = 14 × 106 and V0 =0 = 1:5 × 10−2 .
satisfactory performance depending on the signal correlation which is measured here by the value of a1 . It is important to note that the parameters 0 and ak , k = 1; : : : ; p, have been assumed to be known in this section. These parameters are obviously unknown in practical applications and have to be estimated. Their estimation relies on a calibration procedure that is performed independently of the detection using the method described in Section 3.1. However, it is interesting to note that the MLE’s of parameters 0 and ak , k = 1; : : : ; p could also be included in the likelihood ratio. This might signi>cantly modify the GLRD performance.
5. Conclusions This paper studied a new model denoted CGARP (for Constraint Gaussian Autoregressive Process) model for the analysis of high .ux photometric signal. The estimation of CGARP parameters using the maximum likelihood method was addressed. A comparison between the parameter estimates and the corresponding Cram2er Rao lower bounds was provided. CGARP modeling was shown to be a useful tool for the detection of extrasolar planets by occultation. Indeed, the extrasolar planet detection problem was formulated as the detection of abrupt changes in the CGARP parameters. The Neyman–Pearson detector for this detection problem with known
658
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
parameters was studied. Closed form expressions of the Neyman–Pearson detector ROC’s were derived. The generalized likelihood ratio detector was then studied for practical applications in which the parameters are unknown. ROC’s for the generalized likelihood ratio detector were computed by using Monte Carlo simulations and provided the detection performance. The abrupt change detection problem addressed in this paper assumed that the .ux variations due to the presence of a planet were observed on the last data samples. Within the framework of extrasolar planet detection, this problem occurs when photometric signals are recorded during a short time interval. When signals are recorded during a longer time interval, the photometrix .ux decreases when the planet is in front of the star and takes its original value when the planet has moved. These .ux variations can even be observed during successive occultations. It is interesting to note that the proposed detection algorithm can be extended to this situation, by de>ning an appropriate set of indices S1 in the AC detection problem (16). In this case, the NPD test statistics T can be obtained from Eq. (22), where the summation over n is replaced by n ∈ S1 and N − r is replaced by the cardinality of S1 . Consequently, the pdf of T and the ROC curves can be computed very similarly to those obtained in this paper. Of course, the GLRD performance increases when the number of samples under hypothesis H1 increases. Consequently, the detection performance should improves, when successive planet occultations can be observed. This generalization is currently under investigation. Finally, the AC model used in this paper does not take into account the ingress=egress of the planet on the star and the limb-darkening. These eNects can be easily incorporated in the formalism proposed in the article at the price of an increased number of unknown parameters. However, it is worthy to note that they will only have a little eNect on the detection performance [22].
expectations: E[x n ] = E[E[x n =n ]] = E[n ] = ;
(A.1)
cov[xi ; xj ] = E[E[xi xj |i ; j ]] − 2
(A.2)
=
cov[i ; j ]
if i = j;
+ var[i ]
if i = j:
(A.3)
2. Asymptotic normality of Y
In order to prove the asymptotic normality of Y , we study its second characteristic function denoted j t Y .Y ( ) = log /Y ( ), where /Y ( ) = E[e ] and √ = (!1 ; : : : ; !N )t . Denote ’ = V (D + IN )−1 where IN is the N ×N identity matrix. The components of the vector ’ are ’k () =
N vkq !q √ d + q=1 q
(A.4)
using obvious notations. Note that assumption A2 and the orthogonality of V implies that
1 : (A.5) ’k () = O √ A straightforward computation yields .Y ( ) = −j
N
’k + log /X (’ ):
(A.6)
k=1
Using conditional expectations, the >rst characteristic function of X can be written as follows: N j!k j t X (e −1)k |]] = E e /X ( ) = E[E[e k=1
= E[e
N
k=1
k (ej!k −1)
(A.7)
]:
(A.8) t
Appendix Asymptotic normality of Y
1. Mean and covariance of X The >rst and second-order statistics of X can be easily computed using conditional
Consequently, by denoting W (S) = E[e S ] the moment generating function of , Eq. (A.8) reads: /X (’ ) = W (S); with S = (s1 ; : : : ; sN )t and sk = ej’k () − 1.
(A.9)
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
Assuming that satis>es some regularity conditions (see [15, p. 198]) and denoting k () the kth order cumulant of (1 ; : : : ; k ) where = (1 ; : : : ; k ), the second-moment generating function of can be expanded as log W (S) =
N 1 =1
+
1 3!
N 1 s 1 + 2 (1 ; 2 )s1 s2 2 1 ;2 =1
N
3 (1 ; 2 ; 3 )s1 s2 s3 + · · ·
1 ;2 ;3 =1
(A.10)
for ¿ " (" ¿ 0 is de>ned such that the series in Eq. (A.10) represents a function which is regular for ¿ " [15, p. 198]). The second characteristic function of Y can then be computed by replacing (A.10) in (A.6): .Y ( ) = −j
N
’1 () +
1 =1
+
N
(ej’1 () − 1)
1 =1
N 1 2 (1 ; 2 )(ej’1 () − 1) 2
(A.11)
Consequently, by substituting ej’k () by its power series expansion and using Eq. (A.5) and the hypothesis k () = o(k=2 ), k ¿ 2, we obtain: N .Y ( ) = − ’1 ()2 + o(1): 2
(A.12)
1 =1
Finally, the limit oN .Y ( ) can be computed noticing that the sum in the >rst term is the norm of ’ : lim .Y ( ) = lim −
→+∞
→+∞
=−
1 2
N
N
2
!q2 :
q=1
N(0; IN ) when → ∞. Consequently, Y converges in distribution to the multivariate Gaussian distribution N(0; IN ). Appendix B. Asymptotic CRLB’s for the CGARP parameters The unknown parameter vector for the CGARP de>ned by (5) is = (; a). Since there is a one-to-one transformation between = (; a) and = (e2 ; a), the CRLB’s for and are linked by the following relation [13]: CRLB() =
@g( ) @g( ) t ; CRLB( ) @ @
(B.1)
where @g( )=@ is the Hessian of the transformation:
@g( ) =e2 ∇a ()t = : (B.2) 0 Ip @ In this expression =e2 = 1=e2 (a; 1) is the partial derivation of with respect to e2 (according to (8)) and ∇a () is the gradient of with respect to a. B.1. Asymptotic CRLB for = (e2 ; a)
1 ;2 =1
×(ej’2 () − 1) + · · · :
659
!q2 √ (dq + )2
The coeGcients of the Fisher Information Matrix (FIM) for a Gaussian process are the sum of a term that only depends on the covariances of the process and a term that takes into account its mean [19]. The >rst term of the asymptotic FIM for is the well-known asymptotic form given for example in [7]. The mean of the process being u; u=(1; : : : ; 1)t , the coeGcients of the second term are
@ t @ @ @ −1 u CN u = ut CN−1 u ; (B.3) @8k @8l @8k @8l where CN is the order N signal covariance matrix. Consequently:
(2e4 )−1 0 FIM( ) = N 0 (e2 )−1 Cp
(A.13)
+ ut CN−1 u∇ ()∇ ()t ;
q=1
Eq. (A.13) shows that the second characteristic function of Y converges to the second characteristic function of the multivariate Gaussian distribution
where t
∇ () =
t ; ∇a () : e2
(B.4)
(B.5)
660
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
The inversion of (B.4) using the inversion lemma and products of bloc matrices yields:
2 2e2 0 CRLB( ) = e 0 Cp−1 N
e2 2 − (2 ∇a ()t Cp−1 ); N Cp−1 ∇a () N (B.6)
from (8) as follows: @(e2 (a; 1))−1 @et (A1 + A2 )−1 e1 @ = e2 = e2 1 @ak @ak @ak (B.12) = −e2 e1t (A1 + A2 )−1 ×(A1 + A2 )−1 e1 :
@(A1 + A2 ) @ak (B.13)
with N
e2 ut CN−1 u 2 1 + e ut CN−1 = N N
2 −1 2 t −1 ×u + ∇a () Cp ∇a () : e2
References
(B.7)
The inverse of CN is computed using the Gohberg– Semencul formula [19] which yields:
p 2 1 ak ut CN−1 u = 2 (N − p) e k=0
+
p q=1
p−q
2 ak
−
p q=1
k=0
p
2 ak :
k=q
(B.8) Hence ut CN−1 u 1 lim = 2 N →+∞ N e
p
2 ak
:
(B.9)
k=0
Replacing this expression in (B.7) yields: = lim
N →+∞
=
p k=0
×
(B.10)
N
p 2 2 ak 1 + ak k=0
−1
22 + ∇a ()t Cp−1 ∇a () e2
:
(B.11)
B.2. Asymptotic CRLB for = (; a) The substitution of Eqs. (B.2), (B.6) in (B.1) leads to Eqs. (13), (12). Note that ∇a () can be computed
˙ om, Introduction to Stochastic Control Theory, [1] K. AstrY Academic Press, New York, 1970. [2] M. Basseville, A. Benveniste, Detection of Abrupt Changes in Signals and Dynamical Systems, Lecture Notes in Control and Information, Springer, Berlin, 1986. [3] M. Basseville, I. Nikiforov, Detection of Abrupt Changes, Theory and Application, Prentice-Hall, Englewood CliNs, NJ, 1993. [4] P.J. Brockwell, R. Davis, Time Series: Theory and Methods, Springer, Berlin, 2nd Edition, 1990. [5] H. CramZer, Mathematical Methods of Statistics, Princeton University Press, Princeton, 1971. [6] Eddington, Assessment Study Report, ESA-SCI (2000)8, July 2000. [7] B. Friedlander, B. Porat, The Exact Cramer–Rao Bound for Gaussian Autoregressive Processes, IEEE Trans. Aerospace Elect. Systems, 25 (4) 3–8. [8] J. Goodman, Statistical Optics, Wiley Interscience, New York, 1985. [9] J. Grandell, Mixed Poisson Processes, Chapman & Hall, London, 1997. [10] N.L. Johnson, S. Kotz, N. Balakrishnan, Continuous Univariate Distributions, Vol. 2, Wiley, New York, 1994. [11] S.M. Kay, Modern Spectral Estimation, Signal Processing Series, Prentice-Hall, Englewood CliNs, NJ, 1987. [12] S.M. Kay, Fundamentals of Statistical Signal Processing, Detection Theory, Signal Processing Series, Prentice-Hall, Englewood CliNs, NJ, 1998. [13] S.M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory, Signal Processing Series, Prentice-Hall, Englewood CliNs, NJ, 1993. [14] M. Lavielle, Optimal segmentation of random processes, IEEE Trans. Signal Process. 46 (5) (1998) 1365–1373. [15] E. Lukacs, Characteristic Functions, GriGn, London, 2nd Edition, 1970. [16] C. Mehta, Progress in Optics: Theory of photoelectron counting, Vol. 8, North-Holland, Amsterdam, 1970. [17] M.B. Priestley, Spectral Analysis and Time Series, Monographs and Textbooks, 8th Edition, Academic Press, New York, 1994.
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661 [18] J. Schneider, The extrasolar planets encyclopaedia, http://www.obspm.fr/planets. [19] P. Stoica, R. Moses, Introduction to Spectral Analysis, Prentice-Hall, Englewood CliNs, NJ, 1997. [20] F. Sultani, M2ethodes d’Inversion de la Transformation de Poisson, Ph.D. Thesis, Universit2e de Nice-Sophia Antipolis, France, 1995.
661
[21] J.Y. Tourneret, Detection and estimation of abrupt changes contaminated by multiplicative Gaussian noise, EURASIP Signal Process. 68 (3) (1998) 259–270. [22] J.Y. Tourneret, A. Ferrari, A. Swami, Cramer–Rao lower bounds for abrupt changes in additive and multiplicative noise, IEEE Trans. Signal Process., to be submitted.
www.elsevier.com/locate/sigpro
Parametric modeling of photometric signals A. Ferraria; ∗ , J.Y. Tourneretb , G. Alengrina a UMR
6525 Astrophysique, Universite de Nice Sophia-Antipolis, Parc Valrose, 06108 Nice cedex 2, France b ENSEEIHT=TESA, 2 rue Camichel, BP 7122, 31071 Toulouse cedex 7, France
Received 25 February 2001; received in revised form 2 August 2001; accepted 12 September 2001
Abstract This paper studies a new model for photometric signals under high .ux assumption. Photometric signals are modeled by Gaussian autoregressive processes having the same mean and variance denoted Constraint Gaussian Autoregressive Processes (CGARP’s). The estimation of the CGARP parameters is discussed. The Cram2er Rao lower bounds for these parameters are studied and compared to the estimator mean square errors. The CGARP is intended to model the signal received by a satellite designed for extrasolar planets detection. A transit of a planet in front of a star results in an abrupt change in the mean and variance of the CGARP. The Neyman–Pearson detector for this changepoint detection problem is derived when the abrupt change parameters are known. Closed form expressions for the Receiver Operating Characteristics (ROC) are provided. The Neyman–Pearson detector combined with the maximum likelihood estimator for CGARP parameters allows to study the generalized likelihood ratio detector. ROC curves are then determined using computer simulations. ? 2002 Elsevier Science B.V. All rights reserved. Keywords: Abrupt change detection; Photometric signal; Extra solar planet detection
1. Introduction The detection of extrasolar planets is a challenging problem in astronomy (see the extrasolar planets encyclopedia [18]). Among the methods currently pursued to detect extrasolar planets, the transit method (also referred as photometric method or occultation method) may be the only one to >nd earth class planets in the near future. The transit method is based on the detection of photometric .ux variations which results from the transit of a planet in front of a star. Several space-based projects have been proposed to achieve this goal including the US Kepler project, the Eddington european project and the COROT (COnvection ∗ Corresponding author. Tel.: +33-9207-6349; fax: +33-92076321. E-mail address: [email protected] (A. Ferrari).
and ROTation) french project. For a star like the Sun, the typical relative .ux variation for a Jupiter-like giant planet is 10−2 (during 25 h) whereas for a Earth-like planet it is 10−4 (during 13 h). Consequently, the photometric .ux variations caused by terrestrial planets are diGcult to detect using the conventional transit method and new detectors have to be investigated. The >rst problem addressed in this paper is the modeling of the photometric signals received by the satellite. Based on two realistic assumptions, the photometric signals are modeled by Gaussian autoregressive processes having the same mean and variance, denoted Constraint Gaussian Autoregressive Processes (CGARP’s). The second problem addressed in this paper is the detection of extrasolar planets by using standard estimation and detection tools. The theoretical distribution of photometric signals being
0165-1684/02/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 1 ) 0 0 2 1 2 - 2
650
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
generally intractable, most extrasolar planet detectors reduce to measure the brightness drop of a star which results from the transit of one of its planets across its disk [6]. This paper proposes to detect a variation in the received photometric .ux from the CGARP modeling and conventional likelihood ratio based detectors. The paper is organized as follows: • Section 2 studies the theoretical model for photometric signals. The concept of CGARP is introduced from the asymptotic distribution, of the data. • Section 3 studies the maximum likelihood estimation (MLE) of the CGARP parameters and the corresponding Cram2er Rao lower bounds (CRLB’s) for large sample size. These bounds are compared to the mean square errors (MSE’s) of the estimates obtained through Monte Carlo simulations. • Finally, the Neyman–Pearson detector (NPD) for the abrupt change (AC) detection is developed in Section 4. The exact distribution of the test statistic is obtained, allowing computation of Receiver Operating Characteristics (ROC). The practical application where the AC parameters are unknown is then investigated. The generalized likelihood ratio detector (GLRD) is derived and its performance is studied from Monte Carlo simulations.
2. Signal model derivation We assume that the signal is dominated by the photon noise i.e. the read-out noise and the thermal noise for the electronic are negligible. The model derivation relies on the semiclassical theory of photodetection: the light propagation to the detector is described by diNraction theory and the photocount only occurs during the sensor photodetection. When light with a >xed intensity over time is incident on a photodetector, the joint probability of registering the successive photocounts x n ; n = 1; : : : ; N on a single pixel is distributed according to an i.i.d. Poisson distribution: Pideal (X ) =
N e− x n n=1
x n!
;
(1)
where X = (x1 ; : : : ; xN )t and is the light intensity (square of the wave amplitude) integrated between two successive measurements [8, p. 466].
The problem is more complicated when the light wave incident on the photosurface has stochastic attributes. In this case, the distribution (1) is regarded as a conditional probability distribution and the deterministic parameter is replaced by a random variable n with mean E[n ] = [8, p. 467], [16, p. 419]. The distribution of X =(x1 ; : : : ; xN )t conditioned on =(1 ; : : : ; N )t is then an independent sequence with Poisson distribution of parameter . Consequently, the unconditioned distribution of X is ∞ N e−n nx n (2) P(X ) = p() d1 : : : dN ; x n! 0 n=1
where p() is the joint distribution of (1 ; : : : ; N ). This transform relating the photocount probability to the integrated intensity probability is often referred to as the Poisson Mandel transform of p() [20]. This transform was >rst derived in 1958 using classical arguments and rederived in 1964 using a semiclassical method [16]. It is worthy to note that Eq. (2) can also handle the case where the detector is not perfect. If the quantum eGciency of the detector varies with time [16,8], Eq. (1) can be regarded as the ideal sensor signal. As previously, the statistical properties of the observed signal in the real case are described by Eq. (2). In the monovariate case, x1 is known in the actuarial literature as a mixed Poisson process [9]. The ideal case de>ned in Eq. (1) can be obtained from (2) by choosing for p() a product of Dirac delta functions. This case corresponds to a well-stabilized single-mode laser radiation. Theoretical expressions of P(X ) can also be derived in the case of polarized thermal light [8]. Unfortunately, an a priori distribution for is generally very diGcult to choose. Moreover, an analytic expression of (2) cannot be computed for most probability density functions (pdf’s) p(). For these reasons, few studies use the stochastic model described by Eq. (2) for the detection of photometric .ux variation. This paper proposes a simple model for X using two realistic assumptions: A1 : E[n ] = 1;
A2 : var[n ]:
A1 expresses a high .ux assumption and A2 conveys the fact that the variations of the integrated light intensity and the variations of the sensor are small. Based on these assumptions, this paper proposes to model x n
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
as a Gaussian AR process subjected to the constraint E[x n ] = var[x n ] = , denoted Constraint Gaussian AR Process (CGARP). The Gaussian distribution for X is justi>ed by the following proposition. Proposition 1. If the probability distribution of X veri
CXX = C
+ IN ;
(3)
where u is a N ×1 vector of ones, IN is the N ×N identity matrix; the mean of n and C
the covariance matrix of . (2) De
= VD2 V t the eigen decomposition of C
; k () the kth order cumulant of (1 ; : : : ; k ) where = (1 ; : : : ; k ) and the standardized vector: √ Y = (D + IN )−1 V t (X − u); (4) k=2
(E[Y ] = 0; CY Y = IN ). If k () = o( ) for k ¿ 2; Y converges in distribution to the multivariate Gaussian distribution N(0; IN ) when → +∞. Proof. See Appendix A. The vector Y converges in distribution (when → +∞) to the multivariate Gaussian distribution N(0; IN ). Consequently, the distribution of Y can be well approximated for large by its asymptotic distribution (see for instance [4, p. 204]). Hence, Eq. (4) and assumption A1, A2 imply that the distribution of X can be approximated by the multivariate Gaussian distribution N(u; C
+ IN ). Eq. (3) shows that second order stationarity for implies second order stationarity for X . Moreover, assumption A2 implies that the variance of x n , which equals + var[n ], can be approximated by . Finally, parametric AR modeling for X is motivated in the stationary case by the fact that for any continuous spectral density S(f), an AR process can be found with a spectral density arbitrary close to S(f) [4, p. 132]. Classical justi>cations for this model in the stationary Gaussian context such as the Wold decomposition can be also found for example in [17].
651
Based on the previous comments, X is modeled by a pth order CGARP de>ned by xn = −
p
ak x n−k +
k=1
p
ak + e n ;
(5)
k=0
where a0 = 1 and en is an i.i.d. zero mean Gaussian sequence. The variance of en in model (5) is such that var[x n ] = . It is denoted e2 (a; ) in order to take into account its dependence toward the ak and . An analytic expression of e2 (a; ) is diGcult to obtain. However, a formal expression can be obtained by rewriting the Yule Walker equations [19] as a linear system where the unknowns are the signal covariances, c = (; c(1); : : : ; c(p))t . This leads to (A1 + A2 )c = (e2 (a; ); 0; : : : ; 0)t ; where
1 a1 A1 = a2 .. .
0 1 a1 .. .
0 ··· 0 ··· 1 ··· .. .
0 0 0 ; .. .
ap ap−1 · · · a1 1
0 0 A2 = ... 0 0
(6)
ap 0 .. : . ap 0 · · · 0 0 ··· 0 0 a1 · · · ap−1 a2 · · · ap .. .
(7)
De>ning the vector e1 = (1; 0; : : : ; 0)t , we then obtain: e2 (a; ) =
= e2 (a; 1): e1t (A1 + A2 )−1 e1
(8)
Note that e2 (a; ) can be computed using the recursive algorithm studied in [1, p. 117] for the computation of the power of a linear process. It is important to note that in the ideal case, the high .ux assumption implies that the Poisson distribution can be approximated by a Gaussian distribution with same mean and variance [5]. This particular case corresponds to ak = 0; ∀k ¿ 0 in (5). Consequently, a major eNect of the detector imperfections is to correlate the signal measurements.
652
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
3. CGARP parameter estimation 3.1. Maximum likelihood estimation This section is devoted to the maximum likelihood estimation of a and . Denote L(X |a; ) the log-likelihood function of {xp+1 ; : : : ; xN } conditioned on {x1 ; : : : ; xp }: L(X |a; ) N −p log(2 e2 (a; 1)) 2
p 2 p N 1 − 2 ak x n−k − ak : 2e (a; 1)
=−
n=p+1
k=0
k=0
(9) After dropping some constants, the maximization of (9) with respect to for a known AR parameter vector a leads to
p 2 ak 2 + e2 (a; 1)
x n i.e. E[x2n ] = + 2 . This estimator is the positive root of the following second order equation: N 1 x2n = 0: (11) 2 + − N −r n=r+1
It can be easily checked that in the uncorrelated case (ak =0; ∀k ¿ 1), Eq. (10) reduces to Eq. (11). Consequently, the second order moment based estimator does not take into account the correlations between the observed samples x n . Next section derives the CRLB’s for and a which are compared to the MLE mean square errors. 3.2. Cramer Rao lower bounds CRLB’s are convenient tools for determining the achievable accuracy of estimators. Unfortunately, the CRLB’s cannot always be obtained in a simple closed form expression. In these cases, asymptotic expressions can be used to approximate the CRLB’s for large values of the number of samples. This section derives the asymptotic CRLB’s for the parameters a and of the model (5).
k=0 N 1 − N −r
n=p+1
p
2 ak x n−k
= 0:
(10)
k=0
The two roots of this second order polynomial being obviously of opposite sign, an analytic expression ˆ of (a) is given by the positive root of (10). ˆ In order to obtain the MLE’s of (; a), (a) is replaced in (9) and the resulting criterion is maximized over a using a classical optimization algorithm. A critical point of this step is the optimization initialization. We propose to take as initial condition the estimates obtained with classical AR identi>cation algorithms (e.g. [19]) after removing the estimated mean. These results suggest the following remarks: • In the Poisson i.i.d. case, the MLE of equals the sample mean. Here, could obviously also be estimated by the sample mean (method of moments) or the sample variance. However, the MLE of has to be preferred because of its asymptotical good properties [11]. • In order to take into account the constraint E[x n ] = var[x n ] = , an estimator of could also be constructed by estimating the second order moment of
Proposition 2. The asymptotic Cramer Rao lower bounds for the parameters of a CGARP de
N →+∞
=
22 + ∇a ()t Cp−1 ∇a () e2 (a; )
−
22 + ∇a ()t Cp−1 ∇a () e2 (a; )
2
e2 (a; ); (12)
lim N:CRLB(a)
N →+∞
=[Cp−1 − Cp−1 ∇a ()∇a ()t Cp−1 ]e2 (a; );
(13)
where is given by (50); the components of ∇a () by (52) and Cp−1 is computed using the Gohberg– Semencul formula [19]. Proof. See Appendix B. The theoretical expressions of the CRLB’s derived in Eqs. (12) and (13) have been compared to the
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
653
0.7
64
62
a = 0.8 1 a = 0.6 1 a1=0.6 a =0.8
0.6
60
1
0.5
56
1
ρ(a )
−10Log(MSE), −10Log(CRLB)
58
54
MLE(a1)
52
0.4
CRLB(a ) 1
50
0.3
48
46
44 20
0.2 40
60
80
100
120
140
160
0
10
20
30
40
Fig. 1. CRLB of a1 and estimated MSE of the MLE of a1 = 0:8 for p = 1 and = 1000.
−4 −6 −8
−10Log(MSE), −10Log(CRLB)
−10 −12 −14 −16
MLE(λ) CRLB(λ)
−18 −20 −22 −24 20
40
60
80
100
120
140
160
N
Fig. 2. CRLB of and estimated MSE of the MLE of = 1000 for p = 1 and a1 = 0:8.
MSE’s of the parameter MLE’s. For this purpose, 1000 independent realizations of the signal (5) with p = 1, a1 = 0:8, = 1000 have been generated for different values of N . The parameter MLE’s have been determined for each realization and the corresponding MSE’s have been computed. A comparison between the estimated MSE’s and the CRLB’s is depicted in Figs. 1 and 2. These >gures suggest the following comments: • for large values of N , a perfect adequacy between (12), (13) and the MSE’s is observed: expressions
50
60
70
80
90
100
2
σe
N
Fig. 3. Comparison between the CRLB of a1 for the CGARP and the GARP.
(12), (13) are a good approximation for both the parameter CRLB’s and the estimate MSE’s, • for small values of N , the “loss of eGciency” of the MLE and the unveri>ed asymptotic assumption for CRLB’s result in a increasing diNerence between the two curves, • the validity domain for expressions (12), (13) is approximately N ¿ 40, which is in agreement with the results obtained in [7]. Next simulations compare the CRLB’s for the CGARP parameters corresponding to model (5) with the CRLB’s for a standard unconstraint Gaussian AR Process (GARP) with mean , AR parameter vector a and driving noise variance e2 . The later are denoted CRLBu ( ), where u stands for unconstraint and = (; e2 ; a). The unconstraint Fisher Information Matrix for is block diagonal and does not depend of . Consequently, the asymptotic CRLB’s of (e2 ; a) are given by the >rst terms of (45) and the asymptotic CRLB of is CRLBu () =
N(
2 p e
k=0
ak ) 2
:
(14)
The CRLB’s are compared for a >rst order model by means of the following ratios: CRLB(a1 ) CRLB() "(a1 ) = ; "() = : (15) CRLBu (a1 ) CRLBu () Figs. 3 and 4 represent "(a1 ) and "() for diNerent values of a1 . The bounds corresponding to the
654
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661 5
1
10
λ =22 10 5 λ =35 10 4 λ =5.6 104 λ =8.8 103
4
10 a = 0.8 1 a = 0.6 1 a =0.6 1 a =0.8 1
3
10
ρ(λ)
CRLB(λ)
0.95
0.9
2
10
1
10
0
10
0.85 0
10
20
30
40
50
60
70
80
90
100
2
σe
Fig. 4. Comparison between the CRLB of for the CGARP and the GARP.
CGARP parameters are functions of , contrarily to GARP’s parameters. For this reason, the results are plotted as functions of e2 . This choice guarantees that "(a1 ) and "() are computed for two processes having the same variance. The results suggest the following remarks: • The asymptotic CRLB’s are lower for the CGARP’s than for the GARP’s, since "(a1 ) and "() are ¡ 1. This result can be easily generalized for the regression coeGcients noticing that the >rst term of (13) is the asymptotic CRLB of the GARP coeGcients and the second a positive de>nite matrix. • When a21 tends to 1, the variances of the CGARP’s and GARP’s tend to +∞. Consequently, CRLB(a1 ) and CRLBu (a1 ) decrease to 0. Fig. 3 reveals that in this case CRLB(a1 ) decreases faster than CRLBu (a1 ). • Fig. 4 shows that the limit of "() when e2 tends to +∞ is not a function of a1 , contrarily to "(a1 ). These results prove the better identi>ability of CGARP’s with respect to unconstraint GARP’s (specially for the regression coeGcients). Indeed, it is well known that the MLE is asymptotically eGcient. In other words, the variances of ML estimates are close to the corresponding CRLB’s for large number of samples. Consequently, the qualitative behavior of CRLB’s for CGARP’s and GARP’s is similar to the qualitative behavior of CGARP and GARP parameter
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
a1
Fig. 5. CRLB of for the >rst order model (N = 1000).
estimates. Regarding the detection algorithm, the performance of the GLRD is closely related to the variance of the estimated value of 1 , which is smaller for the CGARP model than for the GARP model. Consequently, the GLRD based on CGARP modeling has to be preferred to the GLRD based on GARP modeling. Finally, Fig. 5 depicts the behavior of CRLB() as a function of a1 for diNerent values of . The objective is to evaluate the achievable precision on the star .ux estimation. The parameters of the simulation have been computed from Eddington speci>cations [6]: a sampling time of 30 s, a collecting area U an optical transof 55 cm2 , a bandwidth of 3250 A, mission of 0:923 and a detector quantum eGciency of 0:84. The diNerent star magnitudes (degrees of brightness) used in the simulation are Mv = 8; 10; 12 and 14 and the observation time has been chosen equal to 8.3 hours (N = 1000). These values correspond to a .ux of = 22 × 105 , 35 × 104 , 5:6 × 104 and 8:8 × 103 . Fig. 5 shows that CRLB() increases when increases and when a1 decreases. Eq. (14) proves that this last property is also veri>ed by the GARP’s.
4. Abrupt change detection This section is devoted to the major application of the model under scope: the detection of a decrease in the photometric .ux. In the search for terrestrial planets by occultation this photometric .ux variation
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
results from the decrease of the star diNracted intensity during the transit of a planet. This can be modeled using the notations of Section 2 as AC (at instant r) in the parameter : ∀n ∈ S0 = {1; : : : ; r}
= 0 ;
∀n ∈ S1 = {r + 1; : : : ; N } = 1 (¡ 0 ):
under hypothesis H1 can be written: L(X |H1 ) = − (N − p) log e2 (a; 1) − (r − p) log 0 − (N − r) log 1 r N 1 1 1 − 2 e2n; 0 + e2n; 1 ; e (a; 1) 0 1
(16)
The planet detection problem can then be formulated as the following AC detection problem:
H1 : S1 = ∅ (jump):
(17)
In the case of i.i.d. Poisson distributed data x n , the NPD for problem (17) can be easily derived and yields: H0 rejected if
N 1 x n ¡ ': N −r
n=p+1
(18)
en; i =
p
ak x n−k − i
k=0
4.1. The Neyman–Pearson detector After dropping the constant terms, the log-likelihood function of {xp+1 ; : : : ; xN } conditioned on {x1 ; : : : ; xp }
p
ak
:
(20)
k=0
The log likelihood function under hypothesis H0 is readily obtained from (19). After dropping the constant terms, we obtain: L(X |H0 ) = −(N − p) log e2 (a; 1) − (N − p) log 0
n=r+1
The problem is obviously more complicated when the observations x n are correlated. The next section derives the NPD for problem (17), when x n is the CGARP de>ned in (5). AC detection and estimation for AR processes have been studied for long time (see [3,14] and references therein for an overview). The new contribution here is the development of a detection scheme in the particular case where the jumps occurs on the mean and variance of a CGARP. The study is restricted to oN-line change point detection algorithm [3]. A similar problem was studied in [21], for nonzero mean AR processes multiplied by a sigmoidal function modeling a jump. However, even if this model represents a jump in the mean and the power of an AR process, it can be easily checked that it cannot handle the case where mean and variance are equal. This constraint, as it will be shown below, simpli>es substantially the test statistic.
n=r+1
(19) where
H0 : S1 = ∅ (no jump);
655
−
N 1 e2n; 0 : e2 (a; 1)0
(21)
n=p+1
Using the hypothesis 1 ¡ 0 , the NPD reduces to:
p 2 N H0 1 T= ak x n−k ?' (22) N −r H1 n=r+1
k=0
which suggests the following comments: • the Neyman Pearson test is a uniformly most powerful test, since the test statistic T does not depend on 0 and 1 , • in the constraint i.i.d. case (ak = 0, ∀k ¿ 0), T reduces to the estimated power of x n , whereas in the unconstraint i.i.d. Poisson case, T is the estimated mean (18). When a jump occurs independently on the mean and the variance of the process as in [21], the test statistic T is the diNerence between two positive de>nite quadratic forms which is generally inde>nite. The exact distribution of T is then very diGcult to study and has been approximated by a Gaussian distribution in [21]. In this paper, because of the constraint (8), T reduces to a single positive de>nite quadratic form whose distribution can be determined under hypothesis H0 and H1 : under hypothesis Hi , T is the sum
656
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661 1
1
∆ λ0/λ0=10 −4
0.9
−5
∆ λ0/λ0=6 10
0.7
0.7
0.6
0.6
0.5
0.4
0.3
0.3
0.2
0.2
0.1
0.1 0
0 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
λ =5.6 104 0
0.5
0.4
0
4
0
0.8
∆ λ0/λ0=4 10−5
Pd
Pd
0.8
λ =35 10
0.9
0
0.1
0.2
0.3
0.4
of the square of N − r independent Gaussian random variables with: √ p • mean i ( k=0 ak )= N − r, • variance e2 (a; i )=(N − r). Consequently, (N − r)T=e2 (a; i ) is distributed as a noncentral *2 distribution with N − r degrees of freedom and noncentrality parameter: p i2 ( k=0 ak )2 = (N − r) i e2 (a; i ) p i ( k=0 ak )2 = (N − r) : (23) e2 (a; 1)
t ¿ 0;
0.7
0.8
0.9
1
1 0.9
a1=0.6 0.8
a =0.3 1
a1=0.05
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Pfa
(24)
where fi (t) is a mixture of central *2 pdf’s: √ 1 fi (t) = (t= i )(N −r−2)=4 I(N −r−2)=2 ( i t)e−( 2
0.6
Fig. 7. ROC curves for the NPD as a function of 0 for a >xed value of V0 =0 (V0 =0 = 10−4 , N − r = 1000). Continuous line: p = 1, a1 = 0:2. Dashed line: p = 2, a1 = 0:2, a2 = 0:05.
Pd
Fig. 6. ROC curves for the NPD as a function of V0 =0 for a >xed value of 0 (0 = 35 × 104 ; N − r = 1000). Continuous line: p = 1, a1 = 0:2. Dashed line: p = 2, a1 = 0:2, a2 = 0:05.
The pdf of T can be expressed as
N −r N −r fi −t ; pT (t|Hi ) = 2 e (a; i ) e2 (a; i )
0.5
Pfa
Pfa
i +t)=2
;
(25)
and I, (x) is the modi>ed Bessel function of the >rst kind of order , [10]. Eq. (24) allows us to plot ROC curves, in order to evaluate the NPD performance and the in.uence of the various parameters. Figs. 6 and 7 show the ROC
Fig. 8. ROC curves for the NPD as a function of a1 (p = 1, a1 = 0:2, N − r = 1000, 0 = 5:6 × 104 , V0 =0 = 10−4 ).
curves computed for a transit of an earth-like planet in front of a magnitude 10 star, as a function of the ratio V0 =0 where V0 = 0 − 1 and 0 . As can be expected, the NPD performance increases when V0 =0 and 0 increase. These >gures also show that the qualitative behavior of the NPD is very similar for a >rst order CGARP (p = 1, continuous line) and an higher-order CGARP (p = 2, dashed line). Based on these comments, next simulations have been carried out for >rst order CGARP’s for simplicity. Fig. 8
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
studies the eNect of the signal correlation on the detection performance. The detector performance increases with the signal correlation, as it could be predicted.
657
1
a1=0.8
0.9
a1=0.6
0.8
a1=0.3
0.7
4.2. Generalized likelihood ratio detector The optimal NPD provides a reference to which suboptimal detectors can be compared. However, it requires a priori knowledge of the abrupt change parameters 1 and r (the parameters ak , k = 1; : : : ; p and 0 are assumed to be known). The GLRD has received much attention in practical applications, where the abrupt change parameters are unknown [2,12]. The GLRD is the ratio of the supremum of the likelihood function with respect to the unknown parameters under both hypothesis. The resulting likelihood ratio is compared to a suitable threshold which depends on the Probability of False Alarm (PFA). The GLRD requires the computation of the MLE’s of 1 and r. When r is known, the MLE of 1 , denoted ˆ1 (r), is de>ned as in Section 3.1 by the positive root of
p 2 ak 12 + e2 (a; 1)1 k=0 N 1 − N −r
n=r+1
p
2 ak x n−k
= 0:
(26)
k=0
In order to obtain the MLE of r, the expression of ˆ1 (r) obtained from (26) is replaced in (19) and the resulting criterion is evaluated for r ∈ {p + 1; : : : ; N − 1}. The global maximizer of this criterion de>nes the MLE of r denoted rˆ and the MLE of 1 is ˆ1 (r). ˆ The GLRD statistics is then obtained by replacing r and 1 in L(X |H1 ) − L(X |H0 ) by their MLE’s. A closed form expression of the GLRD statistic distribution is clearly diGcult to derive. Consequently, ROC curves have been computed using Monte Carlo simulations. Fig. 9 shows the ROC curve of the GLRD for p=1, a1 = 0:2, N = 6000, r = 3000, 0 = 14 × 106 (Mv = 6) and V0 =0 = 1:5 × 10−2 . Note that this >gure has been obtained from 100 independent signal realizations. A comparison between Figs. 6 and 9 shows the well-known loss of performance of the GLRD compared to the NPD. However, the GLRD shows
Pd
0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Pfa
Fig. 9. ROC curves for the GLRD. p = 1, N − r = 6000, r = 3000, 0 = 14 × 106 and V0 =0 = 1:5 × 10−2 .
satisfactory performance depending on the signal correlation which is measured here by the value of a1 . It is important to note that the parameters 0 and ak , k = 1; : : : ; p, have been assumed to be known in this section. These parameters are obviously unknown in practical applications and have to be estimated. Their estimation relies on a calibration procedure that is performed independently of the detection using the method described in Section 3.1. However, it is interesting to note that the MLE’s of parameters 0 and ak , k = 1; : : : ; p could also be included in the likelihood ratio. This might signi>cantly modify the GLRD performance.
5. Conclusions This paper studied a new model denoted CGARP (for Constraint Gaussian Autoregressive Process) model for the analysis of high .ux photometric signal. The estimation of CGARP parameters using the maximum likelihood method was addressed. A comparison between the parameter estimates and the corresponding Cram2er Rao lower bounds was provided. CGARP modeling was shown to be a useful tool for the detection of extrasolar planets by occultation. Indeed, the extrasolar planet detection problem was formulated as the detection of abrupt changes in the CGARP parameters. The Neyman–Pearson detector for this detection problem with known
658
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
parameters was studied. Closed form expressions of the Neyman–Pearson detector ROC’s were derived. The generalized likelihood ratio detector was then studied for practical applications in which the parameters are unknown. ROC’s for the generalized likelihood ratio detector were computed by using Monte Carlo simulations and provided the detection performance. The abrupt change detection problem addressed in this paper assumed that the .ux variations due to the presence of a planet were observed on the last data samples. Within the framework of extrasolar planet detection, this problem occurs when photometric signals are recorded during a short time interval. When signals are recorded during a longer time interval, the photometrix .ux decreases when the planet is in front of the star and takes its original value when the planet has moved. These .ux variations can even be observed during successive occultations. It is interesting to note that the proposed detection algorithm can be extended to this situation, by de>ning an appropriate set of indices S1 in the AC detection problem (16). In this case, the NPD test statistics T can be obtained from Eq. (22), where the summation over n is replaced by n ∈ S1 and N − r is replaced by the cardinality of S1 . Consequently, the pdf of T and the ROC curves can be computed very similarly to those obtained in this paper. Of course, the GLRD performance increases when the number of samples under hypothesis H1 increases. Consequently, the detection performance should improves, when successive planet occultations can be observed. This generalization is currently under investigation. Finally, the AC model used in this paper does not take into account the ingress=egress of the planet on the star and the limb-darkening. These eNects can be easily incorporated in the formalism proposed in the article at the price of an increased number of unknown parameters. However, it is worthy to note that they will only have a little eNect on the detection performance [22].
expectations: E[x n ] = E[E[x n =n ]] = E[n ] = ;
(A.1)
cov[xi ; xj ] = E[E[xi xj |i ; j ]] − 2
(A.2)
=
cov[i ; j ]
if i = j;
+ var[i ]
if i = j:
(A.3)
2. Asymptotic normality of Y
In order to prove the asymptotic normality of Y , we study its second characteristic function denoted j t Y .Y ( ) = log /Y ( ), where /Y ( ) = E[e ] and √ = (!1 ; : : : ; !N )t . Denote ’ = V (D + IN )−1 where IN is the N ×N identity matrix. The components of the vector ’ are ’k () =
N vkq !q √ d + q=1 q
(A.4)
using obvious notations. Note that assumption A2 and the orthogonality of V implies that
1 : (A.5) ’k () = O √ A straightforward computation yields .Y ( ) = −j
N
’k + log /X (’ ):
(A.6)
k=1
Using conditional expectations, the >rst characteristic function of X can be written as follows: N j!k j t X (e −1)k |]] = E e /X ( ) = E[E[e k=1
= E[e
N
k=1
k (ej!k −1)
(A.7)
]:
(A.8) t
Appendix Asymptotic normality of Y
1. Mean and covariance of X The >rst and second-order statistics of X can be easily computed using conditional
Consequently, by denoting W (S) = E[e S ] the moment generating function of , Eq. (A.8) reads: /X (’ ) = W (S); with S = (s1 ; : : : ; sN )t and sk = ej’k () − 1.
(A.9)
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
Assuming that satis>es some regularity conditions (see [15, p. 198]) and denoting k () the kth order cumulant of (1 ; : : : ; k ) where = (1 ; : : : ; k ), the second-moment generating function of can be expanded as log W (S) =
N 1 =1
+
1 3!
N 1 s 1 + 2 (1 ; 2 )s1 s2 2 1 ;2 =1
N
3 (1 ; 2 ; 3 )s1 s2 s3 + · · ·
1 ;2 ;3 =1
(A.10)
for ¿ " (" ¿ 0 is de>ned such that the series in Eq. (A.10) represents a function which is regular for ¿ " [15, p. 198]). The second characteristic function of Y can then be computed by replacing (A.10) in (A.6): .Y ( ) = −j
N
’1 () +
1 =1
+
N
(ej’1 () − 1)
1 =1
N 1 2 (1 ; 2 )(ej’1 () − 1) 2
(A.11)
Consequently, by substituting ej’k () by its power series expansion and using Eq. (A.5) and the hypothesis k () = o(k=2 ), k ¿ 2, we obtain: N .Y ( ) = − ’1 ()2 + o(1): 2
(A.12)
1 =1
Finally, the limit oN .Y ( ) can be computed noticing that the sum in the >rst term is the norm of ’ : lim .Y ( ) = lim −
→+∞
→+∞
=−
1 2
N
N
2
!q2 :
q=1
N(0; IN ) when → ∞. Consequently, Y converges in distribution to the multivariate Gaussian distribution N(0; IN ). Appendix B. Asymptotic CRLB’s for the CGARP parameters The unknown parameter vector for the CGARP de>ned by (5) is = (; a). Since there is a one-to-one transformation between = (; a) and = (e2 ; a), the CRLB’s for and are linked by the following relation [13]: CRLB() =
@g( ) @g( ) t ; CRLB( ) @ @
(B.1)
where @g( )=@ is the Hessian of the transformation:
@g( ) =e2 ∇a ()t = : (B.2) 0 Ip @ In this expression =e2 = 1=e2 (a; 1) is the partial derivation of with respect to e2 (according to (8)) and ∇a () is the gradient of with respect to a. B.1. Asymptotic CRLB for = (e2 ; a)
1 ;2 =1
×(ej’2 () − 1) + · · · :
659
!q2 √ (dq + )2
The coeGcients of the Fisher Information Matrix (FIM) for a Gaussian process are the sum of a term that only depends on the covariances of the process and a term that takes into account its mean [19]. The >rst term of the asymptotic FIM for is the well-known asymptotic form given for example in [7]. The mean of the process being u; u=(1; : : : ; 1)t , the coeGcients of the second term are
@ t @ @ @ −1 u CN u = ut CN−1 u ; (B.3) @8k @8l @8k @8l where CN is the order N signal covariance matrix. Consequently:
(2e4 )−1 0 FIM( ) = N 0 (e2 )−1 Cp
(A.13)
+ ut CN−1 u∇ ()∇ ()t ;
q=1
Eq. (A.13) shows that the second characteristic function of Y converges to the second characteristic function of the multivariate Gaussian distribution
where t
∇ () =
t ; ∇a () : e2
(B.4)
(B.5)
660
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661
The inversion of (B.4) using the inversion lemma and products of bloc matrices yields:
2 2e2 0 CRLB( ) = e 0 Cp−1 N
e2 2 − (2 ∇a ()t Cp−1 ); N Cp−1 ∇a () N (B.6)
from (8) as follows: @(e2 (a; 1))−1 @et (A1 + A2 )−1 e1 @ = e2 = e2 1 @ak @ak @ak (B.12) = −e2 e1t (A1 + A2 )−1 ×(A1 + A2 )−1 e1 :
@(A1 + A2 ) @ak (B.13)
with N
e2 ut CN−1 u 2 1 + e ut CN−1 = N N
2 −1 2 t −1 ×u + ∇a () Cp ∇a () : e2
References
(B.7)
The inverse of CN is computed using the Gohberg– Semencul formula [19] which yields:
p 2 1 ak ut CN−1 u = 2 (N − p) e k=0
+
p q=1
p−q
2 ak
−
p q=1
k=0
p
2 ak :
k=q
(B.8) Hence ut CN−1 u 1 lim = 2 N →+∞ N e
p
2 ak
:
(B.9)
k=0
Replacing this expression in (B.7) yields: = lim
N →+∞
=
p k=0
×
(B.10)
N
p 2 2 ak 1 + ak k=0
−1
22 + ∇a ()t Cp−1 ∇a () e2
:
(B.11)
B.2. Asymptotic CRLB for = (; a) The substitution of Eqs. (B.2), (B.6) in (B.1) leads to Eqs. (13), (12). Note that ∇a () can be computed
˙ om, Introduction to Stochastic Control Theory, [1] K. AstrY Academic Press, New York, 1970. [2] M. Basseville, A. Benveniste, Detection of Abrupt Changes in Signals and Dynamical Systems, Lecture Notes in Control and Information, Springer, Berlin, 1986. [3] M. Basseville, I. Nikiforov, Detection of Abrupt Changes, Theory and Application, Prentice-Hall, Englewood CliNs, NJ, 1993. [4] P.J. Brockwell, R. Davis, Time Series: Theory and Methods, Springer, Berlin, 2nd Edition, 1990. [5] H. CramZer, Mathematical Methods of Statistics, Princeton University Press, Princeton, 1971. [6] Eddington, Assessment Study Report, ESA-SCI (2000)8, July 2000. [7] B. Friedlander, B. Porat, The Exact Cramer–Rao Bound for Gaussian Autoregressive Processes, IEEE Trans. Aerospace Elect. Systems, 25 (4) 3–8. [8] J. Goodman, Statistical Optics, Wiley Interscience, New York, 1985. [9] J. Grandell, Mixed Poisson Processes, Chapman & Hall, London, 1997. [10] N.L. Johnson, S. Kotz, N. Balakrishnan, Continuous Univariate Distributions, Vol. 2, Wiley, New York, 1994. [11] S.M. Kay, Modern Spectral Estimation, Signal Processing Series, Prentice-Hall, Englewood CliNs, NJ, 1987. [12] S.M. Kay, Fundamentals of Statistical Signal Processing, Detection Theory, Signal Processing Series, Prentice-Hall, Englewood CliNs, NJ, 1998. [13] S.M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory, Signal Processing Series, Prentice-Hall, Englewood CliNs, NJ, 1993. [14] M. Lavielle, Optimal segmentation of random processes, IEEE Trans. Signal Process. 46 (5) (1998) 1365–1373. [15] E. Lukacs, Characteristic Functions, GriGn, London, 2nd Edition, 1970. [16] C. Mehta, Progress in Optics: Theory of photoelectron counting, Vol. 8, North-Holland, Amsterdam, 1970. [17] M.B. Priestley, Spectral Analysis and Time Series, Monographs and Textbooks, 8th Edition, Academic Press, New York, 1994.
A. Ferrari et al. / Signal Processing 82 (2002) 649 – 661 [18] J. Schneider, The extrasolar planets encyclopaedia, http://www.obspm.fr/planets. [19] P. Stoica, R. Moses, Introduction to Spectral Analysis, Prentice-Hall, Englewood CliNs, NJ, 1997. [20] F. Sultani, M2ethodes d’Inversion de la Transformation de Poisson, Ph.D. Thesis, Universit2e de Nice-Sophia Antipolis, France, 1995.
661
[21] J.Y. Tourneret, Detection and estimation of abrupt changes contaminated by multiplicative Gaussian noise, EURASIP Signal Process. 68 (3) (1998) 259–270. [22] J.Y. Tourneret, A. Ferrari, A. Swami, Cramer–Rao lower bounds for abrupt changes in additive and multiplicative noise, IEEE Trans. Signal Process., to be submitted.