oscillating stars as clocks

On a walk up the mountain to MPIA, Rix and I discussed quantitatively the question of whether stellar oscillations could be used as clocks, to measure kinematics, pulsar timing, stellar companions, and so on. We got results that make me optimistic. I started to write some words about doing this in practice using the Kepler data. It all hinges on the quality factor—the coherence of the oscillation modes—and the signal-to-noise ratio at which they are detected in the full data set.

Later in the day, Fadely showed early results from his machine learning applied to SDSS colors (to predict imaging morphology, in service of star–galaxy separation). He is basically going to throw everything he has at that. Mykytyn finished making more robust estimates of quasar lightcurve covariance properties and we started to discuss a hierarchical inference on the hyperparameters of those.

interatively reweighted least squares: Gaussian Process edition

Yesterday and today, Mykytyn and I worked out the generalization of iteratively reweighted least squares (a method for down-weighting outliers that is slightly, slightly better than sigma-clipping) to the case of a Gaussian Process noise model. We take the standard algorithm and replace the residual scaled by the noise with the deviation from the conditional mean prediction scaled by the sum in quadrature of the observational noise and the GP conditional variance. We implemented this simple algorithm on Mykytyn's quasar time-series data and it looks like it works beautifully. Now he is going to see if our inferences look better with the re-weighted uncertainties.

predict data, not latent variables

At Galaxy Coffee today, Cristina Garcia (MPIA) spoke about the quasar–galaxy cross-correlation at redshift four. The quasars have huge clustering amplitude, so the cross-correlation is expected to have large amplitude too. She uses a clever galaxy selection technique and finds consistency between the data and expectations. Also at Galaxy Coffee, Željko Ivezić (UW) showed that there are simple situations in which the mean of the data is not the best estimator of the location of a distribution function. In one example, not only was there a better estimator, but it improved with more data as 1/N (rather than 1/sqrt(N)). He strongly advocated using likelihood functions to generate estimators.

Fadely showed up in Heidelberg today, and we discussed star–galaxy classification improvements that could help PanSTARRS and LSST. As our loyal reader will recall, one issue with supervised methods for the problem of star–galaxy clasification is that we don't have any good sets of labels, even when we have HST data, or spectroscopy; there is a lot of label noise at the faint end, and good labels only exist for very small slices of the total population. We realized today that we could try to predict not the labels, but the things that go into making the labels, like "psf minus model" in SDSS or roundness and sharpness and so on. We vowed to give it a try. Key idea: Predict data, not latent variables!

log space, various

We decided to simplify Ben Johnson's calibration code by switching the spectroscopic likelihood function to log space, where the calibration vector is additive. This makes the noise non-Gaussian, of course, but at very high signal-to-noise (where he is working), the Gaussian approximation is not severe. At Milky Way group meeting, Wilma Trick (MPIA) talked about observational signatures of spiral structure in the kinematics of gas and stars, and Alexia Lewis (UW) talked about the star-formation histories of patches of M31 from PHAT data. Her method is very simple (fitting the distribution of main-sequence stars); it could be applied to the Hipparcos data on the local neighborhood, which would be super-cool.

understanding the behavior of complex code

Three of my conversations today were the following: Zhitai Zhang (MPIA) is working out, for a given stellar position in three-space and radial velocity, what kinds of orbits might that star be on, given a Milky Way gravitational potential model and the unknown proper motion. Wilma Trick (MPIA) is generating toy stellar position and velocity data on a toy Milky Way disk and censoring it with a toy selection function and then trying to infer the toy model parameters from the simulated data. Ben Johnson (UCSC) is burning in his simultaneous fit of star-cluster model and spectrophotometric calibration vector. In all three cases, the code is doing something unexpected, despite passing lots of local sanity checks. This is like the difference between unit testing and functional testing: Sometimes the whole system is hard to understand. Especially when it is a complex combination of physics, statistics, and code approximations. Is the puzzling behavior a problem with the code, or is it a problem with the customer (me), who just can't handle the truth.

