Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://lnfm1.sai.msu.ru/~marat/articles/eng/Err_Estimated_Par_2_eng.pdf Äàòà èçìåíåíèÿ: Tue Jul 21 12:07:16 2009 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 23:29:34 2012 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: m 63

ISSN 1063-7729, Astronomy Reports, 2009, Vol. 53, No. 8, pp. 722­749. c Pleiades Publishing, Ltd., 2009. Original Russian Text c M.K. Abubekerov, N.Yu. Gostev, A.M. Cherepashchuk, 2009, published in Astronomicheski Zhurnal, 2009, Vol. 86, No. 8, pp. 778­806. i

Estimation of Parameter Errors in Inverse Problems. Determining Limb-Darkening Coefficients in Classical Eclipsing Systems
M. K. Abubekerov, N. Yu. Gostev, and A. M. Cherepashchuk
Sternberg Astronomical Institute, Moscow State University, Moscow, Russia
Received December 15, 2008; in final form, December 20, 2008

Abstract--The estimation of parameters and their errors is considered using of the eclipsing binary system YZ Cas as an example. The error intervals differential-correction and confidence-region methods. The error intervals and are investigated, and the reliability of limb-darkening coefficients derived from analyzed. A new method for calculating parameter errors is proposed. PACS numbers: 97.80.Hn, 95.75.Pq DOI: 10.1134/S1063772909080046

the observed light curve are calculated using the reliability of the methods the observed light curve

1. INTRODUCTION A comparative analysis of various methods for estimating parameter errors for one-, two-, three-, and four-parameter functions and the observed light curve of the eclipsing binary system YZ Cas was carried out in [1] using numerical simulations. The numerical simulations yielded estimates of the unknown parameters of the problem and their error intervals, which were calculated using the differential-correction, Monte Carlo, and confidence-region methods. The error intervals and reliability of the various methods were studied. The current study is a continuation of the analysis begun in [1]. In contrast to the current study, in [1], the trustworthiness of the methods used to estimate the error intervals were evaluated using an incomplete set of parameters. The analysis of the observed light curve for the binary system YZ Cas was based on three parameters: the radius of the primary and secondary components r1 and r2 and the orbital inclination of the binary i. In addition, the residual used depended on these three parameters and was minimized over the remaining, linear parameters, which included the limb-darkening coefficients for the primary and secondary stars x1 and x2 and the ratio of the brightnesses at the centers of the stellar disks, (1) (2) I0 /I0 . It is important that minimizing over linear parameters does not change the statistics of the problem, only decreases the number of degrees of freedom. This enabled us to analyze the errors in the non-linear parameters r1 , r2 , and i in a mathematically rigorous way. The error intervals were found as projections onto the r1 , r2 , and i parameter axes of the threedimensional confidence region containing the values

of r1 , r2 ,and i for which the residual was less than the quantile corresponding to a specified confidence level. In the current study, in addition to the parameters listed above, the limb-darkening coefficients for the primary and secondary stars are unknown. We obtain fits to the observed light curve of YZ Cas using both linear and quadratic limb-darkening laws. We focus on the reliability of the limb-darkening coefficients and their error intervals. This problem is especially relevant in connection with the appearance of high-accuracy photometric observations from satellites such as COROT and KEPLER, which will soon lead to the availability of massive databases of precise photometric observations. Therefore, it is important to know the reliability of estimates of photometric elements of both classical eclipsing binary systems and eclipsing systems with exoplanets at specified levels of accuracy. As was emphasized in [1], in both the classical least-squares method (or the differential-correction method for non-linear problems) and the Monte Carlo method, the errors of unknown parameters are estimated in the framework of a normal distribution, which are obeyed by the central values of parameters, and searches for the true values of the parameters are carried out using statistics that arise from a normal distribution (such as 2 statistics, where M is the M number of points in the observed light curve). However, it is most reasonable to calculate the errors of parameters using the same statistics that were used to search for the central values of the parameters themselves (for example, the statistics for a residual distributed as 2 ). Although the resulting parameter M errors will be larger than the errors found using the
722


PARAMETER ERRORS IN INVERSE PROBLEMS

723

least-squares or Monte Carlo method, based on the use of confidence regions most trustworthy fits to observational avoids the need to separately consider " "external" errors in parameters.

an approach provides the data; it also internal" and

Here, is the polar distance from the center of the stellar disk and r the radius of the stellar disk. The brightness was assumed to be zero at any point in the disk of the planet. We denote I0 and I0 to be the brightnesses at the disk centers, x1 and x2 the limb-darkening coefficients, and y1 and y2 the quadratic limb-darkening coefficients of the primary and secondary stars. The unknown parameters of the model for the two (1) (2) stars are r1 , r2 , i, I0 , I0 , x1 , x2 , and, in the case of the non-linear limb-darkening law, y1 and y2 . The unit length in our models is the distance between the stellar centers (or the centers of the star and planet), a = 1. There is no "third light" in the model. The light curve of the binary system is determined by the following three equations. 1. The total luminosity of the components, which describes the out-of-eclipse brightness:
r
1

(1)

(2)

2. MODEL FOR A CLASSICAL BINARY SYSTEM AND A BINARY SYSTEM WITH AN EXOPLANET In view of the methodological character of our work, we consider a simple model with two spherical stars with thin atmospheres in circular orbits without any effects due to the proximity of the components. This model can easily be realized on modern computers, and makes it possible to obtain a large number of solutions of the problem with a comparatively small amount of computer time. This type of sphericalstar model for a binary system is physically based on the case when the degree of Roche-lobe filling is small, µ < 0.5. This occurs in the YZ Cas system, and also in most cases of observed eclipsing stars with exoplanets. Our model considers the motion of the stellar disks projected onto the plane of the sky, i.e., perpendicular to the line of sight. Figure 1 shows the geometry of the stellar disks during eclipse. Here, r1 and r2 are the radii of the primary and secondary star (or of the star and planet), is the distance between the centers of the stellar disks, and and are polar coordinates of an arbitrary point on the surface of the disk of the primary (the coordinate origin is located at the geometrical center of the disk). The distance between the centers of the stellar disks is given by the expression 2 = cos2 i +sin2 i sin2 , (1) where i is the orbital inclination of the binary system and is the current value of the orbital phase angle (the orbit of the system is circular with radius a = 1). For the brightness distribution over the disk of each star, we used a linear limb-darkening law, I () = I
0

r

2

2
0

I

(1)

()d +2
0

I

(2)

( )d = L

full

.

(4)

2. The loss of brightness in the first (secondary, in the case of YZCas) minimum, due to the eclipse of the smaller star by the larger companion: L
full

-L

(1)

( ) =
S ()

I

(2)

( )dS,

(5)

where S () is the area of overlap of the two disks. 3. The loss of brightness in the second (primary) minimum, due to the eclipse of the larger star by the smaller star:

2 1-x+x 1- 2 r

L (2)

full

-L

(2)

( ) =
S ()

I

(1)

()dS.

(6)

and also a quadratic limb-darkening law, which differs from the linear law in the presence of a non-zero quadratic limb-darkening coefficient y : I () = I0 1 - x 1 - -y 1 - 2 1- 2 r 1-
2

Equations (1), (4), (5), and (6) fully describe the observed light curve, and contain a set of parameters (1) that depends on the model considered: r1 , r2 , i, I0 , I0 , x1 , x2 , y1 , y2 . Substituting the brightness distribution function approximated by the corresponding limb-darkening law (2) or (3) into the integral and carrying out the integration yields a system of nonlinear, algebraic equations in the corresponding parameters.
(2)



2 r2

(3)

.

ASTRONOMY REPORTS

Vol. 53 No. 8 2009


724

ABUBEKEROV et al.

r1 r2



d

Fig. 1. Model for two eclipsing spherical stars, projected onto the plane of the sky.

3. DERIVATION OF THE COMPUTATIONAL FORMULAS In the linear limb-darkening law, the brightness at a point in the stellar disk is I () = I0 1-x+x 1- 2 r2 ,

limb-darkening coefficients, and r the radius of the star. Introducing the new parameters X0 = I0 (1 - x - 2y ) (9) X1 = I0 (x +2y ) , X2 = I0 y we obtain 2 2 + X2 2 r2 r

where is the distance of the point from the disk center, I0 the brightness at the disk center, x the limb-darkening coefficient, and r the radius of the star. Introducing the new parameters X0 = I0 (1 - x) , (7) X =I x
1 0

I () =

X0 + X1

1-

.

we find that I () = X0 + X1 1- 2 r2 . (8)

In these variables, the brightness at a point on the stellar disk for the non-linear limb-darkening law differs from the corresponding brightness for the linear law in the presence of an additional term containing X2 . This enables us to obtain a linear model by setting X2 = 0 in the non-linear model. At orbital phase = 0, component 1 eclipses component 2. The total brightness of the star in the non-linear model is
r

L

(s)

= 2
0

I

(s)

2 (s) 2 (s) 2 ()d = X0 rs + X1 rs 3 s = 1, 2.

In the non-linear limb-darkening law, the brightness at a point in the stellar disk is I () = I
0

1-x 1- 2 1- 2 r

1-
2



r2

2

1 (s) 2 + X2 rs , 2

The total brightness of the system of two stars in the non-linear model outside of eclipse is
2 Lfull = L(1) + L(2) = X0 r1 2 (1) 2 1 (1) 2 (2) 2 + X1 r1 + X2 r1 + X0 r2 3 2 2 (2) 2 1 (2) 2 + X1 r2 + X2 r2 . 3 2 (1)

-y 1 -

,

(10)

where is the distance of the point from the disk center, I0 the brightness at the disk center, x, y the

ASTRONOMY REPORTS Vol. 53 No. 8

2009


PARAMETER ERRORS IN INVERSE PROBLEMS

725
dec 0

To obtain a model for a system consisting of a star and planet, it is sufficient to set L(2) = 0 (the planet does not radiate). To universalize the external form of the computational formulas and the light-curve minima, and to decrease the number of equations needed, the eclipsing component (near component relative to a terrestrial observer) is ascribed the subscript n and the eclipsed component (farther component relative to a terrestrial observer) the subscript f . When directly computing the light-curve minima at orbital phases -/2 < < /2 (or cos > 0), we must replace the variable r1 with rn ,and thevariable r2 with rf .Atorbital phases cos < 0, we must carry out the opposite substitution, replacing the variable r1 with rf ,and the variable r2 with rn . In the new notation, the drop in brightness during eclipse is L
dec

and the expressions for L L
dec 0

and L

dec 1

obtained in [1]: (14)

2 (,rf ,rn ) = (,rf ,rn )rf 1 2 +(,rn ,rf )rn - Q (,rf ,rn ) , 2

L
2 rf

dec 1

(,rf ,rn )

(15)

=
0

(, , rn ) 1 - 2 d. rf
dec 0

Analogous to the way in which L obtained in [1], we obtain for Ldec : 2 rf L
dec 2

and L

dec 1

were

(,rf ,rn ) = 2
0

3 (,,rn )d

(16)

(

(f ) ,rf ,rn ,X0

(f ) ,X1

(f ) ,X2

)

(11)

=
S ()

I

(f )

()dS,

4 rf r2 2 = (,rf ,rn ) + n 22 + rn (,rn ,rf ) 2 2 1 2 2 2 +5rn + rf Q (,rf ,rn ) . - 8

where is the distance between the disk centers and S () the area of overlap of the disks. To compute the integral (11), as was done in [1], we introduce the functions , x < -1, Ax arccos x, (12) -1 x 1, 0, x > 1 and Qx x, x 0, 0, x < 0, x2 +2 - y 2x
2

The partial derivatives of L L
dec 2

dec 2

have the form (17)

(,rf ,rn ) 2 = 2rn (,rn ,rf ) 2 2 2 + rn + rf Q(,rf ,rn ), - 2

L

dec 2

(,rf ,rn ) 3 = 2rf (,rf ,rn )rf , rf
dec 2

(18)



L

(,rf ,rn ) 2 = 2rn 2 + rn rn â (,rn ,rf ) - 2rn Q(,rf ,rn ).

(19)

(,x,y ) A

,

For a circular orbit, the distance between the disk centers depends on the phase and orbital inclination i as (, i) = cos2 i +sin2 i sin2 . (20) , The brightness of the binary system depends on the phase as follows: (21) L(, i, r1 ,r2 ,X0 ,X1 ,X0 ,X1 ) (1) (1) (1) Ldec ((, i),r1 ,r2 ,X0 ,X1 ,X2 ), cos < 0, - (2) (2) (2) Ldec ((, i),r2 ,r1 ,X0 ,X1 ,X2 ), cos > 0.
(1) 0 (1) (1) (2) (2)

Q (,rf ,rn ) Q
2 rf

- ( - rn )

2

( + rn ) -

2

2 rf

together with a polar-coordinate system with its origin at the center of the disk of the eclipsed star and the polar angle measured from the center of the disk of the eclipsed component "f" to the center of the disk of the eclipsing component "n" (Fig. 1). We then have L
(f ) dec

=L

full

