integrating out nuisances

Further insipired by yesterday's post about binary fitting, I worked today on the treatment of nuisance parameters that have known distributions. These can be treated as noise sometimes. Let me explain:

If I had to cartoon inference (or measurement) in the face of nuisance parameters, I would say that frequentists profile (optimize) over the nuisances and Bayesians marginalize (integrate) over the nuisances. In general frequentists cannot integrate over anything, because there is no measure in any of the parameter spaces. But sometimes there is a measure. In particular, when there is a compact symmetry:

We know (or very strongly believe) that all possible orientations of a binary-star orbit are equally likely. In this model (or under this normal assumption) we have a distribution over two angles (theta and phi for that orbit pole, say); it is the distribution set by the compact group SO(2). Thus we can treat the orientation as a noise source with known distribution and integrate over it, just like we would any other noise source. So, in this case (and many cases like it) we can integrate (marginalize) even as frequentists. That is, there are frequentism-safe marginalizations possible in binary-star orbit fitting. This should drop the 12-parameter fits (for ESA Gaia data) down to 8-parameter, if I have done my math right.

binary stars with periods of exactly one year

On Friday, Kareem El-Badry (Caltech) gave a seminar about looking for (and finding!) stars in binary orbits around dark or much darker companions, like black holes, neutron stars, and white dwarfs. He showed results that involve ESA Gaia astrometry, where he noted that the Gaia Mission has no sensitivity to periods right at (or within an inverse mission-length frequency difference of) one-year periods (inverse year frequencies). After the talk I objected that these are not exactly degenerate; El-Badry said that the inferences blow up there.

I spent some time on the weekend thinking about this point, and I now understand it: There is a particular one-year orbit that a star can have (around a darker companion) such that the photocenter of the system makes a motion that is identical to the apparent parallax motion. Thus there is an exact degeneracy between the parallax and a certain one-year orbit.

Does that mean that we can't measure orbits at one year (or, for that matter, parallaxes)? No, it does not. After all, the parallax ellipse has a particular celestial (angular) shape and phase. But it might require some kind of reparameterization of orbits near one-year periods. I think I know how to do that. Should we find the missing binaries? (Oh and by the way, this degeneracy means that, in a strict frequentist sense, Gaia can't measure parallaxes at all without additional information.)

how significant is your anomaly?

So imagine that you have a unique data set Y, and in that data set Y you measure a bunch of parameters θ by a bunch of different methods. Then you find, in your favorite analysis, your estimate of one particular parameter is way out of line: All of physics must be wrong! How do you figure out the significance of your result?

If you only ever have data Y, you can't answer this question very satisfactorily: You searched Y for an anomaly, and now you want to test the significance. That's why so many a posteriori anomaly results end up going away: That search probably tested way more hypotheses than you think it did, so any significances should be reduced accordingly.

The best approach is to use only part of your data (somehow) to search, and then use a found anomaly to propose a hypothesis test, and then test that test in the held-out or new data. But that often isn't possible, or it is already too late. But if you can do this, then there is usually a likelihood ratio that is decisive about the significance of the anomaly!

I discussed all these issues today with Kate Storey-Fisher (Stanford) and Abby Williams (Chicago) today, as we are trying to finish a paper on the anomalous amplitude of the kinematic dipole in quasar samples.

finding emission lines (and other oddities) in hot stars

I showed my robust spectral decomposition (dimensionality reduction) and residuals to the MPIA Binaries group today. There was much useful feedback (including that my H-gamma was actually H-delta; embarassing!). One comment was that the model isn't truly a causal separation between star and lines, so there will be some mean lines in the star model; lines aren't entirely outliers. That's true! The group suggested that I iterate to remove stars with lines from the training set.

After the meeting, I implemented some of that, but problems like this have a pathology: If you carefully remove stars with high residuals at some wavelength, then the training data will be deficient, or low, at that wavelength. And then the model will go lower, and then more stars will have excess at that wavelength and: Disaster. So when I implemented, I required a 2-sigma deviation, and I removed both high and low outliers. I don't know if this will work, but I am testing now.