Zhang sometimes finds that the true stellar orbit is low probability given the data, even when there are no observational errors on the four observed phase-space coordinates. We think this has something to do with the orbital phase being such that the (unobserved) transverse velocity is unusually large (or small). Trick is finding biased disk structure parameters, even when she is generating and fitting with the same model family and the data are perfect; we suspect either the toy-data generation or else some treatment of the (trivial) censoring. Johnson is finding that the cluster amplitude or mass is unconstrained by the spectral data, even when the calibration uncertainty is set to something finite. In each case, we can't quite tell whether the behavior is evidence for problems with the code or problems with our concepts (that is, is it a bug or a think-o?). All I know how to do in these cases is come up with sensible sanity checks on the code. We suggested sampling a larger range in transverse velocity for Zhang, making even simpler toy data for Trick, and looking at every stage in the likelihood calculation for Johnson. The latter lasted well into the night and I believe in the end it was a think-o not a bug.

don't calibrate your data!

Ben Johnson (UCSC) showed up in Heidelberg today and we discussed his project to obviate spectrophotometric calibration. He said, however, that we are not permitted to give the paper the same title as this blog post. However, I think it is apt, because he has a working system that suggests that inferences on spectra can proceed with a tiny bit of photometric data and no serious spectrophotometric calibration just as well as they can proceed when fully calibrated. Really the no-calibration method is better, in fact, because it models the spectrograph simultaneously with the spectral information; that is, it is makes less rigid assumptions about the hardware in general. The idea is to fit the calibration "vector" along with the parameters of interest, and take up the fiddly bits with a Gaussian Process. It seems to work. We spent part of the day moving a bit of furniture, to make the GP contribute multiplicatively (rather than additively) to the signal.

LSST exposure time

My very busy week ended with a seminar on Gaussian Processes, that I gave on the blackboard, supported with some projected slides showing Foreman-Mackey's demo plots. I tried to go slowly, but there is a lot to do in a one-hour lecture on something that can fill an entire semester's course. I got great questions from the crowd.

In the morning, Željko Ivezić (UW) and I discussed LSST cadence and exposure-time issues. He showed me his "conservation laws" slide, which shows that the exposure time flows down to all sorts of constraints on the survey. We talked about the plan to split the individual (probably 30-sec) exposures into two parts. What's the difference between 15+15 and 5+25 (and so on)? The expectations about the noise and the point-spread function are both related to the exposure time, and there are cosmic rays and moving objects. We tried to spec out some simple projects to analyze the relevant problems. Simplest: What are the implications for point-source position and flux measurements, in the limit that the sky is static? Even this question is not trivial.

Gaussian Processes for astronomers

I spent my research time today (much of it on planes and trains) working on my seminar (to be given tomorrow) about Gaussian Processes. I am relying heavily on the approach advocated by Foreman-Mackey, which is to start with weighted least squares (and the Bayesian generalization), then show how a kernel function in the covariance matrix changes the story, and then show how the Gaussian Process is latent in this solution. Not ready! But Foreman-Mackey made me a bunch of great demonstration figures.

DDD meeting, day 2

On the second day of the Moore Foundation meeting, I gave my talk (about flexible models for exoplanet populations, exoplanet transits, and exoplanet-discovering hardware calibration). After my talk, I had a great conversation with Emmanuel Candès (Stanford), who asked me very detailed questions about my prior beliefs. I realized in the conversation that I have been violating all my own rules: I have been setting my prior beliefs about hyper-parameters in the space of the hyper-parameters and not in the space of the data. That is, you can only assess the influence and consistency of the prior pdf (consistency with your actual beliefs) by flowing the prior through the probabilistic model and generating data from it. I bet if I did that for some of the problems I was showing, I would find that my priors are absurd. This is a great rule, which I often say to others but don't do myself: Always sample data from your prior (not just parameters). This is a rule for Bayes but also a rule for those of us who eschew realism! More generally, Candès's expressed the view that priors should derive from data—prior data—a view with which I agree deeply. Unfortunately, when it comes to exoplanet populations, there really aren't any prior data to speak of.