(,rf ,rn ,X0 ,X1 ,X2 ) (,rf ,rn )+ X1 L
(f ) X2 (f ) dec 1

(f )

(f )

(f )

(13)

= X0 L

dec 0

(,rf ,rn )

+

L

dec 2

Introducing the functions L (, r1 ,r2 ,i) (22)

(,rf ,rn ),
Vol. 53 No. 8 2009

ASTRONOMY REPORTS


726
2 = r1 -

ABUBEKEROV et al.

Ldec ((, i),r1 ,r2 ), cos < 0, 0 0, cos > 0, L
(1) 1

(,rn ,rf ) 4rn 2 2 rf - (rn - )2 - 2 2 3rf (2 - rn ) 9

(, r1 ,r2 ,i) cos < 0, cos > 0, + âF + (, r1 ,r2 ,i)

=

22 r - 31

Ldec ((, i),r1 ,r2 ), 1 0, L
(1) 2

2 2 [rf - (rn - )2 ][(rn +)2 - rf ] rf

(, r1 ,r2 ,i) cos < 0, cos > 0,

12 = r1 - 2

Ldec ((, i),r1 ,r2 ), 2 0, L
(2) 0

(,rn ,rf ) 4rn 2 2 rf - (rn - )2

2 2 2 2 Q rf - (rn - )2 (7rn +2 - 4rf ) 9rf âE (,rn ,rf ) 4rn 2 2 rf - (rn - )2 .

2 = r2 -

Ldec ((, i),r2 ,r1 ), cos > 0, 0 0, cos < 0, L
(2) 1

(, r1 ,r2 ,i) cos > 0, cos < 0,

22 = r2 - 3

Ldec ((, i),r2 ,r1 ), 1 0, L
(2) 2

4. GENERAL RELATIONS FOR THE LEAST-SQUARES METHOD Since our goal is to elucidate the relationship between various methods for estimating parameter errors, let us recall some relationships from mathematical statistics that are relevant to our calculations. Let us consider a linear model specified by the sequential set of functions g0 ( ) ... gP ( ) and the linear functions expressed in terms of them1 f lin (, 1 ... P ), defined for real 1 ... P and from the set {1 ... M }:
P

(, r1 ,r2 ,i) cos > 0, cos < 0,

=

12 r - 22

Ldec ((, i),r2 ,r1 ), 2 0,

we can write the total brightness as the linear combination L(, r1 ,r2 ,i,X0 ,X1 ,X0 ,X1 ) = + +
(1) (1) X0 L0 (, r1 ,r2 ,i)+ (1) (1) X2 L2 (, r1 ,r2 ,i)+ (2) (2) X1 L1 (, r1 ,r2 ,i)+ (1) (1) (2) (2)

(23)

(1) (1) X1 L1 (, r1 ,r2 ,i) (2) (2) X0 L0 (, r1 ,r2 ,i) (2) (2) X2 L2 (, r1 ,r2 ,i).

f

lin

(, 1 ... P ) = g0 ( )+
p=1

gp ( )p .

(25)

The partial derivatives of a function of the total brightness are expressed in terms of the partial derivatives of the functions Ldec , Ldec , and Ldec with 0 1 2 respect to the variables , rf , and rn and the partial derivatives of (, i) with respect to and i. We will show that the function Ldec can be ex1 pressed in terms of the elliptical functions , E ,and F and the function , for whose calculation there exist efficient algorithms [2]: L
dec 1

(,r1 ,r2 ) = +Q
2 rf

2 2 (rn - )rf 3

(24)

1 - (rn - )2

Here, {1 ... M } corresponds to the set of M points at which observations are conducted. Let us specify the vector (1 ... P )T , corresponding to the set of true values of the physical quantities (here, M is the total number of observational points), and the vector (w1 ... wM )T representing the weighting coefficients, assuming that the matrix Aqp = M m=1 gq (m )gp (m )wm is non-degenerate. We also specify a vector representing the random observed quantities = (1 ... M )T , assumed to be statistically independent, i.e., cov (i ,j ) = 0 for i = j , and to be normally distributed with the mathematical expectation values M(k ) = f lin (k , 1 ... P ). Moreover, we assume that the measurements of the random quantities 1 w1 ... M wM are equally accurate; i.e., 2 (1 w1 ) = 2 (2 w2 ) = ... = 2 (M wM ) = 2 , 0 where 2 (·) denotes the operation of calculating the
1

â

3 2rf ( + rn ) 4rn - ; 3( - rn ) (rn - )2

The form f lin [(25)] in which the first term is independent of the linear parameter is chosen here for convenience in practical use. 2009

ASTRONOMY REPORTS Vol. 53 No. 8


PARAMETER ERRORS IN INVERSE PROBLEMS

727

dispersion and 2 is called the unit-weight dispersion. 0 2 2 Denoting 2 (m ) = m , we find that wm = 2 /m . 0 We specify the residual functional for a given model by the expression R
M lin

the statistical independence of 1 ... M , and the fact that cov (i ,i ) = 2 (i ) = 2 /wi , we find the covari0 ance matrix of the parameters c ... c : 1 P cov(c ( ),c ( )) p q
M P 2 ( 2 (m )wm m=1 P i,j =1 M

(30)

(1 ... P ,1 ... M )
lin 2

(26) = Ainv Ainv gi (m )gj (m )) ip jq

=
m=1 M

m - f

(m ,1 ... P ) wm
P p=1



2

=
m=1

m - g0 (m ) -

gp (m )p wm .

=

2 0

(Ainv ip i,j =1

Ainv jq
m=1

wm gi (m )gj (m ) = 2 Ainv ). 0 pq

The values of the parameters c ( ) ... c ( ) yielding 1 P the minimum of the residual functional (26) for fixed 1 ... M (which we will call their central values), are obtained from the solution of the system of P linear equations R
P lin

The dispersions of the central values of c ( ) ... c ( ) are the diagonal elements of the 1 P covariance matrix: 2 (c ( )) = 2 Ainv . p 0 pp (31)

(1 ... P ) = 2Bq q q = 1 ... P,

(27)

In the model obtained via the linear substitution of i parameters using the non-degenerate matrix Cp ,
P

-2
p=1

Aqp p = 0,

p = C0 +
i=1

i Cp i ,

where
M

the dispersions of the new parameters 1 ...P are found from the formula gq (m )gp (m )wm , (28)
P

Aqp =
m=1 M

2 ( p ) =
i,j =1

c

ij Cp Cp cov (c ,c ). i j

(32)

Bq =

[m - g0 (m )]gq (m )wm .
m=1

These values are equal to
P

c p
M

( ) =
q =1

Ainv Bq qp
P

(29)

Knowing the dispersion of the central value of a parameter, we can construct the interval within which the true value of the parameter falls with a specified probability. For this, it is sufficient to note that it follows from the normal distribution for the central value of the parameter that ¯ P c ( ) - p ( ) (c ( )) = , p p (33)

=
m=1

[m - g0 (m )]wm
q =1

Ainv gq (m ), qp

p = 1 ... P, where Ainv are elements of the inverse matrix to A: qp Ainv (A-1 )qp . qp exThus, the central values of ( ) ... and pressed as linear combinations of the are consequently normally distributed. Their mathematical expectation values are 1 ... P [which is obtained if we take the mathematical expectations of both parts of (29) and substitute them into (25)]. Using the linearity of the covariance operation cov (·, ·),
ASTRONOMY REPORTS Vol. 53 No. 8 2009

where P denotes the probability that a condition is satisfied; depends on a chosen probability (confidence level) , and is given by the root of the equation


2
0

exp -

t2 2

dt = .

c 1

c ( ) are P 1 ... M ,

For example, the confidence level = 0.6827, 0.9545, and 0.9973 for = 1, 2, and 3, respectively (corresponding to 1, 2, or 3 ). When 0 is not known (for example, as in the case of a real observed light curve), in place of (31),


728

ABUBEKEROV et al.

we must find the errors using formulas obtained by replacing 0 in (31) and (30) with the quantity wq [4]:
2 0 ( ) =

R

lin

(c ( ) ... c ( ),1 ... M ) 1 P , M -P

is non-degenerate for (1 ... P )T B . We will also denote the elements of the inverse matrix to (37) as Ainv (1 ... P ). qp The residual functional is specified by the expression R(1 ... P ,1 ... M )
M

(34)

called the rms estimated unit-weight dispersion. The resulting rms estimates for the dispersions (which are random quantities) are
2 est 2 (c ( )) 0 ( )Ainv . p pp

(38)

(35)

=
m=1

[m - f (, 1 ... P )]2 wm ,

