The AstroStat Slog » Equation of the Week http://hea-www.harvard.edu/AstroStat/slog Weaving together Astronomy+Statistics+Computer Science+Engineering+Intrumentation, far beyond the growing borders Fri, 09 Sep 2011 17:05:33 +0000 en-US hourly 1 http://wordpress.org/?v=3.4 Blackbody Radiation [Eqn] http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-blackbody/ http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-blackbody/#comments Wed, 27 Aug 2008 17:00:22 +0000 vlk http://hea-www.harvard.edu/AstroStat/slog/?p=507 Like spherical cows, true blackbodies do not exist. Not because “black objects are dark, duh”, as I’ve heard many people mistakenly say — black here simply refers to the property of the object where no wavelength is preferentially absorbed or emitted, and all the energy input to it is converted into radiation. There are many famous astrophysical cases which are very good approximations to perfect blackbodies — the 2.73K microwave background radiation left over from the early Universe, for instance. Even the Sun is a good example. So it is often used to model the emission from various objects.

The blackbody spectrum is
$$B_{\nu}(T) = \frac{2 h \nu^3}{c^2} \frac{1}{e^{h \nu / k_B T} – 1} ~~ {\rm [erg~s^{-1}~cm^{-2}~Hz^{-1}~sr^{-1}]} \,,$$
where ν is the frequency in [Hz], h is Planck’s constant, c is the speed of light in vacuum, and kB is Boltzmann’s constant. The spectrum is interesting in many ways. Its shape is characterized by only one parameter, the radiation temperature T. A spectrum with a higher T is greater in intensity at all frequencies compared to one with a lower T, and the integral over all frequencies is σ T4, where $$\sigma \equiv \frac{2\pi^5k_B^4}{15 c^2 h^3}$$ is the Stefan-Boltzmann constant. Other than that, the normalization is detached, so to speak, from T, and differences in source luminosities are entirely attributable to differences in emission surface area.

The general shape of a blackbody spectrum is like a rising parabola at low ν (which led to much hand-wringing in the late 19th century about the Ultraviolet Catastrophe) and an exponential drop at high ν, with a well-defined peak in between. The frequency at which the spectrum peaks is dependent on the temperature, with
$$\nu_{\rm max} = 2.82 \frac{k_B T}{h}$$,
or equivalently,
$$ \lambda_{\rm max} = \frac{2.9\cdot10^4}{T} ~~ {\rm[\AA]} \,,$$
where T is in [degK].

]]>
http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-blackbody/feed/ 0
Magnitude [Eqn] http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-magnitude/ http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-magnitude/#comments Wed, 20 Aug 2008 17:00:34 +0000 vlk http://hea-www.harvard.edu/AstroStat/slog/?p=447 I still remember my first class as a new grad student. As a cocky Physics graduate, I was quite sure I knew plenty of astronomy. Astro 301, class 1, and it took all of 20 minutes of talk about stellar magnitudes to put that notion to permanent rest. So, for the sake of our stats colleagues, here’s a brief primer on one of the basic building blocks of astronomy.

For historical reasons, astronomers measure the brightness of celestial objects in rank order. The smaller the rank number, aka magnitude, the brighter the object. Thus, a star of the first magnitude is much brighter than a star of the sixth magnitude, and it would take exceptionally good eyes and a very dark sky to see a star of the seventh magnitude. Now, it turns out that the human eye perceives brightness on a log scale, so magnitudes are numerically similar to log(brightness). And because they are a ranking list, it is always with reference to a standard. After some rough calibration to match human perception to true brightness of stars in the night sky, we have a formal definition for magnitude,
$$m = – \frac{5}{2}\log_{10}\left(\frac{f_{object}}{f_{standard}}\right) \,,$$
where fobject is the flux from the object and fstandard is the flux from a fiducial standard. In the optical bands, the bright star Vega (α Lyrae) has been adopted as the standard, and has magnitudes of 0 in all optical filters. (Well, not exactly because Vega is not constant enough, and as a practical matter there is nowadays a hierarchy of photometric standard stars that are accessible at different parts of the sky.) Note that we can also write this in terms of the intrinsic luminosity Lobject of the object and the distance d to it,
$$m = – \frac{5}{2}\log_{10}\left(\frac{L_{object}}{4 \pi d^2}\frac{1}{f_{standard}}\right) \,.$$