wrote like the wind; frequentist vs Bayes on sparsity

My goal this year in Heidelberg is to move forward all writing projects. I didn't really want to start new projects, but of course I can't help myself, hence the previous post. But today I crushed the writing: I wrote four pages in the book that Rix (MPIA) wants me to write, and I got more than halfway done with a Templeton Foundation pre-proposal that I'm thinking about, and I partially wrote up the method of the robust dimensionality reduction that I was working on over the weekend. So it was a good day.

That said, I don't think that the iteratively reweighted least squares implementation that I am using in my dimensionality reduction has a good probabilistic interpretation. That is, it can't be described in terms of a likelihood function. This is related to the fact that frequentist methods that enforce sparsity (like L1 regularization) don't look anything like Bayesian methods that encourage sparsity (like massed priors). I don't know how to present these issues in any paper I try to write.

young stars in SDSS-V

Over the last two weeks, I built a new robust dimensionality-reduction method called Robust-HMF. This method is a hammer, looking for a nail. This week, Hans-Walter Rix (MPIA) suggested that I use the method on young stellar objects observed in SDSS-V BOSS spectra. I did that, and I found hundreds of young stars with narrow H-alpha emission lines. It turns out that the Robust-HMF method does really well on these data; it can fit all the absorption lines in the stars, plus all the wacky continuum shapes generated by a combination of instrumental effects and dust attenuation. That is, I can run the method on more-or-less raw data, provided that it has been shifted to rest frame (and the raw SDSS-V pipelines do that pretty well on stars like these).

Anyway, I don't know what to say about my results quantitatively yet, but I find hundreds of emission-line stars, and my equivalent-width sensitivity is excellent. Who wants this sample?

SPHEREx data

Dustin Lang (Perimeter) and I spoke today about many things, but the conversation got de-railed when Lang showed me his visualizations of the brand-new NASA SPHEREx data. Oh. My. Goodness. First of all, the data are being released daily, before the team has done its analysis, so that anyone in the world can do anything with it. Talk about open! Talk about international! Talk about everything I believe about astrophysics! Also, the data are well documented, high in signal to noise, and released with good and useful metadata. This is killer. Of course Lang has already made a viewer for it and can do all comparisons to the DESI data. Get in touch with him if you want tools.

should I write a book?

Is it research to have a set of conversations about whether to write a book on data analysis? Hans-Walter Rix (MPIA) thinks I should put together my arXiv-only submissions plus a lot more into a book about data analysis. His point of view is that the most important thing is how to convert an ill-posed question about the Universe into a well-posed operation on data.

robust matrix factorization

There is a very nice algorithm and set of methods called “Robust PCA,” originating in a paper by Candès. This method makes use of ideas from convex optimization to simultaneously learn a low-rank representation of the data plus a sparse representation for the outliers. This kind of situation comes up in astronomy all the time.

Way back, Tsalmantza and I made a replacement for PCA called HMF that deals with heteroskedastic data (variable error bars or variable data weights) and also missing data; we used it to build low-rank models of quasar spectra. This weekend I built a robust version of HMF (Robust HMF maybe?) that uses ideas from iteratively reweighted least squares to mimic the algorithm behind Robust PCA. It works! And it works well. Unfortunately, right now the investigator has to tune the rank of the low-rank part and also the soft outlier cutoff used in the IRLS. I would love to figure out principled ways to choose both. If you want to follow along, development is happening here for now.

is it surprising that there are high-redshift supermassive black holes?

A nice talk at MPIA by Hanna Ũbler (MPE) about very high-redshift (redshifts 8 to 14 even) galaxies and their black-hole contents started some nice discussions in the audience and afterwards about the formation of black holes. Because of the Eddington limit (which is a limit on luminosity), black-hole growth is probably limited. The limit is on luminosity, not mass accretion rate, and the relationship between these is the radiative efficiency. As the efficiency goes down, the mass accretion rate of an Eddington accretor goes up. So the question: When can you first form a super-massive (10 million solar masses, say) black hole is a question simultaneously about seed black holes and about radiative efficiency. Anyway, this is all unfortunate, because if the efficiency couldn't get very low, then the black holes we find with NASA JWST would already be putting very strong pressure on fundamental physics in the early Universe.