Here, 0 is a random quantity, and (c ( ) - p p )/est (c ( )) obeys a Student distribution with ¯ p M - P degrees of freedom. However, this is already close to a normal distribution for sufficiently large ¯ M -P 10, and we can assume P(|c ( ) - p | p c ( ))) P(|c ( ) - p | est (c ( ))); i.e., ¯ (p p p the probability that the true value falls in an interval constructed by multiplying the rms estimate of the dispersion2 by the corresponding coefficient ( ) will be fairly close to . Let us now consider a model specified by an arbitrary, in general non-linear, function f (, 1 ... P ), defined for {1 ... M } and for the vectors (1 ... P )T B , where B is some region of real, Euclidean space. We assume that f (, 1 ... P ) is differentiable with respect to 1 ... P over this entire region. As for the linear model, we specify the 1 ... P , w1 ... wM , the normally distributed random quantities 1 ... M , the unit-weight dispersion 0 , and the residual functional. The random quantities 1 ... M display a normal distribution, and M(k ) = f (k , 1 ... P ),
2 2 2

and we assume that, for fixed 1 ... M , it is convex in the variables 1 ... P and reaches its minimum in the region of B . In the differential-correction method, the function f is replaced by a Taylor series expansion up to a linear term at the minimum of the residual functional. Dispersions found from a least-squares solution for the corresponding linear model are used as estimates of the dispersions of the optimal values 1 ... P .
c c We denote 1 ( ) ... P ( ) to be the parameter values (which we will call central values) providing the minimum of the residual function R(1 ... P ,1 ... M ) for specified 1 ... M , and set in (25), and then (28) and (29), c p = p - p ( ),

(39)

g0 ( ) = f (, ( ) ... ( )), f c c (, 1 ( ) ... P ( )), gp ( ) = c p where p = 1 ... P . Then, the quantities
c c c c covo (q ( ),p ( )) 2 Ainv (1 ( ) ... P ( )) (30 ) 0 qp

c 1

c P

(36)

(1 w1 ) = (2 w2 ) = ... = (M wM ) = 2 , 0 where M(k ) denotes the mathematical expectation of k and 2 (·) the operation of finding the dispersion. We also assume that the matrix Aqp (1 ... P )
M

and
2c c c o (p ( )) 2 Ainv (1 ( ) ... P ( )), 0 pp

(31 )

(37)

=
m=1
2

f f (m ,1 ... P ) (m ,1 ... P )wm q p

obtained by substituting (39) into (30) and (31) for the covariances and dispersions of the central values in the linear model are taken as estimates of the c c covariances and dispersions 1 ( ) ... P ( ). In the case of a real observed light curve, when the unitweight dispersions are unknown, analogous to the linear case, we use in place of 0 the rms estimate of the unit-weight dispersion
2 0 ( ) = c c R(1 ( ) ... P ( ),1 ... M ) M -P

Note that variation of the rms estimates of the dispersions (which are random quantities) from sample to sample, i.e., their scatter for various measurements of the same curve, can be appreciably larger than the difference between the Student and Gaussian distributions. This scatter is inversely proportional to the square root of the number of points in the light curve.

(34 )

and the corresponding approximate rms estimates of the parameter dispersions
2 est c 2 c c (p ( )) 0 ( )Ainv (1 ( ) ... P ( )). pp

(35 )

ASTRONOMY REPORTS Vol. 53 No. 8

2009


PARAMETER ERRORS IN INVERSE PROBLEMS

729

Using this approximation in the computations assumes that we can neglect variations in the derivatives of the function f in (37) calculated with the central parameter values when is varied in the vicinity of its mathematical expectations. It is obvious that, if the function f (, 1 ... P ) is linear in 1 ... P , then (31 )­(35 ) coincide with (31)­(35); i.e., in the case of a linear function, the differential-correction method is identical to the least-squares method. Note that the described methods make the assumption that the model used is fully correct, and statistics for a normal distribution are used to estimate the parameter errors. In the model obtained via the reversible parameter substitution (40) p = p (1 ... P ), where are smooth functions and the matrices p (1 ... P ) are non-degenerate, we can obi tain an expression for approximate estimates of the central-value dispersions for the new parameters c ( ) ... c ( ) by replacing the right-hand side 1 P of (40) with a Taylor expansion up to the linear c c term at the point (1 ( ) ... P ( ))T , setting in forc c p (1 ( ) ... P ( )) i and replacing mula (32) Cp = c i there the covariances c ( ) ... c ( ) with the corre1 P sponding approximate estimates of the covariances c c (1 ( ) ... P ( )):
2 o ( p ( )) P c

(41)

=
i,j =1

c c c c p (1 ( ) ... P ( )) p (1 ( ) ... P ( )) c c i j c c â covo (i ( ),j ( )).

We should make an important comment about the described least-squares method (and, in the case of a non-linear model, the differential-correction method). In these methods, the statistical hypothesis H of the adequacy of the model to describe the observational data is usually not verified, and instead it is immediately assumed that the model used is correct. In other words, when adopting the hypothesis H in this case, the possibility that the researcher could make an error of the second kind (the model is not correct, but is adopted in accordance with some statistical criterion) is excluded. Precisely because of this, and also due to the use of a simple statistical normal distribution when estimating the uncertainties in the unknown parameters [see (31), (35), (31 ), (35 )], it is possible to obtain relatively small parameter errors in the least-squares and differential-correction methods (and also in the Monte Carlo method).
ASTRONOMY REPORTS Vol. 53 No. 8 2009

In reality, it is not known a priori how appropriate the assumption of the correctness of the model is, since, as a rule, the number of observed points in the light curve is much larger than the number of unknown parameters, so that the search for the parameters of the problem is strongly overdetermined. Therefore, the search for the parameters and their uncertainties should be carried out jointly with a statistical verification of the hypothesis H of the adequacy of the model. Further, as was noted in [5, 8, 12], the acceptance or rejection of a hypothesis H according to some statistical criterion does not provide a final, logical basis for its ultimate acceptance or rejection (because of the statistical nature of the criterion). Four cases are possible. 1. The hypothesis H is correct and is accepted by the criterion. 2. The hypothesis H is incorrect and is rejected by the criterion. 3. The hypothesis H is correct, but is rejected by the criterion (an error of the first kind). 4. The hypothesis H is incorrect, but is accepted by the criterion (an error of the second kind). There exists for a selected statistical criterion a non-negative probability to admit an error of the first kind, denoted 0 0. The quantity 0 is called the significance level of the criterion. It is obvious that, in the case of an adequate model, the probability of case (1) is = 1 - 0 . For example, we can use the following criterion to verify a hypothesis H (for more detail, see [8, 12]). Let us choose some random quantity (statistic) (u) that ~ depends on some experimental data u, with the dis~ tribution (u) known (for example, this is a 2 dis~ M tribution, where M is the number of observed points in a light curve). We specify a priori some significance level 0 and calculate the number M (quantile) such that P {(u) > 0 } = 0 for all parameter values of ~ the model if the hypothesis H is correct. If we obtain that (u) > M for any values of parameters ~ for a specific realization of the experimental data (a specific light curve) u, the hypothesis H is rejected. ~ On the contrary, if ((u) M ), the hypothesis H is ~ accepted. At the same time, we must not forget that the model is accepted, not because it is necessarily correct, but because there is no basis for it to be rejected. In the best cases, the criterion that the minimum reduced chi-squared residual be close to unity is used when applying the least-squares, differentialcorrection, or Monte Carlo method to find parameters 2 1, where P is the and their uncertainties: M -P M -P number of unknown parameters. However, here, no quantitative measure of the adequacy of the model to


730

ABUBEKEROV et al.

R(, ) M() ­ R(, )
­ R(, ) ­ Rmin()

2M (, ) 21() 1()

Rmin() c()

Fig. 2. A one-parameter problem. The dependence of the residual R(, ) on the parameter for a fixed vector comprising the random quantities and the main statistical quantities. Rmin ( ) is the minimum residual in terms of , R(, ); P ( ) (P = 1) the quantile of the 2 (P = 1) distribution; M ( ) the quantile of the P 2 distribution; 21 ( ) and 2M (, ) the widths of the M confidence intervals obtained from methods using the statistics R(, ) - Rmin ( ), R(, ) distributed as 2 P (P = 1) and 2 , respectively (see text for detail). M

rigorous analytical relation between the error intervals for the above methods for a linear model. Namely, we derive the function ap (t) for the distribution of the ratio of the confidence interval for the 2 method to M its value for the 2 method. In the one-dimensional P problem (P = 1), the 2 confidence interval coinP cides with the corresponding interval for the leastsquares method, while, for arbitrary P , these intervals are related by a well defined dependence (which will be presented below). We also study the behavior of ap (t) as a function of the number of observational points M and the selected confidence level . Let us first consider the one-dimensional case. We consider the residual in a one-parameter linear model that depends quadratically on the one parameter . We also assume for convenience that 0 = 1 (the statements below will remain valid for any 0 ). This residual can be written in the form (as follows from (26), where the terms that are quadratic in do not contain the random quantities ) R
lin lin (, ) = C [ - c ( )]2 + Rmin ( ),

(42)

the observational data is considered, namely the significance level 0 at which the model can be rejected. We have investigated this question here (see below). 5. ANALYTICAL RELATION BETWEEN THE ERROR INTERVALS OBTAINED USING STATISTICS DISTRIBUTED AS 2 AND 2 , M P AND THE DIFFERENTIAL-CORRECTION METHOD In our previous study [1], we established via numerical simulations that, on average, the error interval obtained using the confidence-region method based on a statistic distributed as 2 , where M is M the number of measurement points (we will call this the 2 method), exceeds the corresponding interval M obtained using the differential-correction method by a factor of three to five. The error interval obtained using the 2 method also exceeds the corresponding M interval obtained using the confidence-region method based on a statistic distributed as 2 , where P is the P number of parameters being searched for (we will call this the 2 method). The sizes of the confidence inP tervals obtained using the least-squares method with exact values of the dispersions and the 2 method P do not depend on the sample of random quantities, in contrast to the 2 confidence intervals, which are M themselves random quantities. We present below a

lin where Rmin is its minimum value in terms of , reached for c ( ),and the coefficient C is expressed in awell defined way in terms of the specified parameters of the model (w and g), and does not depend on the sample of random quantities . The quantile n ( ) for (0, 1) is given by

2 (n ( )) = ; n i.e., as an inverse function to 2 . n Let us briefly recall the residual distribution. We showed in [1] that the residuals are distributed as R
lin lin (, ) - Rmin ( ) 2 , 1

(43) (44) (45)

R

lin

(, ) 2 , M
2 M -1

lin Rmin ( )

,

where is the true value of the parameter , and the symbol denotes "distributed as." Figure 2 depicts the dependence of the residual R(, ) on the parameter for a fixed sample of random quantities and its dependence on the main statistical parameters. The parabola representing this dependence in a linear model can shift in the upper half-plane for different , while maintaining its shape, since the coefficient of the quadratic term in does not depend on . The error interval obtained using the 2 method for 1 aspecified confidence level is defined as the distance between the roots of the quadratic equation in :
lin C [ - c ( )]2 + Rmin ( )

(46)

=

lin Rmin

( )+ 1 ( ).
2009

ASTRONOMY REPORTS Vol. 53 No. 8


PARAMETER ERRORS IN INVERSE PROBLEMS

731

We denote half this interval as 1 ( ). It is obvious that 1 ( ) obeys the relation
2 C1 ( ) = 1 ( )

= It follows M (, ) > 0 zero for a 2 -1 (M ( M

2 M -1

(M ( ) - 1 ( )t) .

(47)

and is independent of the sample of random quantities . Wethushave 1 ( ) . (48) C= 2 1 ( ) The error interval obtained using the method for a specified confidence level is the distance between the roots of the quadratic equation in :
lin C [ - c ( )]2 + Rmin ( ) = M ( ). 2 M

from (45) that the probability that (i.e., that the confidence interval is nonspecified confidence level ), is )).

The corresponding conditional probability is P = M (, ) > 0|
2 M -1

M (, ) >t 1 ( )
2

(54)

(49)

We set M (, ) equal to half this interval when this last equation has positive roots (i.e., the corresponding confidence set is not empty) and equal to zero otherwise. For positive values of M (, ),
2 CM

M ( ) - 1 ( )t 2 M -1 (M ( ))

or P M (, ) > 0|
2 M -1

(, )+

lin Rmin

( ) = M ( ).

(50)

M (, ) t 1 ( )

(55)

The quantity M (, ) depends on the random quantities and is itself a random quantity. Using (45) we obtain 2 (51) M ( ) - CM (, ) 2 -1 . M This last expression indicates that
2 P((M ( ) - CM (, )) < t) = 2 M -1

(t).

(52)

We obtain from this last expression via straightforward linear and shift transformations in t and using (48) P
2 M (, ) >t 2 1 ( )

(53)

M ( ) - 1 ( )t2 . 2 -1 (M ( )) M The right-hand side of this last expression is a function of the distribution of positive values of M (, ) , which is the ratio of the confidence interval 1 ( ) for the 2 method to the corresponding interval M obtained when the main statistic is distributed as 2 , 1 for a specified confidence level . The corresponding density of the distribution a1 (t), obtained by differentiating the distribution function (with accuracy to within a normalization coefficient A1 ), is equal to =1-

e 1 (1 ( )t2 -M ( )) t ( ) - ( )t 2 M 1 0 < t < M ( ) , a1 (t) = A1 1 ( ) M ( ) 0, . t 0 or t 1 ( ) This function reaches its maximum when t
max1

2

M -3 2

,

=

2 - M +M ( )+

(2 - M )2 +M ( )(-2M +M ( )+8) . 21 ( )

(56)

Since the problem is one-dimensional, the confidence interval 1 ( ) coincides with the interval obASTRONOMY REPORTS Vol. 53 No. 8 2009

tained using the least-squares method for the corresponding confidence level. The above can be generalized to the case of a


732

ABUBEKEROV et al.

model with P parameters by obtaining the distribution of the ratio of the projections of the confidence regions determined using the 2 method and the M 2 method. The distribution of the corresponding P residuals is (57) R(1 ... P , ) - Rmin ( ) 2 , P R(1 ... P , ) 2 , M R
min

Replacing throughout the previous derivations except in (43) and (44) R(, ) with the P -parameter residual minimized over all parameters except for one; and 2 , 2 -1 , 1 ( ),and 1 ( ) with 2 , 2 -P , 1 M P M P ( ),and P ( ), we find that the distribution function is equal to P M ( ) > 0|
2 M -P

(58) (59)

( )

2 M -P

.

M (, ) t P ( )
2

(60) ,

Recall that, in the 2 and 2 methods, by conP M fidence interval, we mean the projection of the confidence region obtained from the main corresponding statistic onto the axis for the corresponding parameter. This projection is found as the distance between the roots of the equation obtained by equating the residual minimized over all parameters except for one with the quantile for the selected confidence level. 1 2( e aP (t) = Ap 0,

=1-

M ( ) - P ( )t 2 M -P (M ( ))

while the density of the distribution of positive values M (, ) (with accuracy to within a normalization P ( ) coefficient Ap )is

P ( )t2 -

M

( ))

t M ( ) - P ( )t M ( ) . P ( )

2

M -P -2 2

,

0
M ( ) , P ( )

t0

or t

The maximum of this function is reached when

t

maxP

=

1+ P - M +M ( )+

(1 + P - M )2 +M ( )(2P - 2M - 2+ M ( )+8) . 2P ( )

(61)

Figure 3 shows the behavior of the density function a1 (t) as a function of the confidence level for the 2 M method. Recall that, in the one-parameter problem, the projection of the confidence region obtained using a statistic that is distributed as 2 coincides with the P confidence interval obtained using the least-squares method. Figure 3 shows that the most probable value ratio of the confidence intervals obtained using the 2 method and the least-squares method depends M on the chosen confidence level. For example, with = 0.68 and = 0.99, the most probable confidence interval obtained using the 2 method exceeds the M confidence interval for the least-squares method by

about a factor of four and a factor of 2.1, respectively. These calculations were carried out for M = 101. Further, we investigated the behavior of the maximum of aP (t) as a function of and the dimension of the 2 distribution. Figure 4 presents a set of these M dependences for M = 50, 100, 500, and 1000. Note that the dependence has a universal character for a linear model; i.e., it is independent of the coefficients of the linear function considered. Figure 4 shows that the excess of the 2 confidence intervals over the M least-squares confidence intervals grows (comparatively weakly) with increasing dimension M of the 2 M distribution (i.e., with the number of observational points M ). It also follows from Fig. 4 that the ratio of
ASTRONOMY REPORTS Vol. 53 No. 8 2009


PARAMETER ERRORS IN INVERSE PROBLEMS

733

a

1

0.99

tmax 8

1000 500

0.8 0.6 0.4 0.2
0.95 0.86

6 4
100 50 10

0.6827 0.5

2 8 10 12 14 M/1

0

2

4

6

0.5

0.6

0.7

0.8

0.9

1.0

Fig. 3. Differential density of the distribution function for the ratio of the confidence interval (2 ) (obtained usM ing the 2 method with an adequate model) to the confiM dence interval obtained using the differential-correction method. The case of a one-dimensional function with M = 101 and = 0.5, 0.6827, 0.86, 0.95, 0.99 is shown.

the most probable 2 confidence interval to the leastM squares confidence interval decreases with increasing . This dependence is shown in Fig. 5. This raises the question of why the confidence intervals in the differential-correction method (which are identical to the intervals for the least squares method for linear models) and Monte Carlo method are appreciably smaller than those obtained for the 2 method. In the one-dimensional case (P = 1), M the problem of projecting the confidence region onto the parameter axis does not arise, and the difference between the confidence intervals is clearly associated with two circumstances. 1. When using the 2 statistic (residual), there is M no initial assumption that the model used is correct. In practice, this means that the surface of the residual functional (Fig. 2) for the 2 statistic for various M realizations of the observational data can be shifted both along the parameter axis (in the "horizontal" direction) and along the axis for the residual functional (in the "vertical" direction). For some realizations of the observational data, the confidence region can even degenerate into an empty set (the number of such cases should be close to the adopted significance level, = 1 - ). The surface of the residual functional (Fig. 2) for the 2 statistic for various P realizations of the observational data can only be shifted along the parameter axis (in the "horizontal" direction). The corresponding confidence region never degenerates into an empty set. 2. For equal conditions, the distribution of 2 M confidence intervals, obtained assuming that the model is adequate and using the 2 method is M
ASTRONOMY REPORTS Vol. 53 No. 8 2009

Fig. 4. Case of a one-parameter problem. Dependence of the maximum of the differential distribution function of the ratio of the confidence interval obtained using a main statistic distributed as 2 to the confidence interval M obtained using the differential-correction method on the confidence level for numbers of observational points M = 10, 50, 100, 500, and 1000. These data were obtained for a one-dimensional linear model function (the form of this function is not important; the dependence shown is universal).

broader (Fig. 3) than the normal distribution of the central parameter values used in the differentialcorrection method (in which the model is assumed to be adequate by definition). The characteristic parameters of these distributions are, in the former case, the most probable 2 error interval and, in M the latter case, the root-mean-square (rms) estimate of a standard deviation multiplied by a coefficient corresponding to the selected confidence level (which coincides with the 2 interval in the one-dimensional P case). The parameter tmax is the ratio of the first to the second of these parameters. In the multi-parameter case (P > 1), the effect of having to project a multi-dimensional (in the spatial parameters) confidence region onto the parameter axes is added to the reasons indicated above. This must be done in order to make the transition from an "exact" confidence region to the values of the confidence intervals. Since the projection of the confidence region is replaced by its volumization as a parallelopiped whose volume is larger, the corresponding probability for overlapping the exact solution additionally grows. Thus, when choosing the confidence level for the confidence intervals, we must understand that the probability P that they encompass the exact solution of the corresponding confidence region is P < in the differential-correction or Monte Carlo methods, while the probability of encompassing the exact solution by the corresponding parallelopiped volumizing the


734

ABUBEKEROV et al.
0.5 0.6827 0.86 0.95 0.99

tmax 8 6 4 2

ciably affect the results, as is indicated by the fact that the estimates of the dispersions of the parameters obtained using the differential-correction and Monte Carlo methods do not substantially differ. The appreciable difference between the intervals obtained using the different methods is due precisely to the use of different statistics, as well as differences in the assumptions concerning the correctness of the model.

5.1. New Method for Estimating the Parameter Errors
200 400 600 800 1000 M

Fig. 5. The one-parameter problem. Dependence of the maximum of the differential distribution function for the ratio of the 2 and differential-correction confidence inM tervals on the number of points M for specified confidence levels = 0.5, 0.6827, 0.86, 0.95,and 0.99.

kptmax 8 6 4 2

0.5 0.6827 0.86 0.95 0.99

200

400

600

800

1000 M

Fig. 6. Multi-parameter case. The quantity kP tmax for the multi-parameter problem as a function of the number of observational points M and of for P = 10. The confidence level is indicated next to the corresponding curve. It is striking that the probability that the true value of the parameter falls into the projection of the 2 confiM dence region is higher than the specified confidence level when = 0.6827: it is equal to 0.80. At the same time, the probability that the true solution falls into the error interval for the differential-correction method is equal to the specified confidence level, = 0.6827. See text for details.

"exact" confidence region is P > in the 2 method. M In the 2 method, the probability of encompassing M the exact solution of the corresponding confidence region is P = . In the case of a non-linear model, there is another reason why the error intervals obtained using the differential-correction method can be reduced: the application of a linearization procedure. However, in the model used by us, this factor does not appre-

Thus, when the dependence of a function describing a model on the unknown parameters is linear, then (56) and (61) enable, based on the Monte Carlo or differential-correction method via multiplication by the coefficient tmax (or kP tmax , see below), the derivation of the size of the error interval for an unknown parameter in the framework of the 2 confidence M regions; this distribution is the most probable among the cases in which the model is not rejected. Having obtained the size of the Monte Carlo or differentialcorrection error interval and convinced ourselves that the model is adequate in the framework of the 2 M method (for which it is sufficient to calculate the minimum corresponding residual), we can draw preliminary conclusions of a probabilistic nature about the size of the 2 confidence interval, without turning M directly to the more labor-intensive determination of the sizes of the projection of the confidence region onto the parameter axis. We emphasize again that tmax is the most probable ratio of the projection of the 2 and 2 conM P fidence regions. As was noted above, in the onedimensional case, the projection of the 2 confidence P region is equal to the confidence integral yielded by the differential-correction method [1]. Therefore, the value of tmax for the one-dimensional case is the most probable ratio of the projection of the 2 confidence M region and the differential-correction confidence interval. Values of tmax for the one-dimensional problem are given in Table 1. Table 1 shows that tmax differs appreciably from unity, reflecting the reduction in the parameter errors obtained using the differential-correction method (a normally distributed statistic). Tables 2­5 present values of tmax for the multiparameter case with P = 3, 6, 10, and 50. In the multi-dimensional problem, the projection of the 2 P confidence region is not equal to the confidence intervals obtained using the differential-correction or Monte Carlo methods [1]. Therefore, for convenience in comparing the ratio of the projection of the 2 M confidence region to the differential-correction confidence interval, Table 6 presents the corresponding
ASTRONOMY REPORTS Vol. 53 No. 8 2009


PARAMETER ERRORS IN INVERSE PROBLEMS

735

Table 1. Parameter tmax for the one-parameter problem as a function of M and for P = 1. (The confidence interval obtained using the 2 method coincides with that for the differential-correction method) P M 0.1 10 20 30 40 50 60 70 80 90 100 200 300 400 500 600 700 800 900 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500 10 000 8.49 10.4 11.7 12.7 13.5 14.1 14.7 15.3 15.8 16.2 19.5 21.7 23.3 24.7 25.9 26.9 27.9 28.7 29.5 32.7 35.2 37.2 39.0 40.5 41.9 43.2 44.4 45.4 46.4 47.4 48.3 49.1 49.9 50.7 51.4 52.1 52.8 0.2 5.20 6.21 6.89 7.41 7.84 8.22 8.54 8.84 9.11 9.35 11.1 12.3 13.3 14.0 14.7 15.3 15.8 16.3 16.7 18.5 19.9 21.0 22.0 22.9 23.6 24.4 25.0 25.6 26.2 26.7 27.2 27.7 28.1 28.6 29.0 29.4 29.7 0.3 3.98 4.68 5.16 5.53 5.84 6.10 6.34 6.55 6.74 6.92 8.20 9.06 9.73 10.3 10.8 11.2 11.6 11.9 12.2 13.5 14.5 15.3 16.1 16.7 17.2 17.8 18.2 18.7 19.1 19.5 19.8 20.2 20.5 20.8 21.1 21.4 21.7 0.4 3.32 3.87 4.25 4.54 4.79 5.00 5.18 5.35 5.50 5.64 6.67 7.36 7.89 8.33 8.71 9.05 9.35 9.63 9.88 10.9 11.7 12.4 13.0 13.5 13.9 14.3 14.7 15.1 15.4 15.7 16.0 16.3 16.5 16.8 17.0 17.3 17.5 0.5 2.89 3.35 3.67 3.92 4.12 4.30 4.45 4.59 4.72 4.84 5.70 6.29 6.74 7.11 7.44 7.72 7.98 8.21 8.42 9.31 9.99 10.6 11.0 11.5 11.9 12.2 12.5 12.8 13.1 13.4 13.6 13.9 14.1 14.3 14.5 14.7 14.9 0.6 2.58 2.98 3.25 3.47 3.64 3.80 3.93 4.06 4.17 4.27 5.02 5.53 5.93 6.25 6.54 6.78 7.01 7.21 7.40 8.17 8.77 9.26 9.68 10.1 10.4 10.7 11.0 11.2 11.5 11.7 11.9 12.1 12.3 12.5 12.7 12.9 13.0 0.7 2.33 2.68 2.92 3.11 3.27 3.41 3.53 3.63 3.73 3.82 4.49 4.94 5.29 5.58 5.83 6.06 6.25 6.43 6.60 7.29 7.82 8.26 8.63 8.97 9.27 9.54 9.79 10.0 10.2 10.4 10.6 10.8 11.0 11.2 11.3 11.5 11.6 0.8 2.11 2.42 2.64 2.81 2.95 3.07 3.17 3.27 3.36 3.44 4.03 4.44 4.75 5.01 5.23 5.43 5.61 5.77 5.92 6.53 7.00 7.39 7.73 8.03 8.30 8.54 8.77 8.98 9.17 9.35 9.53 9.69 9.84 9.99 10.1 10.3 10.4 0.9 1.89 2.16 2.35 2.50 2.62 2.73 2.82 2.90 2.98 3.05 3.57 3.93 4.20 4.43 4.63 4.80 4.96 5.10 5.23 5.77 6.19 6.53 6.83 7.09 7.33 7.54 7.74 7.93 8.10 8.26 8.41 8.56 8.69 8.82 8.95 9.07 9.19 0.95 1.76 2.01 2.18 2.31 2.42 2.52 2.60 2.68 2.75 2.81 3.29 3.61 3.87 4.08 4.26 4.41 4.56 4.69 4.81 5.30 5.68 6.00 6.27 6.52 6.73 6.93 7.11 7.28 7.44 7.58 7.72 7.86 7.98 8.10 8.22 8.33 8.43 0.99 1.58 1.79 1.94 2.05 2.15 2.23 2.31 2.37 2.43 2.49 2.90 3.18 3.40 3.58 3.74 3.88 4.00 4.12 4.22 4.65 4.99 5.27 5.50 5.71 5.90 6.08 6.23 6.38 6.52 6.65 6.77 6.89 7.00 7.10 7.20 7.30 7.39

ASTRONOMY REPORTS

Vol. 53 No. 8 2009


736

ABUBEKEROV et al.

Table 2. Parameter tmax for the multi-parameter problem as a function of M and for P = 3.(The quantity kP tmax must be used to translate to the confidence interval obtained using the differential-correction or Monte Carlo methods; the values of kP were taken from Table 6) M 0.1 10 20 30 40 50 60 70 80 90 100 200 300 400 500 600 700 800 900 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500 10 000 1.71 1.93 2.10 2.24 2.36 2.46 2.56 2.64 2.71 2.78 3.30 3.64 3.91 4.13 4.32 4.49 4.64 4.78 4.91 5.43 5.83 6.17 6.45 6.71 6.93 7.14 7.33 7.51 7.67 7.83 7.97 8.11 8.24 8.37 8.49 8.60 8.72 0.2 1.60 1.78 1.92 2.03 2.13 2.21 2.29 2.36 2.42 2.48 2.91 3.20 3.43 3.62 3.78 3.93 4.05 4.17 4.28 4.73 5.07 5.36 5.60 5.82 6.02 6.19 6.36 6.51 6.65 6.78 6.91 7.03 7.14 7.25 7.35 7.45 7.55 0.3 1.55 1.71 1.84 1.94 2.03 2.11 2.18 2.24 2.30 2.35 2.74 3.01 3.22 3.39 3.54 3.68 3.80 3.90 4.00 4.42 4.74 5.00 5.23 5.43 5.61 5.78 5.93 6.07 6.20 6.32 6.44 6.55 6.65 6.75 6.85 6.94 7.03 0.4 1.51 1.67 1.79 1.89 1.98 2.05 2.11 2.17 2.23 2.28 2.65 2.90 3.10 3.27 3.41 3.54 3.65 3.75 3.85 4.24 4.55 4.80 5.02 5.21 5.38 5.54 5.68 5.82 5.94 6.06 6.17 6.28 6.38 6.48 6.57 6.66 6.74 0.5 1.48 1.64 1.76 1.86 1.94 2.01 2.07 2.13 2.18 2.23 2.59 2.84 3.03 3.19 3.33 3.45 3.56 3.66 3.75 4.13 4.43 4.67 4.88 5.07 5.24 5.39 5.53 5.66 5.78 5.90 6.01 6.11 6.21 6.30 6.39 6.47 6.56 0.6 1.45 1.61 1.73 1.83 1.91 1.98 2.04 2.10 2.15 2.20 2.55 2.79 2.98 3.13 3.27 3.39 3.50 3.59 3.68 4.06 4.35 4.58 4.79 4.97 5.14 5.29 5.42 5.55 5.67 5.78 5.89 5.99 6.08 6.18 6.26 6.35 6.43 0.7 1.42 1.59 1.71 1.80 1.88 1.95 2.01 2.07 2.12 2.16 2.51 2.75 2.93 3.09 3.22 3.34 3.44 3.54 3.63 3.99 4.28 4.51 4.71 4.89 5.05 5.20 5.34 5.46 5.58 5.69 5.79 5.89 5.99 6.07 6.16 6.24 6.32 0.8 1.39 1.56 1.68 1.77 1.85 1.92 1.98 2.03 2.08 2.13 2.47 2.70 2.88 3.03 3.16 3.28 3.38 3.48 3.57 3.93 4.20 4.44 4.64 4.81 4.97 5.11 5.25 5.37 5.48 5.59 5.69 5.79 5.88 5.97 6.06 6.14 6.21 0.9 1.35 1.52 1.63 1.72 1.80 1.87 1.93 1.98 2.03 2.07 2.41 2.64 2.81 2.96 3.09 3.20 3.30 3.39 3.48 3.83 4.10 4.33 4.52 4.70 4.85 4.99 5.12 5.24 5.35 5.46 5.56 5.65 5.74 5.83 5.91 5.99 6.06 0.95 1.32 1.48 1.60 1.69 1.76 1.83 1.89 1.94 1.99 2.03 2.36 2.58 2.75 2.90 3.02 3.13 3.23 3.32 3.40 3.75 4.01 4.23 4.42 4.59 4.74 4.88 5.01 5.13 5.24 5.34 5.44 5.53 5.62 5.70 5.78 5.86 5.93 0.99 1.28 1.43 1.54 1.62 1.69 1.75 1.81 1.86 1.90 1.95 2.26 2.47 2.63 2.77 2.89 2.99 3.09 3.18 3.25 3.58 3.84 4.05 4.23 4.39 4.53 4.66 4.78 4.90 5.00 5.10 5.19 5.28 5.37 5.45 5.52 5.60 5.67
2009

ASTRONOMY REPORTS Vol. 53 No. 8


PARAMETER ERRORS IN INVERSE PROBLEMS

737

Table 3. Parameter tmax for the multi-parameter problem as a function of M and for P = 6.(The quantity kP tmax must be used to translate to the confidence interval obtained using the differential-correction or Monte Carlo methods; the values of kP were taken from Table 6) M 0.1 10 20 30 40 50 60 70 80 90 100 200 300 400 500 600 700 800 900 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500 10 000 1.23 1.22 1.26 1.31 1.35 1.40 1.43 1.47 1.50 1.54 1.78 1.94 2.08 2.19 2.28 2.37 2.44 2.51 2.58 2.84 3.04 3.21 3.36 3.49 3.60 3.71 3.81 3.90 3.98 4.06 4.13 4.20 4.27 4.34 4.40 4.46 4.51 0.2 1.22 1.24 1.28 1.32 1.36 1.40 1.44 1.47 1.50 1.53 1.75 1.91 2.03 2.14 2.23 2.30 2.38 2.44 2.50 2.75 2.95 3.11 3.25 3.37 3.48 3.58 3.67 3.76 3.84 3.91 3.99 4.05 4.12 4.18 4.24 4.29 4.35 0.3 1.21 1.25 1.30 1.34 1.39 1.42 1.46 1.49 1.52 1.55 1.76 1.92 2.04 2.14 2.23 2.31 2.38 2.44 2.50 2.74 2.94 3.09 3.23 3.35 3.46 3.56 3.65 3.74 3.82 3.89 3.96 4.03 4.09 4.15 4.21 4.27 4.32 0.4 1.20 1.26 1.32 1.36 1.41 1.45 1.48 1.51 1.54 1.57 1.79 1.94 2.06 2.16 2.25 2.33 2.40 2.46 2.52 2.77 2.96 3.12 3.25 3.37 3.48 3.58 3.67 3.76 3.84 3.91 3.98 4.05 4.11 4.17 4.23 4.29 4.34 0.5 1.19 1.27 1.33 1.38 1.43 1.47 1.50 1.54 1.57 1.60 1.82 1.97 2.09 2.19 2.28 2.36 2.43 2.50 2.56 2.80 3.00 3.16 3.29 3.42 3.53 3.63 3.72 3.80 3.88 3.96 4.03 4.10 4.16 4.22 4.28 4.34 4.39 0.6 1.18 1.27 1.34 1.39 1.44 1.48 1.52 1.56 1.59 1.62 1.84 2.00 2.12 2.23 2.32 2.40 2.47 2.54 2.60 2.85 3.04 3.20 3.34 3.47 3.58 3.68 3.77 3.86 3.94 4.02 4.09 4.16 4.22 4.29 4.35 4.40 4.46 0.7 1.17 1.27 1.35 1.40 1.45 1.50 1.54 1.57 1.61 1.64 1.87 2.03 2.16 2.26 2.35 2.44 2.51 2.58 2.64 2.89 3.09 3.26 3.40 3.52 3.64 3.74 3.84 3.92 4.01 4.08 4.16 4.23 4.29 4.36 4.42 4.48 4.53 0.8 1.16 1.27 1.35 1.41 1.46 1.51 1.55 1.59 1.62 1.65 1.89 2.06 2.19 2.30 2.39 2.47 2.55 2.61 2.68 2.94 3.14 3.31 3.45 3.58 3.70 3.80 3.90 3.99 4.07 4.15 4.23 4.30 4.36 4.43 4.49 4.55 4.61 0.9 1.15 1.27 1.35 1.41 1.47 1.52 1.56 1.60 1.63 1.67 1.91 2.08 2.21 2.32 2.42 2.50 2.58 2.65 2.71 2.98 3.18 3.35 3.50 3.63 3.75 3.86 3.95 4.05 4.13 4.21 4.29 4.36 4.43 4.49 4.56 4.62 4.67 0.95 1.14 1.26 1.34 1.41 1.46 1.51 1.56 1.60 1.63 1.66 1.91 2.08 2.22 2.33 2.42 2.51 2.59 2.66 2.72 2.99 3.19 3.37 3.52 3.65 3.77 3.87 3.97 4.06 4.15 4.23 4.31 4.38 4.45 4.51 4.58 4.64 4.70 0.99 1.13 1.24 1.33 1.39 1.45 1.50 1.54 1.58 1.62 1.65 1.90 2.07 2.20 2.31 2.41 2.49 2.57 2.64 2.70 2.97 3.18 3.35 3.50 3.63 3.75 3.85 3.95 4.04 4.13 4.21 4.29 4.36 4.43 4.49 4.55 4.61 4.67

ASTRONOMY REPORTS

Vol. 53 No. 8 2009


738

ABUBEKEROV et al.

Table 4. Parameter tmax for the multi-parameter problem as a function of M and for P = 10. (The quantity kP tmax must be used to translate to the confidence interval obtained using the differential-correction or Monte Carlo methods; the values of kP were taken from Table 6) M 0.1 10 20 30 40 50 60 70 80 90 100 200 300 400 500 600 700 800 900 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500 10 000 1.16 1.08 1.05 1.06 1.07 1.08 1.09 1.11 1.13 1.14 1.28 1.38 1.46 1.53 1.59 1.64 1.69 1.74 1.78 1.95 2.08 2.20 2.29 2.38 2.46 2.53 2.59 2.65 2.71 2.76 2.81 2.85 2.90 2.94 2.98 3.02 3.06 0.2 1.13 1.10 1.10 1.11 1.12 1.14 1.15 1.17 1.18 1.20 1.33 1.42 1.50 1.57 1.63 1.68 1.73 1.78 1.82 1.99 2.12 2.23 2.33 2.41 2.49 2.56 2.62 2.68 2.74 2.79 2.84 2.89 2.93 2.97 3.01 3.05 3.09 0.3 1.12 1.12 1.13 1.15 1.16 1.18 1.20 1.21 1.23 1.25 1.38 1.48 1.56 1.62 1.68 1.74 1.79 1.83 1.87 2.04 2.18 2.29 2.39 2.47 2.55 2.62 2.69 2.75 2.80 2.86 2.91 2.96 3.00 3.04 3.09 3.13 3.16 0.4 1.10 1.13 1.15 1.17 1.19 1.21 1.23 1.25 1.27 1.29 1.43 1.53 1.61 1.68 1.74 1.79 1.84 1.89 1.93 2.11 2.24 2.36 2.46 2.55 2.63 2.70 2.77 2.83 2.89 2.94 2.99 3.04 3.09 3.13 3.17 3.22 3.26 0.5 1.09 1.13 1.16 1.19 1.22 1.24 1.27 1.29 1.31 1.32 1.47 1.58 1.66 1.73 1.80 1.85 1.91 1.95 2.00 2.18 2.32 2.44 2.54 2.63 2.71 2.78 2.85 2.92 2.98 3.03 3.09 3.14 3.18 3.23 3.27 3.32 3.36 0.6 1.08 1.14 1.18 1.21 1.24 1.27 1.29 1.32 1.34 1.36 1.51 1.62 1.71 1.79 1.85 1.91 1.97 2.01 2.06 2.25 2.39 2.52 2.62 2.71 2.80 2.88 2.95 3.01 3.07 3.13 3.19 3.24 3.29 3.34 3.38 3.43 3.47 0.7 1.08 1.14 1.19 1.22 1.26 1.29 1.32 1.34 1.36 1.39 1.55 1.67 1.76 1.84 1.91 1.97 2.03 2.08 2.12 2.32 2.47 2.60 2.71 2.80 2.89 2.97 3.04 3.11 3.18 3.24 3.29 3.35 3.40 3.45 3.49 3.54 3.58 0.8 1.07 1.14 1.19 1.24 1.27 1.31 1.34 1.36 1.39 1.41 1.59 1.71 1.81 1.90 1.97 2.03 2.09 2.14 2.19 2.39 2.55 2.68 2.80 2.90 2.99 3.07 3.15 3.22 3.28 3.35 3.41 3.46 3.51 3.57 3.61 3.66 3.71 0.9 1.06 1.14 1.20 1.25 1.29 1.32 1.36 1.39 1.41 1.44 1.63 1.76 1.86 1.95 2.03 2.09 2.15 2.21 2.26 2.47 2.64 2.77 2.89 3.00 3.09 3.18 3.26 3.33 3.40 3.47 3.53 3.59 3.64 3.69 3.75 3.79 3.84 0.95 1.05 1.14 1.20 1.25 1.30 1.33 1.37 1.40 1.43 1.45 1.65 1.78 1.89 1.98 2.06 2.13 2.19 2.25 2.30 2.52 2.69 2.83 2.95 3.06 3.15 3.24 3.32 3.40 3.47 3.54 3.60 3.66 3.72 3.77 3.82 3.87 3.92 0.99 1.04 1.13 1.20 1.25 1.30 1.34 1.37 1.40 1.43 1.46 1.66 1.81 1.92 2.01 2.09 2.16 2.22 2.28 2.34 2.56 2.73 2.88 3.00 3.11 3.21 3.30 3.39 3.46 3.54 3.60 3.67 3.73 3.79 3.84 3.90 3.95 4.00
2009

ASTRONOMY REPORTS Vol. 53 No. 8


PARAMETER ERRORS IN INVERSE PROBLEMS

739

Table 5. Parameter tmax for the multi-parameter problem as a function of M and for P = 50. (The quantity kP tmax must be used to translate to the confidence interval obtained using the differential-correction or Monte Carlo methods; the values of kP were taken from Table 6) M 0.1 50 60 70 80 90 100 200 300 400 500 600 700 800 900 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500 10 000 1.03 1.01 1.00 0.991 0.982 0.973 0.913 0.881 0.864 0.856 0.852 0.852 0.854 0.857 0.861 0.886 0.913 0.939 0.964 0.987 1.01 1.03 1.05 1.06 1.08 1.10 1.11 1.13 1.14 1.15 1.17 1.18 1.19 0.2 1.02 1.02 1.01 1.00 1.00 0.996 0.967 0.953 0.947 0.945 0.946 0.948 0.952 0.956 0.961 0.989 1.02 1.04 1.07 1.09 1.12 1.14 1.16 1.17 1.19 1.21 1.22 1.24 1.25 1.27 1.28 1.30 1.31 0.3 1.02 1.02 1.02 1.01 1.01 1.01 1.00 0.998 1.00 1.00 1.01 1.02 1.02 1.03 1.04 1.07 1.10 1.13 1.16 1.18 1.21 1.23 1.25 1.27 1.29 1.31 1.32 1.34 1.36 1.37 1.39 1.40 1.41 0.4 1.02 1.02 1.02 1.02 1.02 1.02 1.03 1.03 1.04 1.05 1.06 1.07 1.08 1.09 1.10 1.14 1.18 1.21 1.24 1.27 1.29 1.32 1.34 1.36 1.38 1.40 1.42 1.44 1.45 1.47 1.48 1.50 1.51 0.5 1.02 1.02 1.02 1.03 1.03 1.03 1.05 1.06 1.08 1.09 1.11 1.12 1.13 1.14 1.15 1.20 1.25 1.28 1.32 1.35 1.37 1.40 1.43 1.45 1.47 1.49 1.51 1.53 1.55 1.56 1.58 1.60 1.61 0.6 1.02 1.02 1.03 1.03 1.03 1.04 1.07 1.09 1.11 1.13 1.15 1.16 1.18 1.19 1.20 1.26 1.31 1.35 1.39 1.42 1.45 1.48 1.51 1.53 1.56 1.58 1.60 1.62 1.64 1.66 1.68 1.70 1.71 0.7 1.02 1.02 1.03 1.03 1.04 1.04 1.08 1.11 1.14 1.16 1.18 1.20 1.22 1.24 1.25 1.32 1.37 1.42 1.46 1.50 1.53 1.57 1.60 1.62 1.65 1.67 1.70 1.72 1.74 1.76 1.78 1.80 1.82 0.8 1.02 1.02 1.03 1.04 1.05 1.05 1.10 1.14 1.17 1.20 1.22 1.25 1.27 1.29 1.30 1.38 1.44 1.49 1.54 1.58 1.62 1.66 1.69 1.72 1.75 1.77 1.80 1.82 1.85 1.87 1.89 1.91 1.93 0.9 1.02 1.03 1.04 1.04 1.05 1.06 1.12 1.17 1.21 1.24 1.27 1.30 1.32 1.34 1.37 1.45 1.52 1.58 1.64 1.68 1.72 1.76 1.80 1.83 1.87 1.90 1.92 1.95 1.98 2.00 2.02 2.05 2.07 0.95 1.01 1.03 1.04 1.05 1.06 1.07 1.14 1.19 1.24 1.27 1.30 1.33 1.36 1.38 1.41 1.50 1.58 1.65 1.70 1.75 1.80 1.84 1.88 1.92 1.95 1.98 2.01 2.04 2.07 2.09 2.12 2.14 2.17 0.99 1.01 1.03 1.04 1.05 1.07 1.08 1.16 1.22 1.27 1.32 1.35 1.39 1.42 1.44 1.47 1.58 1.67 1.74 1.80 1.86 1.91 1.95 2.00 2.04 2.07 2.11 2.14 2.18 2.21 2.24 2.26 2.29 2.32

ASTRONOMY REPORTS

Vol. 53 No. 8 2009


740

ABUBEKEROV et al.

Table 6. Recalculation coefficient kP ( ) for the parameter tmax from Tables 2­5 in the ratio of the projection of the confidence region obtained for the 2 method to the confidence interval obtained for the differential-correction or Monte M Carlo methods P 0.1 1 3 6 10 50 1.00 6.08 11.8 17.6 48.9 0.2 1.00 3.96 6.92 9.81 25.4 0.3 1.00 3.10 5.08 7.00 17.3 0.4 1.00 2.61 4.08 5.49 13.1 0.5 1.00 2.28 3.43 4.53 10.4 0.6 1.00 2.04 2.96 3.85 8.56 0.7 1.00 1.85 2.59 3.31 7.14 0.8 1.00 1.68 2.28 2.86 5.95 0.9 1.00 1.52 1.98 2.43 4.83 0.95 1.00 1.43 1.81 2.18 4.19 0.99 1.00 1.31 1.59 1.87 3.39

P ( ) [this ex1 ( ) pression follows from (47) and is analogous to the expression for P ( )], which must be applied to the values of tmax from Tables 2­5. As is described for the one-dimensional case above, the quantity kP tmax in the multi-dimensional case is the ratio of the projection of the 2 confidence region to the differentialM correction or Monte Carlo confidence interval. Table 7 and Fig. 6 present values of kP tmax for the case P = 10. The character of the numerical values of kP tmax from Table 7 is very close to the values of tmax from Table 1, which shows the ratio of the 2 to the M differential-correction error intervals. A comparison of Fig. 5 and Fig. 6 illustrates this clearly. Thus, in the multi-parameter case, the numerical ratios of the 2 and differential-correction error intervals remain M as before (their exact values are obtained from (56) and (61); see also Tables 2­6). Recall that, in the one-dimensional case, the probabilities that the true parameter value falls in the projection of the 2 confidence region and in the P differential-correction confidence interval are equal. In the multi-dimensional case, the probability that the true parameter value falls into the projection of the 2 P confidence region is higher than the corresponding probability for the differential-correction confidence interval. Note that the results presented above were obtained for a linear model. The use of a 2 statistic P to construct the confidence intervals was based on the fact that the minimum residual is distributed as 2 -P . However, this assertion is valid only for a M model that depends linearly on its parameters. When this dependence is not linear, this assertion can be valid only asymptotically, as M [5, 8]. Therefore, the results obtained can also be applied to nonlinear models, but only approximately, in an asymptotic sense. values of the coefficient kP ( ) =

5.2. Significance Level Corresponding to the Reduced Chi-Squared Statistic
The adequacy of a model to describe observational data is often judged based on the nearness to unity of ^ the minimum value of the residual R R/(M - P ), 2 which is distributed as M -P (as is the random quan^ , if 2 -P ), where P is the number tity M M -P of parameters over which the minimization is carried out (an alternative criterion is nearness of the residual Rmin , distributed as 2 -P , to the value M - P ). This M omits the quantitative characteristic of the adequacy of the model, namely the significance level 0 = 1 - at which the model can be rejected (recall that 0 is the number of errors of the first kind that will occur if we reject the model). Figure 7 presents a graphical depiction of the solution of the equation
M -P

(1 - 0 ) =q M -P

for the significance level 0 as a function of the difference M - P (where M is the number of observational points, P the number of unknown parameters, and M -P ( ) a function of the dependence of the quantile on the confidence level = 1 - 0 ) for a specified quantile q for the 2 -P distribution. The value of q ^M was taken to be 0.5, 0.6, 0.7, 1.0, 1.1, 1.5, and 2.0. This dependence can be written in the explicit form 0 (M - P ) = 1 -
2 M -P

(q (M - P )).

Figure 7 shows that 0 begins to approach its asymptotic limit fairly rapidly (for M 50). For example, for q = 1, the asymptotic limit of 0 is 50%; i.e., in this case, we will make an error of the first kind (rejecting a correct model) in 50% of cases if we reject the model. Consequently, there is no serious basis to reject the model, and the model should therefore be accepted.
ASTRONOMY REPORTS Vol. 53 No. 8 2009


PARAMETER ERRORS IN INVERSE PROBLEMS Table 7. Parameter kP t M 0.1 10 20 30 40 50 60 70 80 90 100 200 300 400 500 600 700 800 900 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500 10 000 20.4 18.9 18.5 18.5 18.7 18.9 19.2 19.5 19.8 20.0 22.4 24.2 25.6 26.8 27.9 28.9 29.7 30.5 31.2 34.2 36.6 38.6 40.2 41.7 43.1 44.3 45.5 46.5 47.5 48.4 49.3 50.1 50.9 51.7 52.4 53.1 53.7 0.2 11.1 10.8 10.8 10.9 11.0 11.2 11.3 11.5 11.6 11.8 13.0 14.0 14.7 15.4 16.0 16.5 17.0 17.4 17.8 19.5 20.8 21.9 22.8 23.7 24.4 25.1 25.7 26.3 26.9 27.4 27.9 28.3 28.8 29.2 29.6 30.0 30.3 0.3 7.82 7.81 7.90 8.01 8.13 8.26 8.38 8.50 8.61 8.72 9.63 10.3 10.9 11.4 11.8 12.2 12.5 12.8 13.1 14.3 15.2 16.0 16.7 17.3 17.8 18.3 18.8 19.2 19.6 20.0 20.3 20.7 21.0 21.3 21.6 21.9 22.1 0.4 6.07 6.18 6.31 6.44 6.56 6.67 6.78 6.88 6.98 7.08 7.83 8.39 8.84 9.22 9.56 9.86 10.1 10.4 10.6 11.6 12.3 13.0 13.5 14.0 14.4 14.8 15.2 15.5 15.8 16.1 16.4 16.7 17.0 17.2 17.4 17.7 17.9 0.5 4.96 5.13 5.27 5.41 5.53 5.63 5.74 5.83 5.92 6.00 6.66 7.14 7.53 7.86 8.15 8.40 8.63 8.85 9.04 9.86 10.5 11.0 11.5 11.9 12.3 12.6 12.9 13.2 13.5 13.7 14.0 14.2 14.4 14.6 14.8 15.0 15.2
max

741

for the multi-parameter problem as a function of M and for P = 10 (see also Fig. 6) 0.6 4.17 4.37 4.52 4.65 4.77 4.87 4.97 5.06 5.14 5.22 5.82 6.25 6.59 6.88 7.13 7.36 7.56 7.75 7.92 8.64 9.20 9.67 10.1 10.4 10.8 11.1 11.3 11.6 11.8 12.0 12.3 12.5 12.6 12.8 13.0 13.2 13.3 0.7 3.56 3.77 3.93 4.06 4.17 4.27 4.36 4.44 4.52 4.59 5.14 5.53 5.84 6.10 6.33 6.53 6.71 6.88 7.04 7.68 8.18 8.60 8.97 9.29 9.58 9.84 10.1 10.3 10.5 10.7 10.9 11.1 11.3 11.4 11.6 11.7 11.9 0.8 3.05 3.26 3.41 3.54 3.65 3.74 3.83 3.90 3.97 4.04 4.55 4.90 5.19 5.42 5.63 5.81 5.98 6.13 6.27 6.85 7.30 7.68 8.00 8.29 8.55 8.79 9.00 9.21 9.40 9.57 9.74 9.90 10.1 10.2 10.3 10.5 10.6 0.9 2.57 2.77 2.92 3.03 3.13 3.22 3.30 3.37 3.44 3.50 3.96 4.28 4.53 4.74 4.93 5.09 5.24 5.37 5.49 6.01 6.41 6.74 7.03 7.29 7.52 7.73 7.92 8.10 8.27 8.43 8.58 8.72 8.85 8.98 9.10 9.22 9.34 0.95 2.29 2.49 2.62 2.73 2.83 2.91 2.98 3.05 3.11 3.17 3.60 3.89 4.13 4.32 4.49 4.64 4.78 4.90 5.02 5.49 5.86 6.17 6.44 6.67 6.88 7.08 7.25 7.42 7.57 7.72 7.86 7.99 8.11 8.23 8.34 8.45 8.55 0.99 1.95 2.12 2.24 2.34 2.43 2.50 2.57 2.63 2.68 2.73 3.11 3.38 3.58 3.76 3.91 4.04 4.16 4.27 4.37 4.79 5.11 5.38 5.62 5.82 6.01 6.18 6.33 6.48 6.61 6.74 6.86 6.98 7.09 7.19 7.29 7.38 7.48

ASTRONOMY REPORTS

Vol. 53 No. 8 2009


742

ABUBEKEROV et al.

0 = 1 ­ 1.0 0.8

0.5 0.6 0.7 0.8

0.6 0.4 0.2

0.9 1.0 1.1 1.2 1.3 1.4 1.5 2.0

20

40 M­P

60

80

100

(1 - 0 ) = q for the significance level 0 for P = 3 as a function of the number M -P of observational points M for a specified quantile q of the reduced statistic 2 -P . The value of q is indicated next to the ^M corresponding curve. When the minimum value of the reduced chi-squared is equal to unity, the corresponding significance level 0 approaches 0.5. Fig. 7. Solution of the equation
M -P

At the same time, when q = 1.5, the asymptotic limit for 0 is already close to 0%. In this case, we will make an error of the first kind in 0% of cases if we reject the model. Thus, here, there is no basis to accept the model, and the model should be rejected. Thus, we conclude that, if the reduced chi-squared for the solution of an inverse parametric problem is 1.5, the model used is clearly inadequate to the observational data. The resulting values of the unknown parameters and their confidence intervals (errors) are "poor", and it is necessary to construct a better model. 6. FITTING THE OBSERVED LIGHT CURVE OF THE YZ Cas BINARY SYSTEM Here, we analyze the normalized light curve of the eclipsing binary system YZ Cas obtained in a red filter ° ( = 6700 A) by Kron [6]. The observed light curve included M = 42 points, 1 ... 42 . The dispersions of each brightness measurement in magnitudes were obs obs taken to be (m )2 = 1.77 â 10-6 [8], or (m )2 = -6 on our (intensity) scale, taking into 1.5015 â 10 account that Lfull = 1. The central value of each point in the light curve was obtained by averaging Nm = 12 points for m = 1 ... 42 [3, 6]. We supplemented the fits of the observed light curve of YZ Cas presented in [1] with fits for the limb-darkening coefficient assuming a linear limbdarkening law. Thus, our fitting of the light curve was carried out over six parameters: r1 , r2 , i, x1 , x2 ,

I0 /I0 . The minimization of the residual R between the observed and theoretical light curves was carried out using a gradient-descent method. The central values and confidence intervals of the unknown parameters are presented in Table 8. The calculation of the central values and error intervals of r1 ,r2 ,i in the differential-correction and confidence-region methods is described in detail in [1], and we do not present this information here. Note that we focus here specifically on the calculation of the central values and confidence intervals of the limb-darkening coefficients, together with verification of the reliability of the resulting confidence intervals. The errors in the limb-darkening coefficients in the differential-correction method were obtained using (41). Since the use of statistics distributed as 2 and P 2 requires knowledge of the unit-weight dispersion M 0 , which is not precisely known for a real observed light curve, we analyzed a model binary system whose phases 1 ...42 coincided with the observed phases for YZ Cas. The true parameter values for the binary system were set equal to the central values obtained from the fit of the observed light curve using a linear limb-darkening model [1], i.e., r1 = 0.14408, ¯ i r2 = 0.07556, ¯ = 88.27 , x1 = 0.2998, x2 = 0.4071. ¯ The unit-weight dispersion was set equal to the rms estimate 0 obtained from this fitting of the observed light curve, 0 = 0.003422904. The weighting coefficients wm and number of measurements at each
ASTRONOMY REPORTS Vol. 53 No. 8 2009

(1)

(2)


PARAMETER ERRORS IN INVERSE PROBLEMS Table 8. Fitting the observed light curve of YZ Cas with a linear limb-darkening law Method Differential correction (est ) Confidence regions, FM,N -M ( = 68.2%) r1 r2 i,deg x1 x2

743

0.14408 ± 0.00023 0.07556 ± 0.00038 88.27 ± 0.090 0.2998 ± 0.02115 0.4071 ± 0.1172 0.1442 ± 0.00217 0.07554 ± 0.00114 88.37 ± 0.57 0.2917 ± 0.1199 0.1959 ± 0.7666

Table 9. Fitting the model "observed" light curve of YZ Cas with a linear limb-darkening law. (The unit-weight rms is taken to be 0 = 0.003422904) Differentialcorrection method 0.14422 ± 0.00023 0.07557 ± 0.00038 88.28 ± 0.091 0.3052 ± 0.01942 0.3759 ± 0.1184 Confidence-region method, 2 P 0.1439 ± 0.000641 0.07538 ± 0.00103 88.42 ± 0.27 0.3013 ± 0.05321 0.3403 ± 0.3280 Confidence-region method, 2 M 0.1439 ± 0.0007767 0.07539 ± 0.001247 88.43 ± 0.33 0.3002 ± 0.06488 0.3213 ± 0.4015 Confidence-region method, Fisher distribution F 0.1445 ± 0.00070 0.07564 ± 0.0011 88.19 ± 0.30 0.3007 ± 0.0589 0.3316 ± 0.3637

Parameter r1 r2 i,deg x1 x2

phase N1 ... N42 we set equal to 12. Further, we call this the model observed system. The differential-correction method [using (31 ))] and confidence-region method based on 2 , 2 , M P and Fisher statistics (FM,N -M in the case with x1 , x2 r1 , r2 , i, I0 /I0 , and FM -3,N -M in the case with r1 , r2 , i) were used to obtain the parameter error intervals for the model observed system, as onedimensional projections of the confidence region at the confidence level = 0.6827. Note that the error intervals of the parameters obtained for the model observed light curve differ from those for the observed light curve from [6] due to the different realizations of the light curve itself (in spite of the fact that the unit-weight dispersion 0 for the model curve was set equal to the rms estimate for the observed curve, 0 = 0.003422904). Table 9 presents the results of our fitting of the model observed light curve of YZ Cas. Note that the error intervals obtained for the confidence-region method depend on the realization of the model observed light curve. The results presented in Table 9 were all obtained using the same realization of the model obsered light curve of YZ Cas. Figures 8 and 9 present the residuals obtained using the 2 method to calculate the limb-darkening M coefficients x1 and x2 (the residuals were minimized (1) (2) over all remaining parameters, r1 , r2 , i,and I0 /I0 ) based on the model observed light curve. The model
ASTRONOMY REPORTS Vol. 53 No. 8 2009
(1) (2)