Because astronomical objects are located at a vast variety of distances, it is useful to define an intrinsic magnitude of the object, independent of the distance. Thus, in contrast to the apparent magnitude m, which is the brightness at Earth, an absolute magnitude is defined as the brightness that would be perceived if the object were 10 parsecs away,
$$M \equiv m|_{d={\rm 10~pc}} = m – \frac{5}{2}\log_{10}\left(\frac{d^2}{{\rm (10~pc)}^2}\right) \equiv m – 5\log_{10}d + 5$$
where d is the distance to the object in [parsec], and the squared term is of course because of the inverse square law.

There are other issues such as interstellar absorption, cosmological corrections, extent of the source, etc., but let’s not complicate it too much right away.

Colors are differences in the magnitudes in different passbands. For instance, if the apparent magnitude in the blue filter is mB and in the green filter is mV (V for “visual”), the color is mB-mV and is usually referred to as “B-V” color. It is the difference in magnitudes, and is related to the log ratio of the intensities.

For an excellent description of what is involved in the measurement of magnitudes and colors, see this article on analyzing photometric data by Star Stryder.

]]>
http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-magnitude/feed/ 0
Differential Emission Measure [Eqn] http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-dem/ http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-dem/#comments Wed, 13 Aug 2008 17:00:13 +0000 vlk http://hea-www.harvard.edu/AstroStat/slog/?p=421 Differential Emission Measures (DEMs) are a summary of the temperature structure of the outer atmospheres (aka coronae) of stars, and are usually derived from a select subset of line fluxes. They are notoriously difficult to estimate. Very few algorithms even bother to calculate error envelopes on them. They are also subject to numerous systematic uncertainties which can play havoc with proper interpretation. But they are nevertheless extremely useful since they allow changes in coronal structures to be easily discerned, and observations with one instrument can be used to derive these DEMs and these can then be used to predict what is observable with some other instrument.

The flux at Earth due to an atomic transition u –> l from a volume element δV at a location ɼ,

Iul(ɼ) = (1/4 π) (1/d(ɼ)2) A(Z,ɼ) Gul(ne(ɼ),Te(ɼ)) ne(ɼ)2 δV(ɼ) ,

where ne is the electron density and Te is the temperature of the plasma, A(Z,ɼ) are the abundance of element Z, Gul(ne,Te) is the atomic emissivity for the transition, and d is the distance to the source.

We can combine the flux from all the points in the field of view that arise from plasma at the same temperature,

Iul(Te) = (1/4 π) ∑ɼ|Te (1/d(ɼ)2) A(Z,ɼ) Gul(ne(ɼ),Te) ne2δV(ɼ) .

Assuming that A(Z,ɼ), ne(ɼ) do not vary over the points in the summation,

Iul(Te) ≈ (1 / 4 π d2) Gul(ne,Te) A(Z) ne2 (ΔV / Δlog Te) Δlog Te ,

and hence the total line flux due to emission at all temperatures,


Iul = ∑Te (1 / 4 π d2) A(Z) Gul(ne,Te) DEM(Te) ΔlogTe

The quantity

DEM(Te) = ne2 (ΔV / Δlog Te)

is called the Differential Emission Measure and is a very useful summary of the temperature structure of stellar coronae. It is typically reported in units of [cm-3] (or [cm-5] if ΔV is written out as area*Δh). Sometimes it is defined as ne2(ΔV/ΔT) and has units [cm-3K-1].