There were many excellent talks again today; again this is an incomplete set of highlights for me: Titus Brown (MSU) explained his work on developing infrastructure for biology and bioinformatics. He made a number of comments about getting customer (or user) stories right and developing with the current customer in mind. These resonated for me in my experiences of software development. He also said that his teaching and workshops and outreach are self-interested: They feed back deep and valuable information about the customer. Jeffrey Heer (UW) said similar things about his development of DataWrangler, d3.js, and other data visualization tools. (d3.js is github's fourth most popular repository!) He showed some beautiful visualizations. Heer's demo of DataWrangler simply blew away the crowd, and there were questions about it for the rest of the day.

Carl Kingsford (CMU) caused me (and others) to gasp when he said that the Sequence Read Archive of biological sequences cannot be searched by sequence. It turns out that searching for strings in enormous corpuses of strings is actually a very hard problem (who knew?). He is using a new structure called a Bloom Filter Tree, in which k-mers (length-k subsections) are stored in the nodes and the leaves contain the data sets that contain those k-mers. It is very clever and filled with all the lovely engineering issues that the Astrometry.net data structures were filled with lo so many years ago. Kingsford focuses on writing careful code, so the combination of clever data structures and well written code gets him orders of magnitude speed-ups over the competition.

Causal inference was an explicit or implicit component of many of the talks today. For example, Matthew Stephens (Chicago) is using natural genetic variations as a "randomized experiment" to infer gene expression and function. Laurel Larson (Berkeley) is looking for precursor events and predictors for abrupt ecological changes; since her work is being used to trigger interventions, she requires a causal model.

Blair Sullivan (NC State) spoke about performing inferences with provable properties on graphs. She noted that most interesting problems are NP hard on arbitrary graphs, but become easier on graphs that can be embedded (without crossing the edges) on a planar or low-genus space. This was surprising to me, but apparently the explanation is simple: Planar graphs are much more likely to have small sets of vertices that split the graph into disconnected sub-graphs. Another surprising thing to me is that "motif counting" (which I think is searching for identical subgraphs within a graph) is very hard; it can only be done exactly and in general for very small subgraphs (six-ish nodes).

The day ended with Laura Waller (Berkeley) talking about innovative imaging systems for microscopy, including light-field cameras, and then a general set of cameras that do non-degenerate illumination sequences and infer many properties beyond single-plane intensity measurements. She showed some very impressive demonstrations of light-field inferences with her systems, which are sophisticated, but built with inexpensive hardware. Her work has a lot of conceptual overlap with astronomy, in the areas of adaptive optics and imaging with non-degenerate masks.

DDD meeting, day 1

Today was the first day of a private finalist meeting for the new Moore Foundation Data Driven Discovery Individual Investigator grants. The format is a shootout of short talks and a few group activities. There were many extremely exciting talks at the meeting; here is just an incomplete smattering of highlights for me:

Ethan White (Utah State) showed amazing ecology data, with most data sources being people looking at things and counting things in the field, but then also some remote sensing data. He is using high-resolution time-series data on plants and animals to forecast changes in species and the ecosystem. He appeared to be checking his models "in the space of the data"—not in terms of reproduction of some external "truth"—which was nice.

Yaser Abu-Mostafa (Caltech) spoke about the practice and pitfalls of machine learning; in particular he is interested in new methods to thwart data "snooping" which is the name he gives to the problem "if you torture your data enough, it will confess".

Carey Priebe (JHU) opened his talk with the "Cortical Column Conjecture" which claims that the cortex is made up of many repeats of small network structures that are themselves, in some sense, computing primitives. This hypothesis is hard to test both because the graphs of neural connections in real brains are very noisy, and because inference on graphs (including finding repeated sub-graphs) is combinatorically hard.

Amit Singer (Princeton) is using micrographs of large molecules to infer three-dimensional structures. Each micrograph provides noisy data on a two-dimensional projection of each molecule; the collection of such projections provides enough information to both infer the Euler angles (three angles per molecule) and the three-dimensional structure (a very high-dimensional object). This project is very related to things LeCun, Barron, and I were talking about many years ago with galaxies.

Kim Reynolds (UT Southwestern Medical Center) is using genetic variation to determine which parts of a protein sequence are important for setting the structure and which are replaceable. She makes predictions about mutations that would not interrupt structure or function. She showed amazing structures of cellular flagella parts, and proposed that she might be able to create new kinds of flagella that would be structurally similar but different in detail.

empirical models for APOGEE spectra

I spent a chunk of the day with Melissa Ness (MPIA), fitting empirical models to APOGEE infrared spectra of stars. The idea is to do a simple linear supervised classification or regression, in which we figure out the dependence of the spectra on key stellar parameters, using a "training set" of stars with good stellar parameters. We worked in the pair-coding mode. By the end of the day we could show that we are able to identify regions of the spectrum that might serve as good metallicity indicators, relatively insensitive to temperature and log-g. The hopes for this project range from empirical metallicity index identification to label de-noising to building a full data-driven (supervised) stellar parameter pipeline. We ended our coding day pretty optimistic.

dust priors and likelihoods

Richard Hanson and Coryn Bailer-Jones (both MPIA) and I met today to talk about spatial priors and extinction modeling for Gaia. I showed them what I have on spatial priors, and we talked about the differences between using extinction measurements to predict new extinctions, using extinction measurements to predict dust densities, and so on. A key difference between the way I am thinking about it and the way Hanson and Bailer-Jones are thinking about it is that I don't want to instantiate the dust density (latent parameters) unless I have to. I would rather use the magic of the Gaussian Process to marginalize it out. We developed a set of issues for the document that I am writing on the subject. At Galaxy Coffee, Girish Kulkarni (MPIA) gave a great talk about the physics of the intergalactic medium and observational constraints from the absorption lines in quasar spectra.

quasar continuum blueward of Lyman alpha, Galactic center

If you go to the blue side of Lyman alpha, at reasonable redshifts (say 2), the continuum is not clearly visible, since the forest is dense and has a range of equivalent widths. Any study of IGM physics or radiation or clustering or ionization depends on an accurate continuum determination. What to do? Obviously, you should fit your continuum simultaneously with whatever else you are measuring, and marginalize out the posterior uncertainties on the continuum. Duh!

That said, few have attempted this. Today I had a long conversation with Hennawi, Eilers, Rorai, and KG Lee (all MPIA) about this; they are trying to constrain IGM physics with the transmission pdf, marginalizing out the continuum. We discussed the problem of sampling each quasar's continuum separately but having a universal set of IGM parameters. I advocated a limited case of Foreman-Mackey and my endless applications of importance sampling. Gibbs sampling would work too. We discussed how to deal with the fact that different quasars might disagree mightily about the IGM parameters. Failure of support can ruin your whole day. We came up with a clever hack that extends a histogram of samples to complete support in the parameters space.

In Milky Way group meeting, Ness (MPIA) showed that there appears to be an over-density of metal-poor stars in the inner one degree (projected) at the center of the Milky Way. She is using APOGEE data and her home-built metallicity indicators. We discussed whether the effect could be caused by issues with selection (either because of dust or a different explicit selection program in this center-Galaxy field). If the effect is real, it is extremely interesting. For example, even if the stars were formed there, why would they stay there?

extinction and dust, H-alpha photons

While "off the grid" for a long weekend, I spent time writing documents for Coryn Bailer-Jones (MPIA) and Dennis Zaritsky (Arizona). The former was about using spatial priors for inference of the three-dimensional dust density constrained by Gaia data. If you use a Gaussian Process spatial prior, you can perform the inference in extinction space (not dust space) and transfer extinction predictions to new points given extinction data without ever explicitly instantiating the dust density field. This is not a genius idea; it flows from the fact that any linear projection of a Gaussian pdf is itself a Gaussian pdf. The whole thing might not be computationally tractable, but at least it is philosophically possible. One issue with using a Gaussian Process here is that it puts support onto negative dust densities. I don't think that is a problem, but if it is a problem, the fixes are not cheap.

The latter document—for Zaritsky—is about finding the H-alpha photons that are coming from the outskirts of low-redshift galaxies by doing cross-correlations between SDSS spectroscopy and nearby galaxy centers. This project is designed to test or constrain some of the ideas in a talk at NYU by Juna Kollmeier a few months ago.

coffee, objectives in calibration

At MPIA Galaxy Coffee, Schmidt (UCSB) showed resolved spectroscopy of some highly magnified high-redshift galaxies to look for spatial variations of metallicity. Beautiful data! He also talked about the photon density at high redshift and the implications for reionization. McConnell (Hawaii) spoke about the need for more good black-hole mass determinations. He (inadvertently perhaps) showed that the BH-sigma (velocity dispersion) relation could easily have zero scatter, when you consider the uncertainties in both directions and the possibility that the sigma uncertainties are under-estimated. In general it is very hard to get a sigma on a sigma! Filed away for future thinking.

In the afternoon, I spoke with Wang, Schölkopf, Foreman-Mackey about Kepler calibration. We tentatively decided to make exoplanet discovery our primary objective. An objective is required, because we can set our "hyper-parameters" (our choices about model and model complexity) to optimize anything from exoplanet discovery to exoplanet characterization to stellar rotation determination. I am still worried about the flexibility of our models.