can be rejected only at a high significance level, 0 = 0.3173 (the corresponding minimum reduced chisquared is 1.1; Fig. 7). Therefore, the model can be considered to be adequate to the observational data, and we can determine the confidence region for the parameters at the confidence level = 0.6827, which is not empty. However, the region of limbdarkening coefficients for the smaller star x2 satisfying the model extends to negative values (the value of x2 is contained in the interval from -0.0802 to 0.7228), in contradiction with the physical range of these values, x2 > 0. Thus, the given accuracy of the observed light curve (0 = 0.003423) does not enable reliable calculation of the linear limb-darkening coefficient for the smaller star based on an analysis of the secondary minimum in the light curve (total eclipse) of the YZ Cas binary system. The limb-darkening coefficient for the larger star (annular eclipse) is determined more certainly: x1 = 0.3002 ± 0.06488. As follows from Table 9, the errors in x1 and x2 found using the differential-correction method (using a normally distributed statistic) are appreciably smaller than the errors found using the 2 method M (in this case, the projections of the confidence region onto the x1 and x2 parameter axes). When choosing the errors derived from the differential-correction method, we must bear in mind that these represent only "internal" errors, which can be appreciably smaller than the "external" errors in parameters found


744

ABUBEKEROV et al.

R 80 70 60 50 ­ 0.2 0.2 0.4 0.6 0.8 M() P() 1.0 x2