The expression for the line flux is an instance of Fredholm’s Equation of the First Kind and the DEM(Te) solution is thus unstable and subject to high-frequency oscillations. There is a whole industry that has grown up trying to derive DEMs from often highly unreliable datasets.

]]>
http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-dem/feed/ 2
I Like Eq http://hea-www.harvard.edu/AstroStat/slog/2008/i-like-eq/ http://hea-www.harvard.edu/AstroStat/slog/2008/i-like-eq/#comments Wed, 13 Aug 2008 16:59:54 +0000 vlk http://hea-www.harvard.edu/AstroStat/slog/?p=379 I grew up in an environment that glamourized mathematical equations. Equations adorned a text like jewelry, set there to dazzle, and often to outshine the text that they were to illuminate. Needless to say, anything I wrote was dense, opaque, and didn’t communicate what it set out to. It was not until I saw a Reference Frame essay by David Mermin on how to write equations (1989, Physics Today, 42, p9) that I realized that equations should be treated as part of the text. You should be able to read them. David Mermin set out 3 rules for writing out equations, which I’ve tried to follow diligently (if not always successfully) since then.

  1. Number or label all displayed equations (Fisher’s Rule):

    The most common violation of Fisher’s rule is the misguided practice of numbering only those displayed equations to which the text subsequently refers back. … it is necessary to state emphatically that Fisher’s rule is for the benefit not of the author, but the reader.
    For although you, dear author, may have no need to refer in your text to the equations you therefore left unnumbered, it is presumptuous to assume the same disposition in your readers. And although you may well have acquired the solipsistic habit of writing under the assumption that you will have no readers at all, you are wrong.

  2. When referring to an equation within the text, identify it by a phrase as well as a number (aka the Good Samaritan Rule):

    A Good Samaritan is compassionate and helpful to one in distress, and there is nothing more distressing than having to hunt your way back in a manuscript in search of Eq. (2.47), not because your subsequent progress requires you to inspect it in detail, but merely to find out what it is about so you may know the principles that go int othe construction of Eq. (7.38).

  3. Punctuate the equation (aka the Math is Prose Rule):

    The equations you display are embedded in your prose and constitute an inseparable part of it.

    Regardless … of how to parse the equation internally, certain things are clear to anyone who understands the equation and the prose in which it is embedded.

    We punctuate equations because they are a form of prose (they can, after all, be read aloud as a sequence of words) and are therefore subject to the same rules as any other prose. … punctuation makes them easier to read and often clarifies the discussion in which they occur. … viewing an equation not as a grammatically irrelevant blob, but as a part of the text … can only improve the fluency and grace of one’s expository mathematical prose.

]]>
http://hea-www.harvard.edu/AstroStat/slog/2008/i-like-eq/feed/ 2
Background Subtraction, the Sequel [Eqn] http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-bkgsubtract-poisson/ http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-bkgsubtract-poisson/#comments Wed, 06 Aug 2008 17:00:39 +0000 vlk http://hea-www.harvard.edu/AstroStat/slog/?p=407 As mentioned before, background subtraction plays a big role in astrophysical analyses. For a variety of reasons, it is not a good idea to subtract out background counts from source counts, especially in the low-counts Poisson regime. What Bayesians recommend instead is to set up a model for the intensity of the source and the background and to infer these intensities given the data.

For instance, suppose as before, that C counts are observed in a region of the image that overlaps a putative source, and B counts in an adjacent, non-overlapping region that is mostly devoid of the source and which is r times larger in area and exposure than the source region. Further suppose that a fraction f of the source falls in the so-called source region (typically, f~0.9) and a fraction g falls in the background region (we strive to make g~0). Then the observed counts can be written as Poisson realizations of intensities,

C = Poisson(φS) ≡ Poisson(f ­ θS + θB) , and
B = Poisson(φB) ≡ Poisson(g ­ θS + r ­ θB)
,

where the subscripts denote the model intensities in the source (S) or background (B) regions.

The joint probability distribution of the model intensities,

p(φS φB | C B) dφSB

can be rewritten in terms of the interesting parameters by transforming the variables,

≡ p(θS θB | C B) J(φS φB ; θS θB) d θS d θB

where J(φS φB ; θS θB) is the Jacobian of the coordinate transformation, and thus

= p(θS θB | C B) (r ­ f – g) d θS d θB .

The posterior probability distribution of the source intensity, θS, can be derived by marginalizing this over the background intensity parameter, θB. A number of people have done this calculation in the case f=1,g=0 (e.g., Loredo 1992, SCMA II, p275; see also van Dyk et al. 2001, ApJ 584, 224). The general case is slightly more convoluted, but is still a straightforward calculation (Kashyap et al. 2008, AAS-HEAD 9, 03.02); but more on that another time.

]]>
http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-bkgsubtract-poisson/feed/ 0
keV vs keV [Eqn] http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-kev-kev/ http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-kev-kev/#comments Wed, 30 Jul 2008 17:00:20 +0000 vlk http://hea-www.harvard.edu/AstroStat/slog/?p=360 I have noticed that our statistician collaborators are often confused by our units. (Not a surprise; I, too, am constantly confused by our units.) One of the biggest culprits is the unit of energy, [keV], which is 1000 electron Volts, for the energy acquired by an electron when it falls through an electric potential of 1 Volt:

1 [eV] ≡ 1.6021892 · 10-19 [Joule] ≡ 1.6021892 · 10-12 [erg] .

