- Email: [email protected]
Discrimination theory of physical processes
Physica A 293 (2001) 100 –114
www.elsevier.com/locate/physa
Discrimination theory of physical processes Mayumi Shingu-Yano ∗ , Fumiaki Shibata Graduate School of Humanities and Sciences and Department of Physics, Faculty of Science, Ochanomizu University, Otuka 2-1-1, Bunkyo-ku, Tokyo 112-8610, Japan Received 21 April 2000; received in revised form 31 October 2000
Abstract Discrimination problems of stochastic processes and time-evolving physical states are treated on the basis of an information theoretical method. From a set of received data we can select a suitably assumed probability distribution. The basic method is applied to discriminate colored stochastic processes of the Brownian motion and laser system with nonlinearity and 0uctuations. Dynamical state change is clearly speci1ed: This is quite di4cult by a direct observation of c 2001 the data. The methodology is further extended to include quantum mechanical systems. Elsevier Science B.V. All rights reserved. PACS: 02.50.Ey; 05.10.Gg; 05.30.−d; 05.40.Jc Keywords: Discrimination of states; State change; Entropy; Information criterion; Quantum 0uctuation
1. Introduction In treating scienti1c problems, we usually extract essential features on an object from available information, which is sometimes probabilistic in nature. Examples of the available information are a set of data of experiments and a received data in communications. On the basis of the information, we seek to 1nd the stochastic nature behind the data. Explicitly, we should determine such a probabilistic structure of the phenomena like a probability distribution. We are sometimes required to make a possible selection among several candidates. ∗
Corresponding author. Tel.: +81-3-5978-5320; fax: +81-3-5978-5326. E-mail address: [email protected] (M. Shingu-Yano).
c 2001 Elsevier Science B.V. All rights reserved. 0378-4371/01/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 0 ) 0 0 6 1 8 - X
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
101
For instance: 1. In digital communications, we detect signals in noise. On the basis of received signal data, we have to make decision on the reality of the data. Namely, we should determine whether the data contain the messages or not. 2. A similar situation also occurs in observing experimental data stochastically 0uctuating in time. From the experimental data, we must infer or make decision on the inherent stochastic process. 3. This is further related to discrimination problem of physical states. That is, through observations of a time-evolving physical system, it is necessary to determine a state of the relevant system. In conventical signal detection theories (see for instance Ref. [1]), methods in statistics are used: (i) Maximum a posteriori probability (MAP) criterion, (ii) Bayes method and (iii) Neyman–Pearson method. A common aspect to (i) – (iii) is the dependence upon a priori unknown quantities like cost functions and level of type I error. In contrast, we used a systematic method based on AIC [2,3], which is free from the above restrictions, to treat the signal detection problem [4]. In this paper, we show usefulness and applicability of the AIC method to treat the above-mentioned discrimination problems of stochastic processes and states of physical systems. We treat nonwhite (colored) stochastic processes, which play an important role in the theory of Brownian motion. We further study dynamics of laser action with nonlinear eJect. As a key quantity for a measure of the hypothesis testing, we introduce a detection probability Pd de1ned by Pd = P(D1 |H1 ) :
(1.1)
We also use an error probability, Pe = P(D1 |H0 )P0 + P(D0 |H1 )P1 :
(1.2)
In these expressions, P(Dj |Hi ) is the probability of making the decision Dj when Hi is true and Pi is the probability of 1nding Hi . These will be used in the following sections.
2. Short derivation of AIC For later convenience, we give a brief derivation of AIC. Let f(x|) be the probability density of a stochastic variable X having a characteristic parameter . On receiving a set of data {x(1) ; x(2) ; : : : ; x(N ) }, we have a logarithmic likelihood function of the form l() =
N l=1
ln f(x(l) |)
(2.1)
102
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
which behaves in an asymptotic limit like N 1 1 N →∞ ln f(x(l) |) → ln f(X |) l() = N N l=1
≡ I2 () ;
(2.2)
according to the law of large number, where we have put ∞ W (x) ln f(x|) d x ; I2 () = −∞
(2.3)
with the (true) probability density W (x) of X . The probability density f(x|) is diJerent from the true distribution function W (x). Then, we assume both functions coincide with each other only when the parameter takes the true value 0 : f(x|0 ) = W (x) :
(2.4)
Thus, there is a true quantity besides I2 (): ∞ W (x) ln W (x) d x = ln W (X ) : I1 = −∞
(2.5)
A “distance” between the two probability distributions is given by the Kullback–Leibler (KL) information I de1ned by I ≡ I1 − I2 () =
∞
−∞
(2.6)
W (x) ln
W (X ) = ln f(X |)
W (x) dx f(x|)
(2.7)
≡ −S ;
(2.8)
S being the Boltzmann relative entropy. The quantity I1 de1ned by (2.5) is independent of the data, and thus, minimizing I is equivalent to maximizing I2 () to 1nd an appropriate value of . The value of maximizing l() is the “maximum likelihood estimate” ˆ determined by @ ˆ l() = 0 : @ˆ
(2.9)
Expanding l() around = ˆ and taking the asymptotic limit (N → ∞), we 1nd 1 1 ∼ ˆ l() − (2.10) I2 (0 ) = N 2 where · is an average over the distribution of the data {x(1) ; x(2) ; : : : ; x(N ) }. ˆ expanded around 0 gives a relation Similarly, an average of I2 () 1 ˆ ∼ I2 () : = I2 (0 ) − 2N
(2.11)
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
103
From (2.10) and (2.11) we eliminate the unknown (i.e., data free) quantity I2 (0 ) to obtain ˆ −1 : ˆ = 1 l() (2.12) I2 () N Introducing a quantity ˆ − 1] ; AIC ≡ −2[l()
(2.13)
we can rewrite (2.12) as ˆ = − 1 AIC : I2 () 2N
(2.14)
ˆ with respect to the Thus, −(1=2N ) AIC is an “unbiased estimator” of I2 () averaging procedure of ·. When we have a vector parameter of k components, AIC is generalized to ˆ − k] : AIC = −2[l()
(2.15)
Consequently, we have only to minimize AIC in oder to minimize I (or maximize S). We already treated problem 1 of the previous section in the previous paper [3]. Then we turn our attention to problems 2 and 3 treating colored stochastic processes and nonlinear eJect. These are sometimes necessary to analyze experimental data. Moreover, a quantum mechanical extension of the methodology is given. In the following we separately discuss these problems. 3. Discrimination of colored stochastic processes As an important problem of 2 in Section 1, we treat two typical stochastic processes, i.e., a diJusion process and a process of a Brownian motion. These are fundamental but important stochastic processes in treating actual physical systems (see for instance Ref. [5]) and thus the best candidates for the discrimination problem. Explicitly, one of the processes is governed by a Langevin equation of the form d x(t) = R(t) ; (3.1) dt where R(t) is assumed to be a Gaussian white noise whose characteristics are completely determined by the correlation function, R(t)R(t ) = 22 (t − t ) :
(3.2)
The corresponding probability density P0 (x; t) to 1nd x(t) is determined by @ @2 P0 (x; t) = 2 2 P0 (x; t) ; @t @x which can be solved to give P0 (x; t|x0 ; 0) = √
(x − x0 )2 : exp − 42 t 42 t 1
(3.3)
(3.4)
104
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
As an alternative stochastic process, we consider a damping diJusion process characterized by a Langevin equation of the form d x(t) = −x(t) + R(t) ; (3.5) dt where R(t) is the same noise as in (3.1). The process x(t) generated by (3.5) is a colored process (the Gaussian–Marko4an process). Let P1 (x; t) be a probability density corresponding to (3.5). Then the Fokker–Planck equation equivalent to (3.5) is seen to be 2 @ @ 2 @ P1 (x; t) = (3.6) x + 2 P1 (x; t) @x @x @t and is solved to give
(x − x0 e−t )2 : exp − 2 2 (1 − e−2t ) 22 (1 − e−2t )
P1 (x; t|x0 ; 0) =
1
(3.7)
We have thus two hypotheses: H0 :
x(t)
with (3:1)
H1 :
x(t)
with (3:5) :
and With the use of (3.4) and (3.7) we 1nd the corresponding AICs:
N (ˆx(l) − xˆ 0 )2 2tˆl (ˆx(l) − xˆ 0 e−tˆl = )2 AIC0 − AIC1 = + ; ln − 2tˆl (1 − e−2tˆl = ) (1 − e−2tˆl = ) l=1 (3.8) where for later convenience, we used the scaled quantities de1ned by x ; tˆ ≡ t · ; xˆ ≡ : (3.9) We can generate a set of simulated data {x(l) } by solving (3.5) numerically with the stochastic Runge–Kutta method [6,7]. In this section, we call Pd of (1.1) the discrimination probability and is plotted in Fig. 1 as a function of . We see that Pd decreases with increasing and that Pd has a larger value when the observation time T increases for a certain . In Fig. 2 we plot as a function of the observation time T to 1nd Pd ¿ 0:8. From these 1gures we 1nd that the method is quite powerful even for the colored stochastic processes. That is, when the strength of the 0uctuations characterized by is small (weak 0uctuation), we can discriminate the processes with high accuracy even for shorter observation times. As the 0uctuations become large, we need longer observation times. It is rather surprising to see that a simple relation ≡
T ˙ ;
(3.10)
holds for a wide range of the variables. The stronger the 0uctuations, the longer the observation times. But the required observation times increase only linearly. Thus, the
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
105
Fig. 1. (a–b) Pd as a function of (≡ =) for diJerent observation time T .
AIC-based discrimination procedure on the colored stochastic processes works very well from a viewpoint of performance e4ciency. 4. Discrimination of physical states In the previous section we could discriminate the time-dependent stochastic processes by applying the AIC method. In this section, we turn our attention to the time evolution of physical quantities. This is a problem of type 3 in Section 1. Making use of the method, we want to infer a state of time-evolving system in certain time intervals. As an important example of the physical systems, we treat laser [8–10]. This system is composed of atoms and electromagnetic 1elds. In addition, the atoms are
106
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
Fig. 2. versus observation time T to 1nd Pd ¿ 0:8.
stimulated by an external pumping 1eld, and constituent subsystems (atoms and 1elds) are surrounded by their respective reservoirs which give rise to damping and 0uctuations. Thus, we have to solve a complex problem of the interacting composite physical system. In an ordinary laser theory, however, a simpli1ed method is employed. That is, the fast decaying atomic variables are adiabatically eliminated and we are only left with the 1eld variables with nonlinearity and 0uctuations. If the strength of the pumping 1eld exceeds a threshold value, laser begins to oscillate with an increasing amplitude, i.e., laser is in a transient state. With a further lapse of time, the amplitude attains to saturate, i.e., laser is in a stationary state. We will examine a discrimination problem of the laser action between the transient state and the stationary state. In a phenomenological laser theory, after the adiabatic elimination mentioned above, we have the following Langevin equations for the amplitude b(t) of the laser 1eld [11]: dbi (t) √ (4.1) = ! d − (b21 (t) + b22 (t)) bi (t) + q#i (t) (i = 1; 2) ; dt where b1 (t) and b2 (t) are the real and the imaginary parts of b(t) and the random noise #i (t) describes eJects due to the spontaneous emission and is characterized by #i (t) = 0 and #i (t)#j (t ) = 2 i; j (t − t ) :
(4.2) √
The parameter d represents the pumping eJect and q the noise strength. It is convenient to introduce the intensity I (t) and the phase (t) of b(t): I (t) = b∗ (t)b(t) = b21 (t) + b22 (t)
(4.3)
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
and
b2 (t) (t) = arctan b1 (t)
107
:
(4.4)
Furthermore, for later convenience, we introduce the following scaled variables: 1 ! ! ! 4 a≡ bi (t); tˆ ≡ !qt ; d; Iˆ(t) ≡ I (t); bˆ i (t) ≡ q q q (4.5) where a plays a role of the pumping parameter. Then the Langevin equations (4.1) become d bˆ i (tˆ) = (a − Iˆ(tˆ))bˆ i (tˆ) + #ˆ i (tˆ) (i = 1; 2) d tˆ and the corresponding Fokker–Planck equation is of @ @ @2 ˆ W (I ; ; tˆ) = [2(Iˆ − a)Iˆ − 4] + 2 4Iˆ + @tˆ @Iˆ @Iˆ
(4.6) the form 1 @2 W (Iˆ; ; tˆ) : Iˆ @ 2
The stationary solution of (4.6) is given by √ 1 Wst = 4 exp − (Iˆ − a)2 : 4
(4.7)
(4.8)
In contrast, a transient solution for large pump parameter (a¿20) can be obtained as follows [11]. Initially, the laser 1eld grows from a noise and the nonlinear saturating eJect can be negligible: d bˆ i (tˆ) = abˆ i (tˆ) + #ˆ i (tˆ) : d tˆ
(4.9)
This equation is valid until the time tˆc where the nonlinear term becomes appreciable. When the pump parameter is large and for tˆ¿tˆc , the nonlinear term is more eJective than the noise: d bˆ i (tˆ) = (a − Iˆ(tˆ))bˆ i (t) : (4.10) d tˆ Connecting the solutions of (4.9) and (4.10) at tˆc , we obtain a transient solution of the form W (Iˆ; tˆ) =
ae−2a(tˆ−tˆc ) 2(e2atˆc − 1)[1 + (e−2a(tˆ−tˆc ) − 1)Iˆ=a]2 aIˆe−2a(tˆ−tˆc ) : × exp − 2(e2atˆc − 1)[1 + (e−2a(tˆ−tˆc ) − 1)Iˆ=a]
(4.11)
Our aim is to decide whether the laser system is in the transient (dynamically timeevolving) state or in the stationary state when we have a set of data.
108
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
Fig. 3. Distribution of sample data at various tˆ for a = 100.
We have thus two hypotheses: H0 :
Iˆ determined by (4:8)
H1 :
Iˆ determined by (4:11) :
and
From (4.8) and (4.11) we calculate the corresponding AICs as in the previous section. The data set is generated from the nonlinear Langevin equations (4.6) with the use of the stochastic Runge–Kutta method. We show in Fig. 3 a set of simulation data for the probability density W from the transient state to the stationary state. An essential
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
109
Fig. 4. Ptr as a function of a pump parameter a for various observation time T .
feature inherent in the present problem is that the true distribution changes when the system evolves in time as is seen in Fig. 3. Namely, the hypothesis H1 may be true for the transient time region while H0 becomes true for a further lapse of time. We show in Fig. 4 the probability Ptr of selecting H1 (i.e., Pd of (1.1)) for an early stage of the time evolution (T = 0:1 in the 1gure) as a function of a, where the system is expected to be in the transient state. Re0ecting the applicability of the solution (4.11), Ptr becomes almost unity for a¿20. It is rather surprising to see that Ptr ¿0:8 even for a ¡ 20. With increasing value of T , the simulated data contains information on the stationary state as well, and therefore, Ptr becomes very low (T = 0:8). Thus, we are led to study details of the early stage time evolution. In Fig. 5, we show Ptr as a function of time for a = 100. Namely, the abscissa gives a measure of the laser system to be in the transient state and thus, the hypothesis H1 is correct when Ptr = 1, whereas H0 is correct when Ptr = 0. From the 1gure, we can clearly specify the two diJerent states separated by a dynamical change in the laser system around tˆ= 0:08. In contrast, it is di4cult to decide the state of laser action when we observe the data of Fig. 3. Thus, the AIC method is quite powerful in discriminating physical states with nonlinearity and 0uctuations: We should emphasize that the dynamical state change is clearly observed by the method. 5. Generalization to quantum systems The foregoing discrimination problems are con1ned within a realm of classical probability theory. In this section, we further generalize the method to include quantum systems.
110
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
Fig. 5. Discrimination between H0 and H1 for 100 times simulation.
To make the problem precise, we treat a simple but nonetheless important physical system. This is composed of a harmonic oscillator in interaction with a reserver. The Hamiltonian of the composite system is written in the form H = HS + HB + H1 ≡ H0 + H1 ;
(5.1)
HS = ˝!0 a† a ;
(5.2)
where
HB =
j
and H1 = ˝
˝!j Bj† Bj
j
(gj Bj a† + gj∗ Bj† a) :
(5.3)
(5.4)
In these expressions, a and a† are the annihilation and the creation operators of the relevant oscillator system while Bj and Bj† are the reserver operators. Then a master equation of the reduced density matrix -(t) of the relevant oscillator system is given by (see for instance Ref. [12]) -(t) ˙ = −i!0 [a† a; -(t)] + [2a-(t)a† − -(t)a† a − a† a-(t)] † + 2n[a-(t)a T + a† -(t)a − a† a-(t) − -(t)aa† ] ;
(5.5)
where T ij Bi† Bj = n
(5.6)
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
and
|gj |2
j
∞
0
d/ exp{i(!0 − !j )/}nT ≡ nT
111
(5.7)
with nT =
1 e˝!0 =kB T
−1
:
(5.8)
In deriving (5.5) we neglected a frequency shift due to the reservoir. In order to apply the AIC method to the quantum mechanical system, we use a Q-function, or a normally ordered c-number equivalent of -(t), de1ned by [13] P1 (z; z ∗ ; t) ≡ z|-(a; a† ; t)|z
(5.9)
where |z is a coherent state satisfying a|z = z|z :
(5.10)
The function P1 (z; z ∗ ; t) is a good candidate for the probability density in the de1nition of AIC. When we use the Wigner function [14,15] instead of the Q-function, we would have negative values which give rise to di4culties in de1ning the K–L entropy. In contrast, the Q-function is equivalent to the Husimi function [16] which is obtained by smoothing procedure of the Wigner function with nonnegative values. Eq. (5.5) is transformed into the following Fokker–Planck-type equation for P1 (z; z ∗ ; t): @ @ @ P1 (z; z ∗ ; t) = ( + i!0 ) (zP1 ) + ( − i!0 ) ∗ (z ∗ P1 ) @z @z @t + 2(nT + 1)
@2 P1 : @z ∗ @z
(5.11)
The corresponding Langevin equation is of the form d z(t) = −(i!0 + )z(t) + R(t) ; dt
(5.12)
where R(t) is the white Gaussian noise, i.e., R(t) = Rx (t) + iRy (t)
(5.13)
Rx (t)Rx (t ) = Ry (t)Ry (t ) = 22 (t − t ) ;
(5.14)
with
being de1ned by 2 = (nT + 1): The solution to (5.11) is given by P1 (z; z
∗
; t|z0 ; z0∗ ; 0)
1 |z − z0 e−(i!0 +)t |2 : = exp − (nT + 1)(1 − e−2t ) (nT + 1)(1 − e−2t ) (5.15)
112
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
Fig. 6. Pd as a function of for !ˆ 0 = 0:1 and 10.
Our choice for H1 is the stochastic process z(t) determined by (5.5). As a hypothesis H0 , we take a diJusion process of a(t). In the same way as for H1 , we 1nd the corresponding diJusion equation and the Langevin equation @ @2 P0 (z; z ∗ ; t) = 2(nT + 1) ∗ P0 @t @z @z
(5.16)
and d z(t) = R(t) : dt The solution to (5.16) is given by P0 (z; z
∗
; t|z0 ; z0∗ ; 0)
(5.17)
1 |z − z0 |2 : = exp − 2(nT + 1)t 2(nT + 1)t
(5.18)
Using (5.18) and (5.15), we have the corresponding AIC0 and AIC1 :
N 2tˆl |ˆz (l) − zˆ0 |2 1 AIC0 − AIC1 = ln + ˆ 2 2tˆl (1 − e−2t l = ) l=1 |ˆz (l) − zˆ0 e−(i!ˆ 0 +1= )tˆl |2 − 1 − e−2tˆl =
;
(5.19)
where we used the scaled variables and the parameters: zˆ = z=; ! ˆ 0 = !0 =; = =, and tˆ=t ·. In Fig. 6, we show Pd as a function of and in Fig. 7 as a function of ! ˆ 0, where it is found that Pd decreases with increasing value of . When the 0uctuation speed becomes slow, the quantum mechanical damping process of the oscillator looks like the diJusion process. We should also note that the K–L information I (or the Boltzmann relative entropy S), (2.8), is de1ned in the framework of a probability theory. Among several candidates
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
113
Fig. 7. Pd as a function of !ˆ 0 for = 0:5 and 1.0.
in generalizing the AIC method to include the quantum mechanical systems, we used the Q-function. The corresponding entropy is sometimes referred to as the Wehl–Lieb entropy [17,18]. This is a well-behaved quantity and gives no di4culty while the other functions like the Wigner function and the P-function show ill behavior in calculating the entropy as mentioned earlier. This is the reason why we adopted the Q-function in generalizing AIC to the quantum systems.
6. Concluding remarks In the previous paper [3], we studied the detection of signals in a Gaussian white noise on the basis of the AIC method. In this paper, we proceeded to treat the discrimination problems of the colored stochastic processes. That is, the Gaussian–Markovian (Ornstein–Uhlenbeck) process is examined in comparison with the diJusion process. The discrimination probability Pd depends on the parameter and the observation time T . The longer the observation time T and the smaller the strength of 0uctuations , we found the higher values of Pd . This result can be intuitively expected. Moreover, we found a simple linear relation between the two requiring Pd greater than a speci1ed value (Fig. 2). The method was further applied to discriminate dynamically time-evolving states of physical systems. As a typical system with dynamical state change, we treated a problem of laser action where the nonlinearity and the 0uctuations play essential roles. The dynamical state change was clearly speci1ed when we calculate AICs. It should be emphasized that this kind of dynamical state change cannot be recognized by a simple observation of the data. Moreover, we showed that the methodology can be extended even for quantum mechanical systems when we employ the Q-function in the de1nition of the K–L entropy. Among several possibilities
114
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
in de1ning quantum mechanical entropy, use of the Q-function (Wehl–Lieb entropy) is reasonable because of its nonsingular behavior in contrast with other functions. In fact, we could apply the method to the quantum mechanical relaxation process of a damped oscillator. The problems studied in this paper, especially those of the dynamical state change and the quantum mechanical system, have not yet been treated by this method. In conclusion, we showed usefulness and applicability of the AIC method beyond a scope of statistics to include the discrimination problems in dynamically time-evolving physical systems. References [1] R.N. McDonough, A.D. Whalen, Detection of Signals in Noise, 2nd Edition, Academic Press, Sandiego, New York, Tokyo, 1995. [2] H. Akaike, in: B.N. Petrov, F. Csaki (Eds.), Proceedings of the Second International Symposium on Information Theory, Academiai Kiado, Budapest, 1973, p. 267. [3] H. Akaike, in: E. Parzen, K. Tanabe, G. Kitagawa (Eds.), Selected Papers of Hirotsugu Akaike, Springer, New York, 1998. [4] M. Shingu-Yano, A. Kobayashi, Y. Nakamura, F. Shibata, J. Grad. School Humanities Sci., Ochanomizu University, Tokyo 2 (2000) 205. [5] R. Kubo, M. Toda, N. Hashitsume, Statistical Physics II, Springer, Berlin, Heidelberg, New York, Tokyo, 1985. [6] C. Uchiyama, F. Shibata, Physica A 198 (1993) 538. [7] F. Shibata, Y. Shimoo, Physica A 215 (1995) 87. [8] H. Haken, in: L. Genzel (Ed.), Handbuch der Physik, Vol. XXV/2c, Light and Matter, Springer, Berlin, Heidelberg, New York, 1970. [9] W.H. Louisell, Quantum Statistical Properties of Radiation, Wiley, New York, London, Sydney, Tront, 1973. [10] M. Sargent III, M.O. Scully, W.E. Lamb Jr., Laser Physics, Addison-Wesley, Reading, MA, 1974. [11] H. Risken, The Fokker–Planck Equation, 2nd Edition, Springer, Berlin, Heidelberg, New York, London, Paris, Tokyo, 1989. [12] F. Haake, Springer Tract in Modern Physics, Vol. 66, Springer, Berlin, Heidelberg, New York, 1973, p. 98. [13] R.J. Glauber, Phys. Rev. 131 (1963) 2766. [14] E. Wigner, Phys. Rev. 40 (1932) 749. [15] R. Kubo, J. Phys. Soc. Jpn. 19 (1964) 2127. [16] K. Husimi, J. Phys. Math. Soc. Jpn. 22 (1940) 264. [17] A. Wehrl, Rept. Math. Phys. 16 (1979) 353. [18] E.H. Lieb, Commun. Math. Phys. 62 (1978) 35.
www.elsevier.com/locate/physa
Discrimination theory of physical processes Mayumi Shingu-Yano ∗ , Fumiaki Shibata Graduate School of Humanities and Sciences and Department of Physics, Faculty of Science, Ochanomizu University, Otuka 2-1-1, Bunkyo-ku, Tokyo 112-8610, Japan Received 21 April 2000; received in revised form 31 October 2000
Abstract Discrimination problems of stochastic processes and time-evolving physical states are treated on the basis of an information theoretical method. From a set of received data we can select a suitably assumed probability distribution. The basic method is applied to discriminate colored stochastic processes of the Brownian motion and laser system with nonlinearity and 0uctuations. Dynamical state change is clearly speci1ed: This is quite di4cult by a direct observation of c 2001 the data. The methodology is further extended to include quantum mechanical systems. Elsevier Science B.V. All rights reserved. PACS: 02.50.Ey; 05.10.Gg; 05.30.−d; 05.40.Jc Keywords: Discrimination of states; State change; Entropy; Information criterion; Quantum 0uctuation
1. Introduction In treating scienti1c problems, we usually extract essential features on an object from available information, which is sometimes probabilistic in nature. Examples of the available information are a set of data of experiments and a received data in communications. On the basis of the information, we seek to 1nd the stochastic nature behind the data. Explicitly, we should determine such a probabilistic structure of the phenomena like a probability distribution. We are sometimes required to make a possible selection among several candidates. ∗
Corresponding author. Tel.: +81-3-5978-5320; fax: +81-3-5978-5326. E-mail address: [email protected] (M. Shingu-Yano).
c 2001 Elsevier Science B.V. All rights reserved. 0378-4371/01/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 0 ) 0 0 6 1 8 - X
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
101
For instance: 1. In digital communications, we detect signals in noise. On the basis of received signal data, we have to make decision on the reality of the data. Namely, we should determine whether the data contain the messages or not. 2. A similar situation also occurs in observing experimental data stochastically 0uctuating in time. From the experimental data, we must infer or make decision on the inherent stochastic process. 3. This is further related to discrimination problem of physical states. That is, through observations of a time-evolving physical system, it is necessary to determine a state of the relevant system. In conventical signal detection theories (see for instance Ref. [1]), methods in statistics are used: (i) Maximum a posteriori probability (MAP) criterion, (ii) Bayes method and (iii) Neyman–Pearson method. A common aspect to (i) – (iii) is the dependence upon a priori unknown quantities like cost functions and level of type I error. In contrast, we used a systematic method based on AIC [2,3], which is free from the above restrictions, to treat the signal detection problem [4]. In this paper, we show usefulness and applicability of the AIC method to treat the above-mentioned discrimination problems of stochastic processes and states of physical systems. We treat nonwhite (colored) stochastic processes, which play an important role in the theory of Brownian motion. We further study dynamics of laser action with nonlinear eJect. As a key quantity for a measure of the hypothesis testing, we introduce a detection probability Pd de1ned by Pd = P(D1 |H1 ) :
(1.1)
We also use an error probability, Pe = P(D1 |H0 )P0 + P(D0 |H1 )P1 :
(1.2)
In these expressions, P(Dj |Hi ) is the probability of making the decision Dj when Hi is true and Pi is the probability of 1nding Hi . These will be used in the following sections.
2. Short derivation of AIC For later convenience, we give a brief derivation of AIC. Let f(x|) be the probability density of a stochastic variable X having a characteristic parameter . On receiving a set of data {x(1) ; x(2) ; : : : ; x(N ) }, we have a logarithmic likelihood function of the form l() =
N l=1
ln f(x(l) |)
(2.1)
102
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
which behaves in an asymptotic limit like N 1 1 N →∞ ln f(x(l) |) → ln f(X |) l() = N N l=1
≡ I2 () ;
(2.2)
according to the law of large number, where we have put ∞ W (x) ln f(x|) d x ; I2 () = −∞
(2.3)
with the (true) probability density W (x) of X . The probability density f(x|) is diJerent from the true distribution function W (x). Then, we assume both functions coincide with each other only when the parameter takes the true value 0 : f(x|0 ) = W (x) :
(2.4)
Thus, there is a true quantity besides I2 (): ∞ W (x) ln W (x) d x = ln W (X ) : I1 = −∞
(2.5)
A “distance” between the two probability distributions is given by the Kullback–Leibler (KL) information I de1ned by I ≡ I1 − I2 () =
∞
−∞
(2.6)
W (x) ln
W (X ) = ln f(X |)
W (x) dx f(x|)
(2.7)
≡ −S ;
(2.8)
S being the Boltzmann relative entropy. The quantity I1 de1ned by (2.5) is independent of the data, and thus, minimizing I is equivalent to maximizing I2 () to 1nd an appropriate value of . The value of maximizing l() is the “maximum likelihood estimate” ˆ determined by @ ˆ l() = 0 : @ˆ
(2.9)
Expanding l() around = ˆ and taking the asymptotic limit (N → ∞), we 1nd 1 1 ∼ ˆ l() − (2.10) I2 (0 ) = N 2 where · is an average over the distribution of the data {x(1) ; x(2) ; : : : ; x(N ) }. ˆ expanded around 0 gives a relation Similarly, an average of I2 () 1 ˆ ∼ I2 () : = I2 (0 ) − 2N
(2.11)
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
103
From (2.10) and (2.11) we eliminate the unknown (i.e., data free) quantity I2 (0 ) to obtain ˆ −1 : ˆ = 1 l() (2.12) I2 () N Introducing a quantity ˆ − 1] ; AIC ≡ −2[l()
(2.13)
we can rewrite (2.12) as ˆ = − 1 AIC : I2 () 2N
(2.14)
ˆ with respect to the Thus, −(1=2N ) AIC is an “unbiased estimator” of I2 () averaging procedure of ·. When we have a vector parameter of k components, AIC is generalized to ˆ − k] : AIC = −2[l()
(2.15)
Consequently, we have only to minimize AIC in oder to minimize I (or maximize S). We already treated problem 1 of the previous section in the previous paper [3]. Then we turn our attention to problems 2 and 3 treating colored stochastic processes and nonlinear eJect. These are sometimes necessary to analyze experimental data. Moreover, a quantum mechanical extension of the methodology is given. In the following we separately discuss these problems. 3. Discrimination of colored stochastic processes As an important problem of 2 in Section 1, we treat two typical stochastic processes, i.e., a diJusion process and a process of a Brownian motion. These are fundamental but important stochastic processes in treating actual physical systems (see for instance Ref. [5]) and thus the best candidates for the discrimination problem. Explicitly, one of the processes is governed by a Langevin equation of the form d x(t) = R(t) ; (3.1) dt where R(t) is assumed to be a Gaussian white noise whose characteristics are completely determined by the correlation function, R(t)R(t ) = 22 (t − t ) :
(3.2)
The corresponding probability density P0 (x; t) to 1nd x(t) is determined by @ @2 P0 (x; t) = 2 2 P0 (x; t) ; @t @x which can be solved to give P0 (x; t|x0 ; 0) = √
(x − x0 )2 : exp − 42 t 42 t 1
(3.3)
(3.4)
104
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
As an alternative stochastic process, we consider a damping diJusion process characterized by a Langevin equation of the form d x(t) = −x(t) + R(t) ; (3.5) dt where R(t) is the same noise as in (3.1). The process x(t) generated by (3.5) is a colored process (the Gaussian–Marko4an process). Let P1 (x; t) be a probability density corresponding to (3.5). Then the Fokker–Planck equation equivalent to (3.5) is seen to be 2 @ @ 2 @ P1 (x; t) = (3.6) x + 2 P1 (x; t) @x @x @t and is solved to give
(x − x0 e−t )2 : exp − 2 2 (1 − e−2t ) 22 (1 − e−2t )
P1 (x; t|x0 ; 0) =
1
(3.7)
We have thus two hypotheses: H0 :
x(t)
with (3:1)
H1 :
x(t)
with (3:5) :
and With the use of (3.4) and (3.7) we 1nd the corresponding AICs:
N (ˆx(l) − xˆ 0 )2 2tˆl (ˆx(l) − xˆ 0 e−tˆl = )2 AIC0 − AIC1 = + ; ln − 2tˆl (1 − e−2tˆl = ) (1 − e−2tˆl = ) l=1 (3.8) where for later convenience, we used the scaled quantities de1ned by x ; tˆ ≡ t · ; xˆ ≡ : (3.9) We can generate a set of simulated data {x(l) } by solving (3.5) numerically with the stochastic Runge–Kutta method [6,7]. In this section, we call Pd of (1.1) the discrimination probability and is plotted in Fig. 1 as a function of . We see that Pd decreases with increasing and that Pd has a larger value when the observation time T increases for a certain . In Fig. 2 we plot as a function of the observation time T to 1nd Pd ¿ 0:8. From these 1gures we 1nd that the method is quite powerful even for the colored stochastic processes. That is, when the strength of the 0uctuations characterized by is small (weak 0uctuation), we can discriminate the processes with high accuracy even for shorter observation times. As the 0uctuations become large, we need longer observation times. It is rather surprising to see that a simple relation ≡
T ˙ ;
(3.10)
holds for a wide range of the variables. The stronger the 0uctuations, the longer the observation times. But the required observation times increase only linearly. Thus, the
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
105
Fig. 1. (a–b) Pd as a function of (≡ =) for diJerent observation time T .
AIC-based discrimination procedure on the colored stochastic processes works very well from a viewpoint of performance e4ciency. 4. Discrimination of physical states In the previous section we could discriminate the time-dependent stochastic processes by applying the AIC method. In this section, we turn our attention to the time evolution of physical quantities. This is a problem of type 3 in Section 1. Making use of the method, we want to infer a state of time-evolving system in certain time intervals. As an important example of the physical systems, we treat laser [8–10]. This system is composed of atoms and electromagnetic 1elds. In addition, the atoms are
106
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
Fig. 2. versus observation time T to 1nd Pd ¿ 0:8.
stimulated by an external pumping 1eld, and constituent subsystems (atoms and 1elds) are surrounded by their respective reservoirs which give rise to damping and 0uctuations. Thus, we have to solve a complex problem of the interacting composite physical system. In an ordinary laser theory, however, a simpli1ed method is employed. That is, the fast decaying atomic variables are adiabatically eliminated and we are only left with the 1eld variables with nonlinearity and 0uctuations. If the strength of the pumping 1eld exceeds a threshold value, laser begins to oscillate with an increasing amplitude, i.e., laser is in a transient state. With a further lapse of time, the amplitude attains to saturate, i.e., laser is in a stationary state. We will examine a discrimination problem of the laser action between the transient state and the stationary state. In a phenomenological laser theory, after the adiabatic elimination mentioned above, we have the following Langevin equations for the amplitude b(t) of the laser 1eld [11]: dbi (t) √ (4.1) = ! d − (b21 (t) + b22 (t)) bi (t) + q#i (t) (i = 1; 2) ; dt where b1 (t) and b2 (t) are the real and the imaginary parts of b(t) and the random noise #i (t) describes eJects due to the spontaneous emission and is characterized by #i (t) = 0 and #i (t)#j (t ) = 2 i; j (t − t ) :
(4.2) √
The parameter d represents the pumping eJect and q the noise strength. It is convenient to introduce the intensity I (t) and the phase (t) of b(t): I (t) = b∗ (t)b(t) = b21 (t) + b22 (t)
(4.3)
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
and
b2 (t) (t) = arctan b1 (t)
107
:
(4.4)
Furthermore, for later convenience, we introduce the following scaled variables: 1 ! ! ! 4 a≡ bi (t); tˆ ≡ !qt ; d; Iˆ(t) ≡ I (t); bˆ i (t) ≡ q q q (4.5) where a plays a role of the pumping parameter. Then the Langevin equations (4.1) become d bˆ i (tˆ) = (a − Iˆ(tˆ))bˆ i (tˆ) + #ˆ i (tˆ) (i = 1; 2) d tˆ and the corresponding Fokker–Planck equation is of @ @ @2 ˆ W (I ; ; tˆ) = [2(Iˆ − a)Iˆ − 4] + 2 4Iˆ + @tˆ @Iˆ @Iˆ
(4.6) the form 1 @2 W (Iˆ; ; tˆ) : Iˆ @ 2
The stationary solution of (4.6) is given by √ 1 Wst = 4 exp − (Iˆ − a)2 : 4
(4.7)
(4.8)
In contrast, a transient solution for large pump parameter (a¿20) can be obtained as follows [11]. Initially, the laser 1eld grows from a noise and the nonlinear saturating eJect can be negligible: d bˆ i (tˆ) = abˆ i (tˆ) + #ˆ i (tˆ) : d tˆ
(4.9)
This equation is valid until the time tˆc where the nonlinear term becomes appreciable. When the pump parameter is large and for tˆ¿tˆc , the nonlinear term is more eJective than the noise: d bˆ i (tˆ) = (a − Iˆ(tˆ))bˆ i (t) : (4.10) d tˆ Connecting the solutions of (4.9) and (4.10) at tˆc , we obtain a transient solution of the form W (Iˆ; tˆ) =
ae−2a(tˆ−tˆc ) 2(e2atˆc − 1)[1 + (e−2a(tˆ−tˆc ) − 1)Iˆ=a]2 aIˆe−2a(tˆ−tˆc ) : × exp − 2(e2atˆc − 1)[1 + (e−2a(tˆ−tˆc ) − 1)Iˆ=a]
(4.11)
Our aim is to decide whether the laser system is in the transient (dynamically timeevolving) state or in the stationary state when we have a set of data.
108
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
Fig. 3. Distribution of sample data at various tˆ for a = 100.
We have thus two hypotheses: H0 :
Iˆ determined by (4:8)
H1 :
Iˆ determined by (4:11) :
and
From (4.8) and (4.11) we calculate the corresponding AICs as in the previous section. The data set is generated from the nonlinear Langevin equations (4.6) with the use of the stochastic Runge–Kutta method. We show in Fig. 3 a set of simulation data for the probability density W from the transient state to the stationary state. An essential
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
109
Fig. 4. Ptr as a function of a pump parameter a for various observation time T .
feature inherent in the present problem is that the true distribution changes when the system evolves in time as is seen in Fig. 3. Namely, the hypothesis H1 may be true for the transient time region while H0 becomes true for a further lapse of time. We show in Fig. 4 the probability Ptr of selecting H1 (i.e., Pd of (1.1)) for an early stage of the time evolution (T = 0:1 in the 1gure) as a function of a, where the system is expected to be in the transient state. Re0ecting the applicability of the solution (4.11), Ptr becomes almost unity for a¿20. It is rather surprising to see that Ptr ¿0:8 even for a ¡ 20. With increasing value of T , the simulated data contains information on the stationary state as well, and therefore, Ptr becomes very low (T = 0:8). Thus, we are led to study details of the early stage time evolution. In Fig. 5, we show Ptr as a function of time for a = 100. Namely, the abscissa gives a measure of the laser system to be in the transient state and thus, the hypothesis H1 is correct when Ptr = 1, whereas H0 is correct when Ptr = 0. From the 1gure, we can clearly specify the two diJerent states separated by a dynamical change in the laser system around tˆ= 0:08. In contrast, it is di4cult to decide the state of laser action when we observe the data of Fig. 3. Thus, the AIC method is quite powerful in discriminating physical states with nonlinearity and 0uctuations: We should emphasize that the dynamical state change is clearly observed by the method. 5. Generalization to quantum systems The foregoing discrimination problems are con1ned within a realm of classical probability theory. In this section, we further generalize the method to include quantum systems.
110
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
Fig. 5. Discrimination between H0 and H1 for 100 times simulation.
To make the problem precise, we treat a simple but nonetheless important physical system. This is composed of a harmonic oscillator in interaction with a reserver. The Hamiltonian of the composite system is written in the form H = HS + HB + H1 ≡ H0 + H1 ;
(5.1)
HS = ˝!0 a† a ;
(5.2)
where
HB =
j
and H1 = ˝
˝!j Bj† Bj
j
(gj Bj a† + gj∗ Bj† a) :
(5.3)
(5.4)
In these expressions, a and a† are the annihilation and the creation operators of the relevant oscillator system while Bj and Bj† are the reserver operators. Then a master equation of the reduced density matrix -(t) of the relevant oscillator system is given by (see for instance Ref. [12]) -(t) ˙ = −i!0 [a† a; -(t)] + [2a-(t)a† − -(t)a† a − a† a-(t)] † + 2n[a-(t)a T + a† -(t)a − a† a-(t) − -(t)aa† ] ;
(5.5)
where T ij Bi† Bj = n
(5.6)
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
and
|gj |2
j
∞
0
d/ exp{i(!0 − !j )/}nT ≡ nT
111
(5.7)
with nT =
1 e˝!0 =kB T
−1
:
(5.8)
In deriving (5.5) we neglected a frequency shift due to the reservoir. In order to apply the AIC method to the quantum mechanical system, we use a Q-function, or a normally ordered c-number equivalent of -(t), de1ned by [13] P1 (z; z ∗ ; t) ≡ z|-(a; a† ; t)|z
(5.9)
where |z is a coherent state satisfying a|z = z|z :
(5.10)
The function P1 (z; z ∗ ; t) is a good candidate for the probability density in the de1nition of AIC. When we use the Wigner function [14,15] instead of the Q-function, we would have negative values which give rise to di4culties in de1ning the K–L entropy. In contrast, the Q-function is equivalent to the Husimi function [16] which is obtained by smoothing procedure of the Wigner function with nonnegative values. Eq. (5.5) is transformed into the following Fokker–Planck-type equation for P1 (z; z ∗ ; t): @ @ @ P1 (z; z ∗ ; t) = ( + i!0 ) (zP1 ) + ( − i!0 ) ∗ (z ∗ P1 ) @z @z @t + 2(nT + 1)
@2 P1 : @z ∗ @z
(5.11)
The corresponding Langevin equation is of the form d z(t) = −(i!0 + )z(t) + R(t) ; dt
(5.12)
where R(t) is the white Gaussian noise, i.e., R(t) = Rx (t) + iRy (t)
(5.13)
Rx (t)Rx (t ) = Ry (t)Ry (t ) = 22 (t − t ) ;
(5.14)
with
being de1ned by 2 = (nT + 1): The solution to (5.11) is given by P1 (z; z
∗
; t|z0 ; z0∗ ; 0)
1 |z − z0 e−(i!0 +)t |2 : = exp − (nT + 1)(1 − e−2t ) (nT + 1)(1 − e−2t ) (5.15)
112
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
Fig. 6. Pd as a function of for !ˆ 0 = 0:1 and 10.
Our choice for H1 is the stochastic process z(t) determined by (5.5). As a hypothesis H0 , we take a diJusion process of a(t). In the same way as for H1 , we 1nd the corresponding diJusion equation and the Langevin equation @ @2 P0 (z; z ∗ ; t) = 2(nT + 1) ∗ P0 @t @z @z
(5.16)
and d z(t) = R(t) : dt The solution to (5.16) is given by P0 (z; z
∗
; t|z0 ; z0∗ ; 0)
(5.17)
1 |z − z0 |2 : = exp − 2(nT + 1)t 2(nT + 1)t
(5.18)
Using (5.18) and (5.15), we have the corresponding AIC0 and AIC1 :
N 2tˆl |ˆz (l) − zˆ0 |2 1 AIC0 − AIC1 = ln + ˆ 2 2tˆl (1 − e−2t l = ) l=1 |ˆz (l) − zˆ0 e−(i!ˆ 0 +1= )tˆl |2 − 1 − e−2tˆl =
;
(5.19)
where we used the scaled variables and the parameters: zˆ = z=; ! ˆ 0 = !0 =; = =, and tˆ=t ·. In Fig. 6, we show Pd as a function of and in Fig. 7 as a function of ! ˆ 0, where it is found that Pd decreases with increasing value of . When the 0uctuation speed becomes slow, the quantum mechanical damping process of the oscillator looks like the diJusion process. We should also note that the K–L information I (or the Boltzmann relative entropy S), (2.8), is de1ned in the framework of a probability theory. Among several candidates
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
113
Fig. 7. Pd as a function of !ˆ 0 for = 0:5 and 1.0.
in generalizing the AIC method to include the quantum mechanical systems, we used the Q-function. The corresponding entropy is sometimes referred to as the Wehl–Lieb entropy [17,18]. This is a well-behaved quantity and gives no di4culty while the other functions like the Wigner function and the P-function show ill behavior in calculating the entropy as mentioned earlier. This is the reason why we adopted the Q-function in generalizing AIC to the quantum systems.
6. Concluding remarks In the previous paper [3], we studied the detection of signals in a Gaussian white noise on the basis of the AIC method. In this paper, we proceeded to treat the discrimination problems of the colored stochastic processes. That is, the Gaussian–Markovian (Ornstein–Uhlenbeck) process is examined in comparison with the diJusion process. The discrimination probability Pd depends on the parameter and the observation time T . The longer the observation time T and the smaller the strength of 0uctuations , we found the higher values of Pd . This result can be intuitively expected. Moreover, we found a simple linear relation between the two requiring Pd greater than a speci1ed value (Fig. 2). The method was further applied to discriminate dynamically time-evolving states of physical systems. As a typical system with dynamical state change, we treated a problem of laser action where the nonlinearity and the 0uctuations play essential roles. The dynamical state change was clearly speci1ed when we calculate AICs. It should be emphasized that this kind of dynamical state change cannot be recognized by a simple observation of the data. Moreover, we showed that the methodology can be extended even for quantum mechanical systems when we employ the Q-function in the de1nition of the K–L entropy. Among several possibilities
114
M. Shingu-Yano, F. Shibata / Physica A 293 (2001) 100 –114
in de1ning quantum mechanical entropy, use of the Q-function (Wehl–Lieb entropy) is reasonable because of its nonsingular behavior in contrast with other functions. In fact, we could apply the method to the quantum mechanical relaxation process of a damped oscillator. The problems studied in this paper, especially those of the dynamical state change and the quantum mechanical system, have not yet been treated by this method. In conclusion, we showed usefulness and applicability of the AIC method beyond a scope of statistics to include the discrimination problems in dynamically time-evolving physical systems. References [1] R.N. McDonough, A.D. Whalen, Detection of Signals in Noise, 2nd Edition, Academic Press, Sandiego, New York, Tokyo, 1995. [2] H. Akaike, in: B.N. Petrov, F. Csaki (Eds.), Proceedings of the Second International Symposium on Information Theory, Academiai Kiado, Budapest, 1973, p. 267. [3] H. Akaike, in: E. Parzen, K. Tanabe, G. Kitagawa (Eds.), Selected Papers of Hirotsugu Akaike, Springer, New York, 1998. [4] M. Shingu-Yano, A. Kobayashi, Y. Nakamura, F. Shibata, J. Grad. School Humanities Sci., Ochanomizu University, Tokyo 2 (2000) 205. [5] R. Kubo, M. Toda, N. Hashitsume, Statistical Physics II, Springer, Berlin, Heidelberg, New York, Tokyo, 1985. [6] C. Uchiyama, F. Shibata, Physica A 198 (1993) 538. [7] F. Shibata, Y. Shimoo, Physica A 215 (1995) 87. [8] H. Haken, in: L. Genzel (Ed.), Handbuch der Physik, Vol. XXV/2c, Light and Matter, Springer, Berlin, Heidelberg, New York, 1970. [9] W.H. Louisell, Quantum Statistical Properties of Radiation, Wiley, New York, London, Sydney, Tront, 1973. [10] M. Sargent III, M.O. Scully, W.E. Lamb Jr., Laser Physics, Addison-Wesley, Reading, MA, 1974. [11] H. Risken, The Fokker–Planck Equation, 2nd Edition, Springer, Berlin, Heidelberg, New York, London, Paris, Tokyo, 1989. [12] F. Haake, Springer Tract in Modern Physics, Vol. 66, Springer, Berlin, Heidelberg, New York, 1973, p. 98. [13] R.J. Glauber, Phys. Rev. 131 (1963) 2766. [14] E. Wigner, Phys. Rev. 40 (1932) 749. [15] R. Kubo, J. Phys. Soc. Jpn. 19 (1964) 2127. [16] K. Husimi, J. Phys. Math. Soc. Jpn. 22 (1940) 264. [17] A. Wehrl, Rept. Math. Phys. 16 (1979) 353. [18] E.H. Lieb, Commun. Math. Phys. 62 (1978) 35.