Fig. 8. Model observed light curve for the binary system YZ Cas. The residual R for the 2 method (M = 42) M for the limb-darkening coefficient x2 of the secondary (star with a smaller radius), minimized over all remaining 1 2 parameters (r1 , r2 , i, x1 , I0 /I0 ). The rms estimated unitweight dispersion of the light curve corresponds to 0 = 0.003423. The straight lines show the quantiles for 2 M and 2 distributions for the confidence level = 0.6827 P (when using 2 statistics, the quantity P ( ) is meaP sured from the minimum residual).

obs m / wm = 0.0003422904 (where the weight of each point is wm = Nm = 12). Further, we fit the resulting perturbed light curve. Figure 10 presents the behavior of the residual obtained using the 2 method to calculate the limbM darkening coefficient x2 . At the chosen confidence level, = 0.6827, the model is adequate to the observational data, with x2 = 0.4205 ± 0.0399. Thus, to reliably determine the limb-darkening coefficients in classical binaries whose smaller components have radii similar to that of the secondary component in YZ Cas, the root square from unit-weight dispersion must be close to 0 10-4 . Such high-accuracy light curves can be obtained with the COROT and KEPLER orbiting photometric stations. Analysis of the main minimum (annular eclipse) indicates that the limb-darkening coefficient of the larger star x1 is found reliably with 0 = 0.0034: x1 = 0.3002 ± 0.06488 for the 2 method, where this indiM cates the projection of the six-dimensional confidence region onto the x1 axis.