The confusion is because the same units are used to denote two separate quantities which happen to have similar magnitudes for a commonly encountered spectral model, Bremsstrahlung emission.

  1. the frequency ν, or wavelength λ, of a photon: As Planck discovered, the energy of a photon is directly related to the frequency ν,

    E = h · ν ≡ h · c / λ ,

    where h=6.6261760 · 10-27 [erg s] is Planck’s constant and c=2.9979246 · 1010 [cm s-1] is the speed of light in vaccum. When λ is given in [Ångström] ≡ 10-8 [cm], we can convert it as

    [keV] = 12.398521 / [Å] ,

    which is an extraordinarily useful thing to know in high-energy astrophysics.

  2. the temperature T of a gas or plasma: Here we look to thermodynamics, which relates the kinetic energy of random motion of particles in a gas to a gross property, the temperature of the gas,

    E = kB · T ,

    where kB = 1.3806620 · 10-16 [erg K-1] is Boltzmann’s constant. Then, a temperature in degrees Kelvin can be written in units of keV by converting it with the formula

    [keV] = 8.6173468 · 10-8 · [K] ≡ 0.086173468 · [MK] .

It is tempting to put the two together and interpret a temperature as a photon energy. This is possible for the aforementioned Bremsstrahlung radiation, where plasma at a temperature T produces a spectrum of photons distributed as e-h ν / kB T and it is possible to tie the temperature to the photon energy at the point where the numerator and denominator have the same numerical value. For example, a 1 keV (temperature) Bremsstrahlung spectrum extends out to 1 keV (photon energy). X-ray Astronomers use this as shorthand all the time, and it confuses the hell out of everybody else.

]]>
http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-kev-kev/feed/ 1
The Banff Challenge [Eqn] http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-banff-challenge/ http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-banff-challenge/#comments Wed, 23 Jul 2008 17:00:48 +0000 vlk http://hea-www.harvard.edu/AstroStat/slog/?p=357 With the LHC coming on line anon, it is appropriate to highlight the Banff Challenge, which was designed as a way to figure out how to place bounds on the mass of the Higgs boson. The equations that were to be solved are quite general, and are in fact the first attempt that I know of where calibration data are directly and explicitly included in the analysis.

The observables are counts N, Y, and Z, with

N ~ Pois(ε λS + λB) ,
Y ~ Pois(ρ λB)
,
Z ~ Pois(ε υ)
,

where λS is the parameter of interest (in this case, the mass of the Higgs boson, but could be the intensity of a source), λB is the parameter that describes the background, ε is the efficiency, or the effective area, of the detector, and υ is a calibrator source with a known intensity.

The challenge was (is) to infer the maximum likelihood estimate of and the bounds on λS, given the observed data, {N, Y, Z}. In other words, to compute

p(λS|N,Y,Z) .

It may look like an easy problem, but it isn’t!

]]>
http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-banff-challenge/feed/ 4
chi-square distribution [Eqn] http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-chisq-distribution/ http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-chisq-distribution/#comments Wed, 16 Jul 2008 17:00:16 +0000 vlk http://hea-www.harvard.edu/AstroStat/slog/?p=342 The Χ2 distribution plays an incredibly important role in astronomical data analysis, but it is pretty much a black box to most astronomers. How many people know, for instance, that its form is exactly the same as the γ distribution? A Χ2 distribution with ν degrees of freedom is

p(z|ν) = (1/Γ(ν/2)) (1/2)ν/2 zν/2-1 e-z/2 ≡ γ(z;ν/2,1/2) , where z=Χ2.

Its more familiar usage is in the cumulative form, which is just the incomplete gamma function. This is where you count off how much area is enclosed in [0,Χ2) to tell at what point the 68%, 95%, etc., thresholds are met. For example, for ν=1,

0Z dx p(Χ2|ν=1) = 0.68 when Z=1.