coherent oscillator injection and recovery

Coherent oscillators—astronomical sources that pulse or oscillate in a phase-stable way over long timescales—have been useful astrophysical tools. For two examples: Stably pulsing pulsars were used to discover gravitational radiation (and are being used to find the stochastic background). Delta-scuti star asteroseismic modes were used to find orbital companions. This led undergrad Nana Miller (NYU) and me and others to look for all the coherent modes we can find among all the stars in the NASA Kepler Mission data. We have have technology to find modes, and we have technology to test for coherence. Now we have to do injection–recovery tests to estimate our detection limits. Today I delivered a simple plan for doing injections.

Our plan is to produce a catalog, not of stars but of modes, every one of which has a coherence time that is longer than the lifetime (4 years) of the Kepler Mission. Then: What do we use them for?

likelihood ratios not posteriors, please

There is an informal meeting at MPIA every Wednesday regarding binary stars, with a bit of a focus on massive binaries. Today there was a very nice presentation by Jakob Stegmann (MPA) about some anomalies among the (only six) black-hole–neutron-star binaries discovered by NSF LIGO. He showed the example of GW 200105, which shows a large eccentricity (0.15-ish). This eccentricity is very hard to explain, given how the inspirals evolve as they radiate. But the analysis of the eccentricity (from perhaps this paper) is Bayesian, so it isn't clear how much the eccentricity result is forced by the data and how much is forced by the prior over nuisance parameters. That's one of the main points of my forthcoming paper on measurement. I think maybe I should just re-analyze this one with a profile likelihood. I hope the data and code are public!

robust dimensionality reductions

Dimensionality reduction (the basic being PCA) is very sensitive to outliers: A single bad pixel can dominate most objectives and thus create a spurious dimension. One of the best and most classic solutions to this is the robust PCA method, which is presented in a (very long) paper with impressive math and beautiful results. Yesterday Hans-Walter Rix (MPIA) and I coded it up and applied it to ESA Gaia RVS spectra, with extensive (and impressive) help from Claude. It looks very promising, especially in capturing oddities in hot stars. Today I worked out that there should be something similar that takes into account data weights (inverses of squared uncertainties), and I wrote down the algorithm (on paper). We'll see.

stellar twins vs synthetic stellar twins

In the Milky Way meeting at MPIA today, a bit of a discussion broke out about using stellar twins, inspired by work by Yuan-Sen Ting (OSU). The idea is: If you have two stars with very similar overall metallicity, and very similar temperature and surface gravity, then it should be possible to measure accurate element abundnace anomalies between the two stars, even in the absence of an extremely accurate spectral synthesis code.

My view, which does not contradict this point, is that an even better way to use this stellar-twin idea is to synthesize a twin for every star, using stars that are similar in (either) parameters or else spectral space. After all, an interpolation to your target star should more accurately represent it than even the most similar individual comparison star. That idea, fundamentally, is the main idea behind The Cannon.

how did the Solar System form?

I saw a very nice talk today by Philippine Griveaud (MPIA) about how the Solar System formed. The idea is that the giant planets formed in an accretion disk. Their formation opened gaps and caused migration (first Type I and then Type II, if you must know :). That migration pulled them into a resonant chain. That is, if the giant planets formed the way we think they formed, they must have been in a resonant chain. But they aren't in such a chain now; what gives?

The idea is that when the gas is expended (or blown out by winds), the remaining planetestimals (think: asteroids, comets, Kuiper Belt objects) interact with the planets such that they get moved from orbit to orbit and eventually ejected. These dynamical interactions break the resonant chain, migrate the giant planets to their current locations, and scatter rocks and ice balls into the interstellar regions.

It was a great talk, but also led to a lot of interesting questions, such as: How does this all fit in with the formation of the rocky planets? And how does this square with our observations (growing rapidly, apparently) of interstellar asteroids? Oh and: How does all this connect to observations of debris disks, which I now (officially) love.