{'fraction_non_alphanumeric': 0.03712751243870458, 'fraction_numerical': 0.018842167657571893, 'mean_word_length': 4.34299424184261, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mathematics of Statistical Parallax and the Local Distance Scale

We present a mathematical analysis of the statistical parallax method. The method yields physical insight into the maximum-likelihood determinations of the luminosity and velocity distribution and enables us to conduct a vigorous Monte Carlo investigation into various systematic effects. We apply our analytic formalism to the RR Lyrae sample of Layden et al. The velocity distribution of RR Lyrae stars is highly non-Gaussian, with kurtoses K_\pi = 2.04, K_\theta = 3.22 and K_z = 4.28 in the three principal directions, but this has almost no effect on either the best fit or the uncertainty of the luminosity determination. Indeed, our principal result is that the statistical parallax method is extremely robust in the face of all systematic effects that we considered. Our analysis, applied to the Layden et al. RR Lyrae sample, strictly confirms the majority of their results. The mean RR Lyrae absolute magnitude is M_V = 0.75 +/- 0.13 at the mean metallicity of the sample<[Fe/H]>= -1.61, compared to M_V = 0.71 +/- 0.12 obtained by Layden et al. Most of the difference is due to Malmquist bias which was not considered in previous studies. We also analyze a semi-independent non-kinematically selected sample of stars with metallicities at the [Fe/H]<= 1.5 taken from Layden et al. and Beers&Sommer-Larsen and obtain M_V = 0.79 +/- 0.12 at<[Fe/H]>= -1.79. Additionally, this analysis yields measurements of the radial bulk motion (4 +/- 10 km/s) and vertical bulk motion (0 +/- 6 km/s) of the halo relative to the Local Standard of Rest.


Introduction
The discrepancy between the Cepheid and RR Lyrae distance estimates to the Large Magellanic Cloud (LMC) is a long-standing problem (e.g., van den Bergh 1995). There are two possible sources of disagreement between the Cepheid and RR Lyrae distance scales: either the compared populations of Cepheids or RR Lyrae stars in the Milky Way and LMC are different or there are substantial errors in the Cepheid or RR Lyrae calibration. In this paper we investigate the statistical and systematic errors associated with the statistical parallax method, which is one of the main calibration methods for RR Lyrae stars. We prove that the statistical parallax method is very robust. Consequently, miscalibration of the absolute magnitude of RR Lyrae stars is not likely to be the solution to the distance scale problem and it is fair to state that RR Lyrae stars constitute reliable distance indicators. We conclude that either the Cepheid distance scale is incorrect or the currently available data do not allow for a reliable comparison between RR Lyrae stars in the Milky Way and LMC.
Historically, the methods of secular parallax and classical statistical parallax (e.g., Trumpler & Weaver 1962) have been used to expand the method of trigonometric parallax to more distant objects by increasing the baseline of the measurement. Secular parallax takes advantage of the fact that the Sun moves at a speed W ∼ 20 km s −1 relative to disk populations and W ∼ 200 km s −1 relative to halo populations, corresponding to 4 AU per year and 40 AU per year, respectively. After a time t of several decades, the accumulated baseline W t is several orders of magnitude larger than the Earth-Sun baseline used in trigonometric parallax. If the speed of the individual star (rather than of the population as a whole) were known to be W , then one could measure the distance with a precision better by a factor W t/AU than using trigonometric parallax for the same astrometric accuracy. The drawback is that the stars have a dispersion σ, so that even in the limit of perfect measurements, the precision of the distance determination is inversely proportional to the "Mach number", κ = W/σ, which is typically of order 1. By applying this method to a sample of N stars, one can beat down the noise and achieve a precision ∝ N −1/2 κ −1 .
The velocity dispersion itself can be used to obtain an independent measurement of the distance using the method of classical statistical parallax. That is, the distance scale of a class of stars can be estimated by forcing their radial velocities and proper motions to reproduce the same velocity dispersion. The precision of this method is ∝ N −1/2 . The two determinations can be combined to achieve a joint precision ∝ [N(1 + gκ 2 )] −1/2 , where g ≃ 1/6 is a geometrical factor that we derive below. Nevertheless, because the method is limited by the number of stars rather than the precision of the measurements, its usefulness has been restricted to ever more distant classes of stars as astrometry has improved.
Only recently (e.g., Murray 1983;Hawley et al. 1986;Strugnell, Reid, & Murray 1986) was it fully realized that secular parallax and classical statistical parallax are actually two aspects of the same generalized method, now simply called "statistical parallax". Secular parallax is based on forcing equality between the three first moments of the velocity distribution (the bulk motion W) as determined from radial velocity and proper motion measurements, while classical statistical parallax is based on forcing equality of the six second moments (the six independent components of the velocity covariance matrix C ij ). In the modern combined version of statistical parallax one simply determines ten parameters simultaneously by applying maximum likelihood. The ten parameters are an overall distance scaling factor η (relative to an initial arbitrary distance scale) plus the nine first and second moments, W and C ij . Hawley et al. (1986) and Strugnell, Reid, & Murray (1986) were the first ones to use a maximum likelihood analysis to determine these 10 parameters. Actually both studies use 11 parameter models of the form introduced by Murray (1983), but hold fixed the 11th parameter, the dispersion in η. As we show in §3.1, current data do not allow one to put any useful constraint on this parameter, and we remove it from our primary analysis. The dispersion in η also gives rise to Malmquist (1920) bias. However, in the analysis conducted by Hawley et al. (1986) and Strugnell et al. (1986), the inclusion of the dispersion in η does not by itself correct for Malmquist bias. Correcting for the Malmquist bias requires additional modeling of the selection criteria of the stars in the sample and of the distribution of absolute magnitudes. It is most effectively done on a star by star basis (for a discussion see e.g. Smith 1987, also Ratnatunga & Upgren 1997). However when the detailed selection criteria are not available and even the magnitude of the scatter is uncertain, such procedures are no better than correcting the final result for the "average" Malmquist bias (e.g. Jung 1970). 2 Ideally, the statistical parallax method should be applied to a group of stars that are: 1. dynamically homogeneous, i.e. all the stars are drawn from a single velocity distribution (not constrained to be Gaussian) regardless of their locations 2. standard candles, i.e. all have the same absolute magnitude.
Theoretically, condition (1) can be met by careful selection of stars in the nearby solar neighborhood. In practice, to obtain a statistically satisfactory sample often requires that stars be collected from a region comparable (∼ 3 − 5 smaller) in size to the distance to Galactic center. In such a case one does not expect condition (1) to be exactly satisfied.
However, one can reasonably assume that in this part of the Galaxy the bulk motion and velocity ellipsoid change monotonically with distance from Galactic center. Hence, it is 2 Actually, for small scatter in absolute magnitude, the two approaches are equivalent because they differ only by negligible second order terms. important to probe parts of the sky towards and away from the Galactic center. Although it is difficult to check whether condition (1) is satisfied, it is possible to investigate some of the possible systematic errors analytically and to quantify others through Monte Carlo simulations. Condition (2) must be relaxed somewhat, and it suffices to have stars with a small scatter around calibrated absolute magnitudes, as we discussed above and as we further discuss in §3.
The RR Lyrae stars of the Galactic spheroid constitute a very good sample for analysis by statistical parallax. They are close to satisfying condition (1), being abundant in the solar vicinity and distributed reasonably evenly in the sky. Their absolute magnitudes can be obtained from the absolute magnitude -metallicity relation (e.g., Carney et al. 1992) with a scatter not likely to exceed 0.15 magnitudes (condition 2). Additionally, RR Lyrae stars can be observed not only in the Milky Way, but also in neighboring galaxies like the LMC.
The organization of this paper is as follows. In §2 we analyze a simple model of a stellar system with an isotropic velocity ellipsoid and we estimate the error in the distance scaling parameter η. We also discuss the role of observational errors and the systematic error caused by the mis-estimation of these observational errors. In §3, we apply the results of §2 to investigate various possible systematic effects that may bias the statistical parallax solution. In §4, we analyze the full maximum likelihood formulation of the problem for the limiting case of negligible measurement errors. This case mimics the actual situation because the intrinsic velocity dispersion, being much larger than the observational errors, dominates the statistical uncertainty. We obtain algebraic expressions for the uncertainty in all 10 parameters, η, W i , and C ij . In particular, the error in η for the general case has the same form as it does for the naive model of §2. That is, as we anticipate in §2, the seemingly complicated "black box" of maximum likelihood can be understood in simple physical terms. In §5 we conduct the complete maximum likelihood analysis of the general case including observational errors in stellar velocities. We then reanalyze the Layden et al. (1996) sample of 162 halo stars found in the solar neighborhood and confirm most of their findings. In §6 we conduct Monte Carlo simulations to confirm our analytic results and find possible biases, e.g., those induced by the anisotropic spatial distribution of the stars in the sample. In §7 we analyze a sample of halo stars constructed by combining the samples of Layden et al. (1996) and Beers & Sommer-Larsen (1995). Finally, in §8 we summarize our results and discuss their implications for the local distance scale. In particular, we conclude that the statistical parallax method is very robust against various systematic effects.

Isotropic velocity dispersion
In this section we conduct the heuristic analysis of the statistical parallax method. It is aimed at explaining the physical foundations of this powerful method and gives the basis for the analysis of possible systematic errors. In §4 and §5, we develop a more rigorous mathematical treatment which is then used in our Monte Carlo simulations described in §6. As mentioned in §1, in the general case there are 10 parameters of the fit: the distance scaling parameter η, 3 components of bulk motion W i , and 6 independent components of symmetric velocity covariance matrix C ij . We initially consider a stellar system with an isotropic velocity ellipsoid. The number of independent components of C ij is thereby reduced from 6 to 1: C ij = σ 2 δ ij .
We define the scaling parameter η by where L true is the actual luminosity of the chosen class of star (e.g., RR Lyrae), and L assumed is its assumed luminosity, which can be chosen arbitrarily for purposes of making the calculation.
Below we estimate analytically the error associated with the determination of the distance scaling parameter η. When dealing with errors we shall apply the general rule that combining errors of statistically independent quantities we add their variances to get the total variance, computing the resultant error of two statistically independent measurements of the same quantity, we add variances harmonically.
We decompose the velocity of each star into its radial and tangential components, u = (u r , u t ). Now we can conveniently separate two sources of information about η. The first comes from forcing equality between the velocity dispersions as determined from radial and transverse velocities ("classical statistical parallax"). The second comes from forcing equality in the bulk motion inferred from radial and transverse velocities ("secular parallax"). We assume that these two estimates of η are independent in the statistical sense. To simplify the analysis we initially assume that there is no bulk motion, which implies that u = 0, where the brackets, here and afterwards, symbolize averaging unless otherwise noted. Thus where u t = |u t |. The inferred transverse speedũ t is related to the proper motion µ of a star where d true and d assumed are the true distance and the distance inferred based on the assumed luminosity, respectively.
We combine equations (2) and (3) to express η as Each component of the velocity, u i , of a given star may be regarded as being drawn from a one-dimensional distribution f (u i ). The fractional error in the velocity dispersion based on a sample of N stars is where in the last step we have used the definition of the kurtosis K for u i = 0: By equation (5) ∆(u 2 r ) u 2 Recall thatũ 2 t in equation (4) designates the sum of the squares of the transverse velocities in two perpendicular directions. Thus Adding the fractional errors from formulae (7) and (8) in quadrature gives the fractional error in η (see eq. (4)), where we have used ∆η/η = (1/2)∆(η 2 )/η 2 , and where in the last step we have evaluated the expression for the case of a Gaussian distribution.
The bulk motion of the whole sample relative to the Sun gives independent information about η. For simplicity we assumed in the above analysis that the bulk motion was zero, but now we relax this assumption. In analogy to equation (4) we can write where the index i indicates the component and the indices r and t mean "as inferred from radial and transverse velocities", respectively. In general, we decompose the bulk motion into components in an arbitrary frame of reference which generally results in all three components having non-zero values. Here for the sake of simplicity, we pick the frame of reference with z-axis aligned with the bulk motion of the sample. We restrict attention to the radial velocities of the stars but the results are representative of each of the components. One can express the radial velocity as the bulk motion along the line of sight with uncertainty equal to the velocity dispersion σ: where W 3,r represents the bulk motion component inferred from the radial velocity. We may construct an estimator of any of the bulk motion components (in this simplified case of W 3,r ) by dividing equation (11) (or its generalized form) by the appropriate angular dependence: For the sample of N stars randomly distributed in the sky, the fractional error in estimating where in the last step we have used the fact that cos 2 θ = 1/3 averaged over all the angles.
The transverse velocities contain twice as much information as the radial ones, so Applying equation (10) and adding the errors from equations (13) and (14) in quadrature, we find that A similar analysis can be made to extract information from W 1 and W 2 . Since these three pieces of information about η are independent, we add variances harmonically and obtain where κ is the "Mach number", the ratio of bulk motion to dispersion. Combining equations (9) and (16) yields The typical values for population II stars in the solar vicinity are W ∼ 200 km s −1 and σ ∼ 100 km s −1 , κ ∼ 2 so (2/9)κ 2 ∼ 8/9. Hence, 60% of the information about the distance scale comes from velocity dispersions and 40% from the bulk motion.

Measurement Errors
Here we investigate how measurement errors affect the statistical accuracy. We adopt the simplified model discussed above with isotropic dispersion σ, bulk motion W , typical measurement errors σ r and σ µ , and typical stellar distances D. We define δ r ≡ σ r /σ and δ µ ≡ Dσ µ /σ. If η is determined solely from the bulk motion, then the accuracy is On the other hand if η is determined from the dispersion, (∆η/η) 2 = 0.5N −1 ([(1 + δ 2 r ) + 0.5(1 + δ 2 µ )]. Thus, the accuracy of the overall measurement is degraded by a fractional amount, where we have adopted δ r ∼ 0.2 and δ µ ∼ 0.4 appropriate for RR Lyrae stars. Hence, realistic measurement errors have an extremely small effect on the precision of the luminosity determination, a result which we confirm numerically in §6. This justifies our approach of ignoring measurement errors in our analytic investigation of various other effects conducted in §3.
Most generally, the uncertainty in the statistical parallax determination is set by: 1. the intrinsic velocity dispersion of the system 2. the bulk motion of the sample 3. the observational errors 4. the size of the sample Once the stellar system is given we have no control over items 1 & 2, but we can still vary items 3 & 4 to minimize the uncertainty. In the specific case of negligible observational errors (meaning small relative to intrinsic velocity dispersion as in the considered case), the only remaining factor influencing the uncertainty is the size of the sample. That is, the only way to substantially reduce the uncertainty is to increase the sample size and not decrease the observational errors.

Mis-estimate of the Observational Errors
As shown in §2.1, observational errors increase only slightly the uncertainty in the estimate of η. Observational errors can also result in a biased estimate of η if they are not properly accounted for. Suppose that the radial velocity errors are σ r , but that they are taken to be zero in the analysis. (For simplicity, we initially assume that there are no proper-motion errors.) If η is estimated by comparing the radial and proper-motion dispersions, the radial dispersion will be overestimated by a factor 1 + (σ r /σ) 2 while the proper-motion dispersion will be properly estimated. Hence η will be overestimated by a fraction (δη/η) = [1 + (σ r /σ) 2 ] 1/2 − 1 ∼ (σ r /σ) 2 /2. On the other hand, there will be no effect on the estimate of η based on comparing the radial-velocity and proper-motion bulk motion.
Hence the total systematic error in the combined determination will be the weighted average of the two: (δη/η) ∼ (σ r /σ) 2 /[2 + κ 2 /3] ∼ 0.3(σ r /σ) 2 . If the true observational error is σ r , but the analysis incorrectly treats the error asσ r , then (δη/η) ∼ 0.3δσ 2 r /σ 2 , where δσ 2 r ≡ σ 2 r −σ 2 r . Finally, we allow for a mis-estimate of the proper-motion errors and define δσ 2 µ ≡ σ 2 µ −σ 2 µ in analogy to the radial-velocity term. We take the typical distance of stars in the sample to be D and estimate a net systematic error of Since σ r /σ ∼ 0.2 and Dσ µ /σ ∼ 0.4, the most likely source of a major effect is mis-estimation of the proper motion errors, but even this is not likely to be large. For most RR Lyrae stars in the sample, the proper-motion error is estimated to be 5 mas yr −1 . Suppose that the true value is 20% higher. Then η would have been underestimated by only (δη/η) ∼ 1.7%.

Dispersion of the Luminosity
Throughout our treatment, we have assumed and will assume in §4 and §5 that the entire population of stars has exactly the same luminosity. This is not customary. One can in principle fit for the dispersion in luminosity (in which case there are 11 fit parameters instead of 10) or one can assume a certain dispersion which then affects the values of the remaining 10 parameters (e.g. Hawley et al. 1986;Strugnell et al. 1986;Layden et al. 1996). Here we show that for the RR Lyrae sample, there are several orders of magnitude too little information to determine the dispersion from the data as was also suggested by previous studies. Moreover, we show that the effect of including the dispersion is one order of magnitude smaller than the statistical errors, and that one is therefore justified in accounting for this effect separately.
We assume that the stars have a range of luminosities which we parameterize as a dispersion in absolute magnitude M V , σ M , but that the analysis is conducted assuming that all stars have the same luminosity. We work within the simplified framework of §2. The inferred distance to each star in the best-fit solution will then deviate from the true distance by O(σ M ). This will increase both the variance in the inferred transverse speeds and the mean. Averaged over all directions, the fractional increase of the dispersion is 2(ln 10/5) 2 [1 + κ 2 /6]σ 2 M . Thus, if η were determined by matching the dispersions of the radial velocities and proper motions, it would be underestimated by On the other hand, if η were determined by matching the bulk motions, it would be underestimated by a fraction (δη/η) = −0.5(ln 10/5) 2 σ 2 M ∼ −0.1σ 2 M . Thus, by comparing the two estimates, one could in principle determine σ M . For the RR Lyrae sample, however, the precision of each method is only ∼ 7%, so that the precision of the difference is ∼ 10%. This means that the data set would have to be increased ∼ 400 fold in order to detect plausible values of σ M ≤ 0.15 mag. If the dispersion is simply ignored and the two methods are averaged (as they automatically are in the maximum likelihood approach), then η will be underestimated by (δη/η) = −(ln 10/5) 2 [12 + 3κ 2 ]σ 2 M /[12 + 2κ 2 ] ∼ −0.25σ 2 M . This is an extremely small correction for plausible values of σ M . The underestimate in η in turn induces changes of the bulk motion and dispersions, which we estimate numerically in §6.

Malmquist bias
Malmquist bias is another correction to η which scales ∝ σ 2 M , but it is of opposite sign. As we discuss in the introduction, the Malmquist bias correction is not automatically taken into account by simply considering the 11-parameter model, with the 11th parameter being the dispersion in η. The assumption intrinsic to most such models (and also to our 10-parameter model) is that the stars in the sample are representative of the whole population. This assumption is not correct for samples that have been selected based on stellar magnitudes. Malmquist bias must therefore be put by hand 3 . Malmquist (1920) showed that the mean distance of a magnitude-limited sample drawn from a homogeneous distribution with a Gaussian scatter in absolute magnitude will be higher than that of one without scatter by (δη/η) ∼ 3(ln 10/5) 2 σ 2 M ∼ (ln 10/5)1.38σ 2 M ∼ +0.64σ 2 M .
Malmquist bias is a particular form of "selection bias" which arises whenever a sample deviates systematically from the underlying population because of the way that it was selected. The magnitude of the selection bias depends, in general, on the form of the selection. The RR Lyrae sample is "local" in the sense that it contains stars that are less than 2 or 3 kpc from the Sun, but does not contain substantially more distant stars (even though such stars are plentiful). Some procedure (whose details are not known to us) had to be applied to select this sample. In principle, one could distinguish between nearby and distant stars by several methods, including apparent magnitude, trigonometric parallax, and proper motion. However, while we do not know exactly how the sample was selected, we do know that of these types of information only apparent magnitudes were available.
Therefore, each star can have been included only because it satisfied some magnitude limit, even if this magnitude limit is different from star to star. The correction calculated by Malmquist (1920) is independent of the magnitude limit, so the same correction applies to all sub-samples of fixed magnitude limit, even samples with only one object. Therefore, it applies to the sample as a whole. For an inhomogeneous underlying population, the numerical pre-factor in 3(ln 10/5) 2 σ 2 M may in general differ from 3. However, we have confirmed numerically that for the distribution of RR Lyrae stars from which the Layden et al. (1996) sample is drawn, the pre-factor is in fact very close to 3.

Velocity-Position Correlations
Statistical parallax makes the implicit assumption that the stars seen in all directions have the same velocity distribution. Since the RR Lyrae stars in the sample have typical distances ∼ 2 kpc, which is a significant fraction of R 0 ∼ 8 kpc, it is possible that this assumption is not valid. Suppose, for example, that the velocity dispersion in the z direction falls with distance from the Galactic plane. If one attempted to fix the luminosity by matching the z-dispersions of the radial velocities and the proper motions, one would underestimate its value. This is because the radial-velocity measurements would be made primarily on stars far from the plane, while the proper-motion measurements would be made primarily on stars near the plane. The high dispersion of the latter would be mistakenly interpreted as indicating that the stars were closer and hence less luminous than they actually are. While this systematic bias would be diluted by unbiased dispersion measurements in the other two directions and by unbiased measurements of the bulk motion, it would still affect the final result. (We focus attention on the z direction because that is the only axis for which there is a physical plane of symmetry. Any correlations along the other two axes would tend to have opposite signs in the positive and negative directions and hence would cancel to lowest order.) To make a quantitative estimate of the size of this effect, we consider a simple model where the stars are observed over a uniform sphere of radius D, and have a z velocity Consider an ensemble of stars at a distance r and angle θ from the north Galactic pole. These will enter into the radial-velocity and proper-motion estimates of the z dispersion with statistical weights cos 2 θ and sin 2 θ, respectively. If the correct luminosity of the stars were adopted, the respective dispersion estimates would be Thus, one would be led to make an underestimate (δη/η) = −(9/64)f . Since the effect is diluted by measurements of the dispersion in the other two directions and of the bulk 25, this would be a small but perhaps not completely negligible correction.
We have therefore reanalyzed the sample including an 11th parameter, the velocitydispersion gradient in the z direction. The best fit scale length for this gradient is 7 kpc, but since the reduction in χ 2 is only ∆χ 2 = 1.15, the gradient is not statistically significant. In any event, even if the gradient is real it introduces a systematic error of only (δη/η) = −0.005.

Analytic Predictions
In this section we use the maximum likelihood method to analyze the errors intrinsic to the statistical parallax method. We again neglect velocity measurement errors. We argued in §2.1 that for the halo stars of the Galaxy the accuracy of the determination of the model parameters scales with the size of the sample of stars and is barely influenced by observational errors. Having selected the group of stars (e.g., RR Lyrae stars), we choose an orthogonal frame of reference and resolve each star's velocity into components designated as u i . We assume that the velocity distribution is a three-dimensional Gaussian, but we do not assume isotropy. At first Gaussianity of the assumed velocity distribution seems to be a rather restrictive condition, especially in view of the fact that the actual distribution is highly non-Gaussian. However, we argue at the end of this section ( §4.1), and confirm numerically in §6 that no bias is introduced by assuming that the distribution is Gaussian.
The probability of finding a star with three velocity components in the ranges where W is the bulk motion relative to the Sun, C is the velocity covariance matrix, and |C| is its determinant. The volume element in three-dimensional true velocity phase space is d 3 u, whereas d 3 u| assumed is the corresponding element in "assumed" velocity phase space.
Hence the probability density of finding a star with the observed velocity components in the ranges (u i , u i + du i ) is given by: The logarithm of the total probability of finding the system in the observed state is: where L k is the probability (23) associated with the k th star. The curvature matrix is given by where α and β range over the 10 parameters of the model. The covariances among the parameters are then given by We designate the inverse of the covariance velocity matrix as which is a special case of the general definition (52). We express in terms of its radial and transverse components. Note that t i is the projection ofũ t , defined in §2, on i-th direction. We also define v i to be the random part of the velocity: We evaluate the curvature matrix in two steps. First we evaluate the contribution of each star to equation (24). Next, we sum over all stars under the assumption that they are isotropically distributed on the sky. We find where the symbol (k ↔ l) means "add the terms with k and l exchanged, but only if k = l", and where we have employed Dirac notation, i.e.
Averaging over positions leads to: The inverse of the matrix B, given by equations (37) -(42), describes the variances and covariances of the parameter estimates (see eq. (26)). Note that no observational errors have been included in this treatment and that the uncertainties in the model parameters scale with the square root of the size of the stellar sample. We now adopt the frame of reference for which the axes are aligned with the axes of the velocity ellipsoid i.e. where C ij is diagonal. We then find that the errors in η, W i , and C ij are given by Here κ is the generalization of the Mach number to an anisotropic distribution.
We have derived fully analytic estimates of the errors in all parameters. Note that the result given in equation (43) can be written as That is, the error in the distance estimate (eqs. (48) and (43)) for the general case has the same form as the error for an isotropic distribution derived using a few basic assumptions (eq. (17)).

Non-Gaussian Velocities
As we show in §6, the actual distribution of RR Lyrae stars is very far from Gaussian.
On the other hand, the formalism presented in §4 explicitly assumes that the velocities are Gaussian. That is, the likelihood given in equations (22) and (23)  If the underlying distribution is non-Gaussian, the Gaussian maximum likelihood procedure will return estimates of the errors that differ from the true errors. Consider for example the simple isotropic model with dispersion σ, kurtosis K, bulk motion W , and negligible observational errors. The maximum likelihood estimate of the error in η is (see eq. (43)) (∆η/η) −2 = (2/9)N[6 + κ 2 ] while the true error is (see eq. (17)) (∆η/η) −2 = (2/9)N[(12/[K − 1]) + κ 2 ]. For a more general velocity distribution that is the product of distributions in the π, θ, and z directions with dispersions (σ π , σ θ , σ z ) and kurtoses 4 To be more precise σ 2 will be an unbiased estimator of (N − 1)/N× the variance, but this slight difference has no practical impact on the discussion here.
(We confirm this estimate numerically and mention some practical complications in §6.) For an arbitrary velocity distribution, the true and maximum likelihood-estimated errors could in principle differ substantially. However, for the actual RR Lyrae population, 4/(K π − 1) + 4/(K θ − 1) + 4/(K z − 1) ≈ 7, close to the Gaussian value of 6. This implies that the maximum likelihood estimated errors are nearly equal to the true errors.

Complete analysis
In this section we obtain the formulae needed to analyze real data and to carry out Monte Carlo simulations. Our analysis here is similar to the one conducted in §4, but we account for a few additional effects.
The probability density of finding a star with velocity components in the ranges (u i , u i + du i ) is now given by: where w i is the bulk motion (defined more precisely below), and s i is the stellar velocity expressed in its local Galactic coordinate frame (the velocity that would be measured by an observer located at the star's position and at rest with respect to the Galactic center) with v ⊙ = (−9, 232, 7) km s −1 .
-25 - The matrix M is the full velocity covariance matrix, defined as Here σ r is the observational error in the radial velocity, σ t is observational error in each of the inferred transverse velocity components, and P ij is the projection operator on the plane of the sky. In equation (50) we explicitly subtract the velocity of the Sun v ⊙ from the star's velocity. The velocity v ⊙ is the sum of the velocity of the Sun's LSR v LSR and the peculiar motion of the Sun relative to its LSR v ⊙,LSR . The completely new element not present in formula (23) (nor in any previous statistical parallax analyses of which we are aware) is the rotation operator R ij that accounts for the fact that the star's velocity is drawn from a certain distribution (e.g., three-dimensional Gaussian) in the Galactic frame of reference centered at the star which is in general a very different distribution when expressed in terms of the Galactic frame centered at the Sun. The matrix elements of the rotation operator are functions of the star's position and the distance from the Sun to Galactic center. We shall assume that the distance to Galactic center is ∼ 8 kpc. The bulk motion, w i , is relative to the local Galactic frame. That is, in the Galactic frame at the Sun's position, the bulk motion is j (R −1 ) ij w j . As in the no-error case, the logarithm of the total probability of finding the stellar system in the observed state is the product of probability functions of the form (49) for all the stars: ln L = N k=1 ln L k .
We generalize equation (27) to We maintain the decomposition (28), u i ≡ u r r i + ηt i , and we introduce the short notation As previously, we use Dirac notation (36) to designate scalar products and the symbol (k ↔ l) means "add the terms with k and l exchanged, but only if k = l".
The first derivatives of ln L are: Second differentiation leads to: In deriving these results we have made use of the identities We apply equations (54) -(62) and Newton's method to find the maximum likelihood solution for the Layden et al. (1996) sample of halo RR Lyrae's. Table 1 gives the basic results and illustrates the process of correcting them for detected biases. Our "assumed" distance scale is the one derived by Layden et al. (1996). As a result, if our determination were in perfect agreement with Layden et al. (1996), the scaling parameter η should be exactly equal to 1. The four rows give the maximum likelihood values of the ten parameters characterizing the sample under four different sets of assumptions. None of the rows corresponds exactly to the results of Layden et al. (1996), but the first row with η increased by ∼ 0.015 does so approximately. The first row is not corrected for rotation or biases.
The values presented in the second row are obtained using a non-unit rotation operator (see discussion following eq. (51)). Including rotation increases η by about 0.4%, only slightly influences the bulk motion, but makes the velocity distribution less triaxial (C 1 2 22 closer to C 1 2 33 ) and more elongated (higher C 1 2 11 ). The third row gives our best estimate for the values of the parameters after correcting for the biases discussed in §6. We also take into account that maximum likelihood underestimates the velocity variances by a factor The stars that enter our analysis are those defined as "Halo-3" population from Table   3 in Layden et al. (1996). Also, we use Layden's best estimate of the proper motion errors (6.5 mas yr −1 ) for stars taken from the catalog compiled by Wan, Mao & Ji (1980) rather than the value (5 mas yr −1 ) artificially adopted by Layden et al. (1996) to match their own sample (although this makes almost no difference). We assume that RR Lyrae stars follow We checked that the solutions for different slopes are statistically indistinguishable from one another. By taking account of the metallicity, we restrict the possible scatter in RR Lyrae absolute magnitudes to the intrinsic scatter at fixed metallicity. To make the corrections discussed in §3, we adopt for the Layden et al. sample of field Galactic RR Lyrae stars σ M = 0.15. To make this estimate, we first inspected color-magnitude diagrams of several globular clusters and found that for these very homogeneous populations the dispersion is typically σ M = 0.08. We take this as a lower limit. The only available inhomogeneous sample at approximately fixed distance is the 28 RRab Lyrae stars of Hazen & Nemec (1992), which are 4 • east of the LMC center and for which σ M ∼ 0.17. We take this as an upper limit because some of the scatter may be due to the dispersion in distance. For example, if the LMC RR Lyrae population has an r −3.5 profile and the flattening c/a = 0.6 (similar to the Galactic population), we find a dispersion due to distance σ M,dist = 0.09, implying an intrinsic dispersion of σ M ∼ 0.14. In any event, those prefering other values of σ M should note that it is straightforward to find the corrected values of all parameters simply by scaling the difference between rows 3 and 4. The correction we apply in Table 1 is where the first term, estimated numerically, is due to the scatter in the absolute magnitude of RR Lyrae stars ( §3.1) and the second term is due to Malmquist bias ( §3.2).
To make connection to previous studies of Galactic structure, note that it is customary to use the following notation: Notice, however, that by using the non-unit rotational operator R we are actually measuring the underlying velocity distribution in the local Galactic frames of the stars in the sample under the assumption of Galactic axisymmetry. We therefore use (π, θ, z) rather than (U, V, W ). The rotation operator R is the appropriate first order correction to the rectilinear solution (e.g., Layden et al. 1996) regardless of whether the velocity distribution in the Galaxy is exactly axisymmetric. Our best estimate of the RR Lyrae absolute magnitude at the mean metallicity of the sample [Fe/H] = −1.61 is M V = 0.75 ± 0.13. The velocity ellipsoid is (σ π , σ θ , σ z ) = (171 ± 10, 98 ± 8, 96 ± 8) km s −1 and the RR Lyrae population is moving in θ direction at −211 ± 12 km s −1 relative to the Sun, where the errors are adopted from YYTA case in the Table 3.

Monte Carlo simulations
Here we present the results of Monte Carlo simulations aimed at checking and fine-tuning our analytic results. There are two main classes of simulations. One class keeps all the star positions from the Layden et al. (1996) sample unchanged giving results closely related to those obtained using the real sample. For the second class, the stars are placed randomly over the celestial sphere. The position-related biases can be determined by comparing the two classes of simulations.
In each simulation we construct 4000 mock samples of 162 halo stars (162 is the size of the Layden et al. 1996 sample). For each sample we generate a set of 162 space velocities drawn from a distribution with specified means, dispersions and kurtoses in the three principal directions. For each star we transform its velocity components from the star's Galactic frame of reference to the Sun's frame and, in some cases, add Gaussian measurement errors in accordance with the values given by Layden et al. (1996).
In the second step we find the most probable parameters describing each of the samples.
To analyze our mock samples we use exactly the same maximum likelihood procedure that was used to obtain the results for the real stars.
We perform several tests to check for the various systematic effects. The results of these investigations are summarized in the Tables 2 and 3. Table 2 gives the biases in η, w i , and C 1 2 ii found in each simulation. The first column says whether the observational errors were included or were set equal to zero. The second column tells whether the velocities were rotationally adjusted (operator R in equation (50)). The third column gives the information on whether the positions of the stars in the sample were the true ones from Layden et al. (1996) or were chosen randomly. The fourth column lists the assumed values of kurtoses in all three directions 5 . The fifth column assigns the name to each case. The name contains the most important information about the case. For example YYTH means that the observational errors were included in the analysis, that velocities were rotationally adjusted, that the true positions of the stars were considered, and that the kurtoses in all directions were higher than Gaussian (here equal to 4). Generally, "T" stands for true, "R" for random, "G" for Gaussian, "H" for high, "L" for low and "A" for actual. The next seven columns give the results for the most interesting parameters of the fit based on 4000 realizations. In all cases, the biases in the off-diagonal elements of the velocity covariance matrix (normalized in the same way as in the Table 1)  The first five rows test the correctness of our implementation of the maximum likelihood method. They contain the cases with Gaussian velocity distributions (K π = K θ = K z = 3) for which we expect the biases to be small. Indeed, the deviations from the input values in the cases with the random positions of the stars (NNRG and YYRG) are extremely small, that is, the obtained solutions are statistically indistinguishable from the input. For the cases with the true positions of the stars (NNTG, YNTG and YYTG) the deviations remain small but now they are statistically significant and should be treated as biases rather than statistical fluctuations.
The next three rows, below the blank line, give the results for the samples with velocities drawn from high-or low-kurtosis, non-Gaussian distributions. One may suspect that our maximum likelihood procedure derived assuming a Gaussian velocity ellipsoid will not be very accurate in such cases. To the contrary, the determination of the parameters is extremely robust and insensitive to the underlying velocity distribution. The case YYTA is constructed to reproduce the input kurtoses of the real sample of the Layden et al. (1996) stars. A glance at Figure 1 shows that the distributions of the stars in the sample are highly non-Gaussian. The measured kurtoses of RR Lyrae stars in three principal directions are: (2.18, 3.14, 3.93). The right panel of Figure 1 shows how close the distribution in the radial direction is to a box distribution. Analyzing our simulations we find that the most probable values of the kurtoses of the true underlying distribution of RR Lyrae stars in the solar neighborhood are (K π , K θ , K z ) = (2.04, 3.22, 4.28). Fortunately, as we see from the lower part of the Table 2, the parameter determination is almost completely independent of the form of velocity distribution.  operator (YYTG) affects the errors only slightly. Furthermore, the same is true for the high-kurtosis (YYTH) or low-kurtosis (YYTL) cases. The errors in the non-Gaussian cases cannot be predicted exactly from the formulae given in §4.1 because the distribution is not isotropic. Nevertheless, the predicted trend is clearly confirmed. We conclude that the analytic estimates of the uncertainties constitute excellent approximations to the realistic case. Both observational errors and higher kurtosis slightly increase the uncertainties.
In Table 4, we compare the kurtoses measured in the three principal directions to those of the underlying velocity distributions. First note that unless the underlying kurtosis is very low (as in the NYRL case), the measured kurtosis is underestimated due to the finite size of the sample (e.g., NNRG and NYRH cases). Additionally, measurement errors tend to Gaussianize the distribution, and as a result reduce kurtoses in high-kurtosis cases (e.g., YYRH vs. NYRH), increase kurtoses in low-kurtosis cases (e.g., YYRL vs. NYRL), and have little effect on a Gaussian distribution (e.g., YYRG vs. NNRG). The measured kurtoses are also biased due to the non-isotropic positions of the stars in the sample on the sky (e.g., NNTG vs. NNRG or YYTH vs. YYRH), but this effect is relatively less important for non-Gaussian distributions. Our best estimate for the underlying kurtoses of the halo RR Lyrae stars, based on the simulation results presented in Table 4, is (2.04, 3.22, 4.28) (YYTA).

Kinematics of Metal-Poor Halo Stars
As we emphasized in §2.1, the factor that fundamentally limits the precision is the size of the sample. This limitation is severe because to expand the sample requires first identifying the new RR Lyrae stars and then measuring their proper motions. Some additional stars could undoubtedly be found in the South with distances similar to those characteristic of the sample. In most cases, however, new proper motion studies would be required. Moreover, the increase in the sample size would be modest. To really increase the size of the sample substantially would require finding stars at greater distances. In the North, at least, proper motions would be available from the Lick Proper Motion study.
However, the size of the errors for these stars would begin to approach those of the proper motions themselves, increasing the possibility of systematic errors. Thus, this route would appear to require substantial additional work, lasting perhaps several decades. In addition, stars that are found too far from the Sun may have different kinematics which would introduce additional systematic errors.
There is, however, a poor man's route to increased statistics which has the side-benefit of yielding new kinematic information about the stellar halo. If one repeats the analysis of §2, but with different numbers of stars with radial velocity (N r ) and proper motion (N t ) measurements, one finds in place of equation (17). If one measured only the radial velocities for a substantial new population of RR Lyrae stars, one could hope to drive down the error by up to a factor 3 1/2 without obtaining any new proper motions. In fact, it is not even necessary that they be RR Lyrae stars. The only requirement is that the radial-velocity stars and the proper-motion stars have the same kinematics. This is trivially satisfied if all the stars are RR Lyrae, but it can also be attained through careful selection of other stars.
For this purpose, we turn to the catalog of 1936 non-kinematically selected metal-poor ([Fe/H]< −0.6) stars of Beers & Sommer-Larsen (BSL, 1995 Layden et al. (1996) sample are repeated from the final rows of Tables 1 and 3, and the corresponding data are presented for the non-kinematic sample. The two results appear quite consistent, but since they are based partly on the same data, a careful statistical analysis is required.
The first point to note is that Table 5 gives only the statistical error, i.e., the error that arises primarily from the finite size of the sample and secondarily from measurement errors. For the non-kinematic sample, there is an additional error arising from the possible differences in the fraction of thick-disk stars in the two sub-samples. Let these fractions be are determined almost entirely from the BSL sample, the thick disk stars will generate an additional error δη/η ∼ gf ∆ where (1 − g) ∼ 0.5 is the ratio of the typical thick-disk to typical halo speeds relative to the Sun. We estimate the thick disk fraction f from the fact that the non-kinematic sample has a rotation speed of 35 km s −1 (see Table 5). Since halo stars have approximately zero rotation (BSL), this net rotation must be produced by thick-disk contamination. Since the rotation speed of the thick disk is of order 170 km s −1 (Casertano, Ratnatunga, & Bahcall 1990), the contamination level is f ∼ 35/170 ∼ 0.2.
The mean distance from the plane is very similar for the two sub-samples, 1210 pc for the Layden et al. stars and 926 pc for the BSL stars. We therefore estimate g ∼ ±0.2, and hence δη/η ∼ ±0.02. Since this error is uncorrelated with the statistical error, we add the two in quadrature to obtain our final estimate for the non-kinematic sample η = 0.953 ± 0.054.
Next, to either compare or combine the two samples, we must evaluate the correlation coefficient of the two determinations, γ ≡ c 12 (c 11 c 22 ) −1/2 . Here c ij is the covariance matrix of the two measurements, with c 11 = (0.057) 2 (from Table 5) and c 22 = (0.054) 2 (derived in the previous paragraph). In the appendix, we find γ = 0.46. This is very close to 0.5, the value one would naively guess because the radial-velocity measurements are almost completely independent while the proper motions are almost completely dependent.
The 1 σ expected difference between the two measurements is therefore Finally, we note that Table 5 contains the first complete solution of the velocity ellipsoid of halo stars as determined from a non-kinematically selected sample. One important (though not unexpected) result is that the stellar halo is not moving relative to the LSR in either the radial or vertical directions. From Table 5, the estimates would appear to be 3 ± 8 km s −1 away from the Galactic center and 0 ± 5 km s −1 toward the north Galactic pole.
In fact, each of these values (and errors) should be augmented by a factor (1 − f ) −1 ∼ 1.25 because the contaminating thick disk stars are known to be bound to the disk and so have zero mean motion in both directions. Hence the best estimates for the halo bulk motion in the two directions are 4 ± 10 km s −1 (radial) and 0 ± 6 km s −1 (vertical). The rotation of the halo (35 ± 8 km s −1 ) is larger than in the solution for the 162 RR Lyrae alone, but this is to be expected because that sample was selected in part kinematically, that is by eliminating stars with the most prograde orbits (Layden et al. 1996).

Discussion and Conclusions
We have investigated many potential sources of systematic error in the RR Lyrae absolute magnitude calibration by statistical parallax using a combination of analytic and Monte Carlo techniques. We find that all corrections to previous results are small and, in particular, that the highly non-Gaussian RR Lyrae velocity distribution does not bias the determination at all even though (or rather, because) the method explicitly assumes a Gaussian distribution. We find that the mean RR Lyrae absolute magnitude is We also analyze a semi-independent non-kinematically selected sample of stars with metallicities [Fe/H]≤ −1.5 taken from Layden et al. and BSL and find similarly M V = 0.79 ± 0.12 at [Fe/H] = −1.79. Additionally, this analysis yields measurements of the radial bulk motion (4 ± 10 km s −1 ) and vertical bulk motion (0 ± 6 km s −1 ) of the halo relative to the LSR.
Our principal result is therefore that the RR Lyrae absolute magnitude calibration by statistical parallax is extremely robust and that the statistical error (0.13 mag) should be taken at face value.
If one assumes that type ab RR Lyrae stars in the LMC are (apart from different mean metallicity) similar to those in the Layden et al. (1996) sample, then the distance modulus to the LMC is where V 0,LMC is the dereddened mean apparent magnitude of RRab's at the distance of the center of mass of the LMC. Ideally, the way to determine this quantity is to measure The error is from the scatter and hence accounts for measurement errors and uncertainty in the mean distance of the clusters relative to the center of mass. We estimate an additional error in the mean extinction of 0.03. Applying equation (68), we find µ LMC = 18.28 ± 0.14.
However, two lines of evidence suggest that cluster RR Lyrae stars may be systematically brighter than those in the field. First, Sweigart (1997) has argued that abundance anomalies seen in cluster giants but not in halo field giants could be due to rotational mixing in the former. Such mixing would dredge up helium in horizontal branch stars, making them up to a few tenths of a magnitude brighter. Second, new main-sequence-fitting distances to Galactic globulars (Reid 1997;Gratton et al. 1997) (1992) which, from their Figure 5b, has an uncertainty of ±0.015 for each individual star. The mean period of the Graham (1977) stars is P = 0.564 days, implying (B − V ) 0 = 0.36.
Using and an axis ratio c/a = 0.6, and taking the inclination of the LMC to be i = 27 • (Bessel, respectively. This leads to two estimates for the LMC center of mass µ N1763 field LMC = 18.13±0.15 and µ N2210 field LMC = 18.31 ± 0.14. Since the errors of these determinations are correlated, the difference is ∆µ = 0.18 ± 0.09, which is uncomfortably large. While both of these determinations are much smaller than the traditional value µ LMC = 18.5 and are smaller yet compared to several new determinations (Reid 1997;Gratton et al. 1997;Feast & Catchpole 1997), there remains the possibility that the underlying measurements of V 0 are affected by systematic errors. A robust estimate of this quantity should therefore await results from the microlensing surveys.

A. Determination of the Correlation Coefficient
We now return to the framework of §2 in order to evaluate the correlation coefficient γ ≡ c 12 (c 11 c 22 ) −1/2 . The error for a single determination of η has the form given by equation (66). For two partially dependent determinations of η, this generalizes to c ij = 3 2 (6 + κ 2 ) −1 (2c r ij + c t ij ) where c ij is the 2 × 2 covariance matrix of the errors in η, and similarly for c t . Here N r,1 = 162 is the number of radial-velocity stars in the Layden et al. sample, N r,2 = 830 is the number in the non-kinematic sample, and N r,12 = 103 is the number of overlap stars. Similarly, N t,1 = 162, N t,2 = 106, and N t,12 = 103. It is immediately clear that the diagonal elements of equation (A1) reduce to equation (66).
The off-diagonal terms take account of the covariance between the two radial-velocity-based velocity ellipsoids and the covariance between the two tangential-velocity-based velocity ellipsoids. However, equation (A1) takes account only of the statistical errors. We now take account of the additional systematic error caused by differering levels of thick disk contamination. Recall that this error affects only the non-kinematic determination, so it is included by adding a matrix to equation (A1) where gf ∆ = 0.02 was evaluated in §7. Hence, γ = c 12 (c 11 c 22 ) −1/2 = 0.46. Right panel shows that the distribution in radial direction is very boxy. Note. -As previously, the name of the case uniquely characterizes statistical properties of velocity input. The first row, just below column description, gives the errors predicted analytically in §4.  NOTE.|The bulk motion and velocity dispersions are given in km s 1 and the o -diagonal coe cients of the velocity correlation matrix are shown in normalized dimensionless form:C ij = C ij (C ii C jj ) 1 2 . |The input values of the parameters for the underlying distribution of stars are given below the descriptions of the column content. The values in parentheses in the rst row give the approximate errors in the determinations of the biases based on the NNRG case. NOTE.|The bulk motion and velocity dispersions are given in km s 1 and the o -diagonal coe cients of the velocity correlation matrix are shown in normalized dimensionless form:C ij = C ij (C ii C jj )