This is the origin of the ΔΧ2=1 method to determine error bars on best-fit parameters.

]]>
http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-chisq-distribution/feed/ 4
Kaplan-Meier Estimator (Equation of the Week) http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-kaplan-meier/ http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-kaplan-meier/#comments Wed, 09 Jul 2008 17:00:54 +0000 vlk http://hea-www.harvard.edu/AstroStat/slog/?p=356 The Kaplan-Meier (K-M) estimator is the non-parametric maximum likelihood estimator of the survival probability of items in a sample. “Survival” here is a historical holdover because this method was first developed to estimate patient survival chances in medicine, but in general it can be thought of as a form of cumulative probability. It is of great importance in astronomy because so much of our data are limited and this estimator provides an excellent way to estimate the fraction of objects that may be below (or above) certain flux levels. The application of K-M to astronomy was explored in depth in the mid-80′s by Jurgen Schmitt (1985, ApJ, 293, 178), Feigelson & Nelson (1985, ApJ 293, 192), and Isobe, Feigelson, & Nelson (1986, ApJ 306, 490). [See also Hyunsook's primer.] It has been coded up and is available for use as part of the ASURV package.

Consider a simple case where you have N observations of the luminosities of a source. Let us say that all N sources have been detected and their luminosities are estimated to be Li, i=1..N, and that they are ordered such that Li < Li+1 Then, it is easy to see that the fraction of sources above each Li can be written as the sequence

{ N-1, N-2, N-3, … 2, 1, 0}/N

The K-M estimator is a generalized form that describes this sequence, and is written as a product. The probability that an object in the sample has luminosity greater than Lk is

S(L>L1) = (N-1)/N
S(L>L2) = (N-1)/N * ((N-1)-1)/(N-1) = (N-1)/N * (N-2)/(N-1) = (N-2)/N
S(L>L3) = (N-1)/N * ((N-1)-1)/(N-1) * ((N-2)-1)/(N-2) = (N-3)/N

S(L>Lk) = Πi=1..k (ni-1)/ni = (N-k)/N

where nk are the number of objects still remaining at luminosity level L ≥ Lk, and at each stage one object is decremented to account for the drop in the sample size.

Now that was for the case when all the objects are detected. But now suppose some are not, and only upper limits to their luminosities are available. A specific value of L cannot be assigned to these objects, and the only thing we can say is that they will “drop out” of the set at some stage. In other words, the sample will be “censored”. The K-M estimator is easily altered to account for this, by changing the decrement in each term of the product to include the censored points. Thus, the general K-M estimator is

S(L>Lk) = Πi=1..k (ni-ci)/ni

where ci are the number of objects that drop out between Li-1 and Li.

Note that the K-M estimator is a maximum likelihood estimator of the cumulative probability (actually one minus the cumulative probability as it is usually understood), and uncertainties on it must be estimated via Monte Carlo or bootstrap techniques [or not.. see below].

]]>
http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-kaplan-meier/feed/ 13
Poisson Likelihood [Equation of the Week] http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-poisson-likelihood/ http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-poisson-likelihood/#comments Wed, 02 Jul 2008 17:00:32 +0000 vlk http://hea-www.harvard.edu/AstroStat/slog/?p=333 Astrophysics, especially high-energy astrophysics, is all about counting photons. And this, it is said, naturally leads to all our data being generated by a Poisson process. True enough, but most astronomers don’t know exactly how it works out, so this derivation is for them.

Suppose N counts are randomly placed in an interval of duration τ without any preference for appearing in any particular portion of τ. i.e., the distribution is uniform. The counting rate R = N/τ. We can now ask, what is the probability of finding k counts in an infinitesimal interval δt within τ?

First, consider the probability that one count, placed randomly, will fall inside δt,

ρ = δt/τ ≡ Rδt/N ≡ ν/N


where ν = R δt represents the expected counts intensity in the interval δt. When N counts are scattered over τ, the probability that k of them will fall inside δt is described with a binomial distribution,

p(k|ρ,N) = NCk ρk (1-ρ)N-k


as the product of the probability of finding k events inside δt and the probability of finding the remaining events outside, summed over all the possible distinct ways that k events can be chosen out of N. Expanding the expression and rearranging,

= N!/{(N-k)!k!} (R δt/N)k (1-(R δt/N))N-k

= N!/{(N-k)!k!} (νk/Nk) (1-(ν/N))N-k

= N!/{(N-k)!Nk} (νk/k!) (1-(ν/N))N (1-(ν/N))-k


Note that as N,τ —> ∞ (while keeping R fixed),

N!/{(N-k)!Nk} , (1-(ν/N))-k —> 1
(1-(ν/N))N —> e


and the expression reduces to

p(k|ν) = (νk/k!) e


which is the familiar (in a manner of speaking) expression for the Poisson likelihood.

]]>
http://hea-www.harvard.edu/AstroStat/slog/2008/eotw-poisson-likelihood/feed/ 1