6.1. Use of a Quadratic Limb-Darkening Law
from an analysis of another realization of the light curve (for example, obtained at another epoch). Our numerical simulations showed that reliably determining the limb-darkening coefficient of the smaller star in the YZ Cas binary system requires increasing the accuracy of the observed light curve by approximately an order of magnitude. We perturbed the light curve constructed with the "true" parameter values, for which we adopted the central ¯ values from Table 8: r1 = 0.14408, r2 = 0.07556, ¯ ¯ = 88.27 , x1 = 0.2998, x2 = 0.4071. Here, the uniti weight dispersion root square was taken to be 0 =
R 65 60 55 50 45 40 0 0.25 0.30 0.35 0.40 x1 M() P()

We also fit the observed light curve assuming a quadratic limb-darkening law [see (3)]. Here, we obtained two different fits. In the first case, the unknown parameters were r1 , r2 , i, x1 , x2 , y1 ,and y2 ; in the second, the coefficients x1 and x2 were fixed. We used the values of Van Hamme [10] for x1 and x2 , and the limb-darkening law for the smaller component was assumed to be linear. The fitting amounted to searching for the parameters r1 , r2 , i,and y1 . The results of this fitting are presented in Table 10. The errors were obtained using the differentialcorrection method (for the confidence level 0.6827). When fitting the observed light curve over the entire collection of parameters, the values of x1 , x2 , y1 , and
R 52.5 50.0 47.5 45.0 42.5 0.36 0.38 37.5 0.42 0.44 x2
M() P()

Fig. 9. Same as Fig. 8 for the primary component (star with the larger radius) and minimization over the param1 2 eters r1 , r2 , i, x2 , I0 /I0 .

Fig. 10. Same as Fig. 8 for the perturbed light curve. See text for details. ASTRONOMY REPORTS Vol. 53 No. 8 2009


PARAMETER ERRORS IN INVERSE PROBLEMS

745

Table 10. Results of fitting the observed light curve of YZ Cas using a quadratic limb-darkening law (0 = 0.003422904) r1 0.144324 ± 0.000308 r2 0.753775 ± 0.000220 i,deg 88.22 ± 0.099 x1 0.331 x2 0.383 y
1

-0.0134 ± 0.0100

Table 11. Half the error intervals for the limb-darkening coefficients for the YZ Cas components derived by fitting the model observed light curve using a linear limb-darkening law, and the number of times the "true" values fall in these intervals, Nx1 and Nx2 .(N is the number of times the exact solution falls in the error interval. 0 = 0.003422904.The number of times the exact solution falls in the confidence region D is indicated in parantheses.) Method Differential-correction Monte Carlo Confidence region, Confidence region,
2 P 2 M N -M

N

x

1

x1 ± x1 0.3052 ± 0.01942 0.3009 ± 0.02011 0.3013 ± 0.05321 0.3002 ± 0.06488 0.3007 ± 0.0589

N

x

2

x2 ± x2 0.3759 ± 0.1184 0.4934 ± 0.1056 0.3403 ± 0.3280 0.3213 ± 0.4015 0.3316 ± 0.3637
(1) (2)

N 4916 4871 9634 (6898) 8166 (6858) 8060 (6827)

6892 6899 9722 8234 8112

6957 6884 9842 8413 8306

Confidence region, FM,

y2 lie outside the physically reasonable range and are thus artefacts, although the best-fit light curve for these values fits the observed light curve well, and is very close to the light curve obtained for a linear limbdarkening law. In the second case, y1 = -0.0134 ± 0.0100. The best-fit light curve is nearly identical (with accuracy to the third digit after the decimal point) to the bestfit light curve obtained using a linear limb-darkening law (Fig. 11). Thus, with 0 = 0.0034, it is fully sufficient to restrict our analysis to a linear limb-darkening law when fitting the observed light curve. The corrections for non-linearity of this law are negligibly small. The test of reliability of the error intervals for the limbdarkening coefficients is satisfied for a linear law.

with x1 , x2 r1 , r2 , i, I0 /I0 , and FM -3,N -M in the case with r1 , r2 , i) by solving the corresponding (1) (2) inequalities for r1 , r2 , i, x1 , x2 ,and I0 /I0 . Table 11 shows that the number of times the exact solution falls in the individual error intervals for the differential-correction and Monte Carlo methods corresponds to the probability = 68.2%, but the number of joint "hits" is less than this probability (N 4900, which is 30% lower than the value N 6800 corresponding to = 0.6827). In this case, the confidence intervals for x1 and x2 are projections of
L 1.00 0.95 0.90 0.85 0.80 0.75 0 2 4 6 8 10 12 14

6.2. Test of Reliability of the Limb-Darkening Coefficients and Their Errors We verified the trustworthiness of the error intervals for x1 and x2 calculated for each method. We carried out 10 000 perturbations of the model observed light curve (for 0 = 0.003422904) calculated with the "true" parameter values, taken to be the cen¯ tral values from Table 8: r1 = 0.14408, r2 = 0.07556, ¯ ¯ = 88.27 , x1 = 0.2998, x2 = 0.4071. We counted i the number of times the "true" values fell in the corresponding confidence regions or their projections (Table 11). The confidence intervals for x1 and x2 calculated as the projection onto the x1 , x2 parameter axes of the six-dimensional confidence regions obtained using 2 , 2 , and Fisher statistics (FM,N -M in the case P M
ASTRONOMY REPORTS Vol. 53 No. 8 2009

Fig. 11. Observed light curve for YZ Cas (annular eclipse) and the best-fit theoretical light curves obtained for linear and quadratic limb-darkening laws (the curves coincide with accuracy to within the third digit after the decimal place). The points show the observational data from [6].


746

ABUBEKEROV et al.

x1 0.06 0.04 0.02 ­ 0.0010 ­ 0.0005 ­ 0.02 ­ 0.04 ­ 0.06 0.0005 0.0010 r1

l 1.00 0.95 0.90 0.85 0.80 0.75 0 2 4 6 8 10 12 14°

Fig. 12. Projection of the six-dimensional confidence re(1) (2) gion for the parameters r1 , r2 , i, x1 , x2 ,and I0 /I0 onto the plane of the parameters r1 , x1 ,obtained by fitting the model observed light curve for YZ Cas using statistics distributed as 2 (small ellipse), the Fisher distribution P FM,N -M (medium ellipse), and 2 (large ellipse). The M rectangle in the center is the confidence region obtained using the differential-correction or Monte Carlo method. The deviations of the parameters from the central values are shown. The adopted unit-weight dispersion square root is 0 = 0.003423. The confidence intervals obtained using the differential-correction and Monte Carlo methods (using a normally distributed statistic) cover the exact values of the parameters with the specified probability, = 0.68, while the corresponding confidence region (central rectangle) covers the exact solution with a probability that is lower than , equal to 0.49. The confidence regions D obtained using 2 , FM,N -M , and P 2 statistics cover the exact solution with probability M = 0.68, while the corresponding confidence intervals (projections of the confidence region onto the r1 and x1 parameter axes) cover the exact parameter values with a probabilities higher than , namely 0.9904 and 0.9722 for the 2 statistic, 0.8456 and 0.8234 for the 2 , and P M 0.8368 and 0.8112 for the FM,N -M statistic.

Fig. 13. Primary minimum (annular eclipse) of the model observed light curve of YZ Cas together with theoretical light curves obtained for the maximum and minimum values of r1 , r2 , x1 , and x2 determined using the 2 method (outer curves), 2 method (midM P dle curves), and differential-correction method (inner curves). The adopted unit-weight dispersion square root is 0 = 0.003423.

the six-dimensional confidence region for r1 , r2 , i, x1 , (1) (2) x2 , and I0 /I0 onto the x1 and x2 parameter axes. Falling into the confidence region D fully corresponds to the specified confidence level = 68.2%.Since the two-dimensional projection onto the (x1 , x2 ) plane is an elongated "ellipse," number of coincidences with the projection of the "ellipse" onto the x1 and x2 parameter axes appreciably exceeds the confidence level = 68.2%. The situation when the 2 method P is used is analogous. Therefore, the parameter errors obtained using the confidence-area method correspond to a confidence level > , and are most preferred from the point of view of reliability. Figure 12 presents a comparison of the projected confidence regions for the parameters r1 and x1 obtained using various methods. Note that in our previous paper [1] small rectangles show, as in Fig. 12 of the given paper,

the confidence regions obtained by the differentialcorrection method. In spite of the fact that the most probable confidence interval obtained using the 2 method exM ceeds the corresponding interval obtained using the 2 method, the 2 confidence interval contains the P P true value of the parameter in 99% of cases, and the 2 confidence interval in 85% of cases (for a oneM dimensional projection). This reflects the fact that the probability of falling into the confidence intervals obtained using the 2 method is determined largely, not M by the most probable value, but by the character of the distribution of these intervals, which are themselves random quantities.

6.3. "Internal" and "External" Errors
Our results clearly demonstrate and explain the qualitative difference between the magnitudes of parameter errors derived using the differential-correction and confidence-region methods. We constructed theoretical light curves using the upper and lower values for r1 , r2 , x1 , and x2 from Table 9 and fixing the orbital inclination of the binary to be i = 88.28 (Fig. 13; see also Fig. 12). The data in these theoretical light curves form an "error corridor", within which the observed values of the given realization of the light curve are located. The error corridor for the differential-correction method is fairly narrow, and some of the "observed" points in the given realization densely press against its walls.
ASTRONOMY REPORTS Vol. 53 No. 8 2009


PARAMETER ERRORS IN INVERSE PROBLEMS

747

The "observed" points in some realizations of the light curve may lie outside the indicated error corridor. The analogous error corridor for r1 , r2 , x1 , and x2 derived using the 2 method (Fig. 13) is appreciably M broader, and clearly allows for the statistical unpredictability of the light-curve realizations. We emphasize again that, in adopting errors found with the differential-correction or Monte Carlo methods, we must bear in mind that these include only "internal" errors, which may not overlap with the "internal" errors of the parameters found from an analysis of another realization of the light curve (for example, obtained at a different epoch). The probability that the true parameters all fall into their error intervals in the differential-correction method is substantially lower (by a factor of 1.3-1.4) than the probability for each error interval separately. The basis for calling the errors determined using the differential-correction and Monte Carlo methods "internal" is that the "internal" statistical distribution of the derived central values of the parameters are used when computing these errors. The calculation of the parameter errors in the confidence-region method uses the "external" statistical distribution of the observed light-curve values, so that the errors yielded by this method can be considered to be "external". 7. DISCUSSION Our computations have shown that trustworthy derivation of both limb-darkening coefficients based on the light curve for the eclipsing system YZ Cas is possible only if the accuracy of the individual measurements is 0 10-4 (on an intensity scale with Lfull = 1). It is not possible to determine the limbdarkening coefficient for the smaller star reliably when 0 10-3 . Not only the absolute accuracy of the individual brightness measurements for the system is important, but also the depth of the eclipse. For example, in the YZ Cas system, the depth of the secondary minimum (full eclipse) is 10% of the total brightness of the system, which is fully sufficient to determine the limb-darkening coefficient with 0 10-4 . However, if the depth of this eclipse were smaller, the required accuracy of the individual measurements would be higher. When fitting the primary minimum of the YZ Cas light curve (annular eclipse), we are able to reliable determine the linear limb-darkening coefficient x1 = 0.2917 ± 0.1199 for 0 0.0034. The light curves obtained using linear and quadratic limb-darkening laws were identical, indicating the smallness of non-linear limb-darkening effects. This makes it possible to fit the light curve reliably using a linear limb-darkening law when the accuracy of the light curve measurements is 10-3 .
ASTRONOMY REPORTS Vol. 53 No. 8 2009

Table 11 shows that the individual error intervals for x1 and x2 derived using the Monte Carlo and differential-correction methods contain the true values in 68% of cases, while the both parameters lie in their error intervals simultaneously in only 49% of cases. This reflects the non-reliability of the differential-correction and Monte Carlo methods in terms of the simultaneous incidence of the unknown parameter values in their error intervals. Recall that the confidence intervals for x1 and x2 in the 2 method are projections of the six-dimensional M region of the parameters r1 , r2 , i, x1 , x2 ,and I0 /I0 onto the x1 and x2 parameter axes. The probability that a parameter falls in the confidence region D corresponds fully to the specified confidence level = 68.2% (Table 11). Since the two-dimensional projection onto the (x1 , x2 ) plane forms an elongated "ellipse," the number of "hits" in this projected "ellipse" appreciably exceeds the confidence level = 68.2%. The situation for the 2 method is analogous. P Therefore, the parameter errors obtained using the confidence-region method correspond to > ,and are preferable in terms of their trustworthiness. For example, Table 1 shows that, in the case of a one-parameter model with = 0.90, with the number of observational points M = 10 - 10 000, the parameter tmax varies from 1.89 to 9.19; i.e., this is the factor by which the error of an unknown parameter is reduced when derived using the Monte Carlo or differential-correction methods. One of the important results of our analysis are the analytical relationships between the probable values of the confidence intervals obtained using the 2 and the errors obtained using the differentialM correction method for one-parameter and multiparameter problems [see (56) and (61)]. These show that tmax (and kP tmax ) decrease with growth in , while this parameter grows comparatively weakly with growth in M (Figs. 3­6). We again note that formulas (56) and (61) (for a linear dependence of the function describing the model on the unknown parameters) can be used to obtain the error interval for an unknown parameter for the 2 confidence-region method that is most M probable among those for adequate models from the corresponding error interval in the Monte Carlo or differential-correction method via multiplication by the coefficient tmax or kP tmax (see also Tables 1­7). After finding the error interval for the Monte Carlo or differential-correction method and verifying that the model is adequate in the 2 method (for which it M is sufficient to calculate the minimum of the corresponding residual), we can draw preliminary conclusions of a probabilistic nature about the size of the 2 confidence interval, without turning directly to M
(1) (2)


748

ABUBEKEROV et al.

the more labor-intensive determination of the sizes of the projections of the confidence region onto the parameter axes. Briefly summarizing the results, the differentialcorrection and Monte Carlo methods yield estimates only of the "internal" parameter errors. At the same time, the confidence-region method provides estimates of the "external" errors of unknown parameters (see in this connection the study of Popper [11]). Let us make one more important comment. The magnitude of the parameter errors obtained in the differential-correction, Monte Carlo, and 2 methods P do not depend on the sample volume (number of points M in the light curve), and are stable in the sense that, the parameter errors for different realizations of a random process (light curve) with other conditions the same will differ only comparatively weakly (the corresponding differences are much less than 100%). These "good" properties of these methods are associated with the fact that they assume a priori that the model used is fully correct (the possibility of an error of the second kind is taken to be negligible). In the case of the 2 and Fisher FM,N -M methM ods, there is no such a priori assumption that the model is fully correct, and the adequacy of the model is verified using the observtional data in the process of constructing the confidence regions. Therefore, in the 2 and FM,N -M methods, the resulting parameter M errors depend on the sample volume M and vary strongly from one realization of the random process to another (for some realizations, the errors can even degenerate into an empty set). The dependence of the parameter errors in the 2 and FM,N -M methods M on the sample volume M is very weak. For example, as follows from Table 1, when the number of points M is varied by a factor of two, from 100 to 200, the most probable value of relative error in the 2 M method varies for = 0.7 by a factor of 17%; when the sample volume is varied by three orders of magnitude (from M = 10 to M = 10 000), other conditions being equal, the most probable value of relative error varies by only a factor of 5. On the other hand, it is possible to obtain reliable estimates of parameter errors when a model is rejected at a fairly high significance level (corresponding to a reduced chi-squared that is close to unity) in the 2 , FM,N -M , and 2 M P methods. 8. CONCLUSION We have shown that determining the errors of parameters in the differential-correction and Monte Carlo methods leads to a substantial reduction in the magnitude of the errors. These methods yield only

the "internal" errors of the parameters for a given realization of an observed light curve. The parameter values for different realizations of a light curve can even not have common values in the error intervals derived using the differential-correction and Monte Carlo methods. The comparatively small parameter errors obtained in the differential-correction and Monte Carlo methods is due to the application of the a priori assumption that the model is fully correct, and also the application of a simple, normally distributed statistic; at the same time, the solution of the inverse problem is searched for using another statistic--2 , which M is only derived from a normal distribution. It is much more natural to search for the parameter errors using the same statistic that is applied to find the central parameter values, for example, the 2 statistic. It M is desirable not to apply the artificial assumption that the model considered is fully correct; i.e., not to exclude the possibility of an error of the second kind. This requirement is satisfied by the 2 statistic [12]. M Therefore, correct calculation of parameter errors requires the use of the confidence-region method with the 2 statistic, which makes it possible to judge the M "external" errors of model parameters, as well as as the adequacy of the model itself. Given the labor-intensive nature of computating the sizes of projections of a confidence region onto the parameter axes, we stress the importance of our proposed new method for estimating parameter errors (Tables 1­7, Figs. 4­6). Namely, in the case of a linear dependence of the model function on the unknown parameters, formulas (56) and (61) (see also Tables 1­7) can be used to make preliminary conclusions of a probabilistic nature about the sizes of the 2 confidence intervals based on the corresponding M intervals for the Monte Carlo or differential-correction methods via multiplication by tmax or kP tmax , without the need for the more labor-intensive determination of the projections of the confidence region onto the parameter axes. If the confidence region is understood in an asymptotic sense, our new method for estimating parameter errors can also be applied to non-linear problems. ACKNOWLEDGMENTS The authors thank Professor A.G. Yagola for valuable comments and discussions. This work was supported by the Russian Foundation for Basic Research (project code 08-02-01220), the Program of State Support for Leading Scientific Schools of the Russian Federation (NSh-1685.2008.2), a President of the
ASTRONOMY REPORTS Vol. 53 No. 8 2009


PARAMETER ERRORS IN INVERSE PROBLEMS

749

Russian Federation Grant to Support Young Russian PhDs in Science and Their Supervisors (MK2059.2007.2), and the Analytical, Departmental, Directed Program "The Development of the Scientific Potential of Higher Education" (RNP-2.1.1.5940). REFERENCES
1. M. K. Abubekerov, N. Yu. Gostev, and A. M. Cherepashchuk, Astron. Zh. 85, 121 (2008) [Astron. Rep. 52, 127 (2008)]. 2. B. C. Carlson, arXiv:math/9409227v1 [math.CA] (1994). 3. A. V. Goncharskii, A. M. Cherepashchuk, and A. G. Yagola, Ill-Posed Problems in Astrophysics (Nauka, Moscow, 1985), p. 53 [in Russian]. 4. B. M. Shchigolev, Mathematical Processing of Observations (Fizmatgiz, Moscow, 1962) [in Russian].

5. S. Wilks, Mathematical Statistics (Wiley, New York, 1962; Nauka, Moscow, 1967). 6. G. E. Kron, Astrophys. J. 96, 173 (1942). 7. T. M. Brown, D. Charbonneau, R. L. Gilliland, et al., Astrophys. J. 552, 699 (2001). 8. A. V. Goncharskii, S. Yu. Romanov, and A. M. Cherepashchuk, Finite Parametric Inverse Problems in Astrophysics (Mosc. Gos. Univ., Moscow, 1991) [in Russian]. 9. J. Southworth, Mon. Not. R. Astron. Soc. 386, 1644 (2008). 10. W. Van Hamme, Astron. J. 106, 2096 (1993). 11. D. M. Popper, Astron. J. 89, 132 (1984). 12. A. M. Cherepashchuk, Astron. Zh. 70, 1157 (1993) [Astron. Rep. 37, 585 (1993)].

Translated by D. Gabuzda

ASTRONOMY REPORTS

Vol. 53 No. 8 2009