Simple Monte Carlo; a noise model

After the successes of yesterday on our custom radial-velocity sampler, currently called The Joker (but not pronounced how you might think), I put some time into writing the method section. One complex point of the sampling, which is fundamentally not Markov but instead just Simple Monte Carlo, is that if the SMC doesn't lead to many surviving samples, we either do more SMC or else switch over to a standard MCMC, initialized by the output of the SMC. That took some design thought; it capitalizes on an important point of problem structure, which is that—given a finite time window of observations—there is a finite resolution to likelihood peaks in the period direction. It remains to be seen if what we have designed will work.

Early in the day, I spoke with Andy Casey (Cambridge) about a possible noise model for the label outputs from The Cannon acting on the RAVE data. As my loyal reader knows, we consider the formal uncertainties coming from The Cannon to be under-estimates. It sounds like Casey has good evidence for a noise floor, which can be added in quadrature and make repeat visits to spectra more-or-less consistent. It's do-or-die because he needs to submit this paper today or tomorrow!

first group meeting at the SCCA

The single best thing today was the first meeting of the new joint group meeting of observational astrophysicists at the brand-spanking-new Simons Center for Computational Astrophysics. In addition to the union of Mike Blanton (NYU), Anthony Pullen (NYU), and my groups, the Director of the Center, David Spergel (SCCA) was there. We had great attendance; introductions alone took more than 90 minutes. Highlights for me included the following: A brief argument broke out about binary stars: Do we really know that both pairs of binaries have the same chemical abundances in detail? Spergel pointed out that the high-velocity stars that Keith Hawkins (Columbia) is finding could have implications for the early universe and the escape fraction of ultraviolet photons. Pullen talked about finding the SZ effect in filaments, Spergel mentioned work on identifying filaments using good machine learning by Shirley Ho (CMU), and MJ Vakili (NYU) talked about the work he has done on halo occupation, which could (in principle) take filament-environment as inputs. Hawkins also talked about a Gaia DR1 zero-day project; Adrian Price-Whelan (Columbia) and I promised to start off next week's group meeting with a visualization of the DR1 data, which at that point will be hours old!

Before and after group meeting, Price-Whelan and I worked on our binary-star (and exoplanet) sampler, building and executing experiments, and writing in the document. One amusing thing is that Megan Bedell (Chicago) gave us some (proprietary) exoplanet radial-velocity data to fit, where she finds one period but thinks there might be more. We found only her one period; confirming her results, but we found that we could get other possible periods if we drop her first data point or blow up her error bars (uncertainties).

One thing that came up in group meeting and I discussed afterwards with Hawkins is the project by Andy Casey (Cambridge), Hawkins, me, and many others to release a reanalysis of the RAVE data that overlap the Gaia DR1 T-GAS sample. Casey and Hawkins have detailed abundances for the red-giant stars, but not the main-sequence stars. I asked Hawkins why the main sequence is so much harder? He stunned me by saying that in fact he thinks the main sequence ought to be easier; we just have red-giant abundances now because people have worked harder on them. There's some medium-hanging (I won't say low-hanging!) fruit right there!

the best exoplanet spectrographs are much better than m/s

Today was the first day of the new academic year, so much of my day was obliterated by the most fun part of my job, which is teaching! That said, I still got in time for conversations with Megan Bedell, Dun Wang, and Boris Leistedt. I presented to Bedell my causal argument that instruments like HARPS are actually delivering much better than one meter per second precision and that the substantial scatter seen is because of the stars, not the instruments (and not the software pipelines). The argument is about the lack of covariance between measured radial velocities and (say) wavelength calibration parameters: Even if the noise contributing to wavelength calibration jitter is uncorrelated with the noise contributing to stellar jitter, it should show up as a covariance between calibration parameters and stellar radial velocity measurements. The argument is subtle, causal, and uncertain (because I am bad at this kind of reasoning). But if I am right, we don't need better instruments, and we don't need better pipelines. We need better models of stars!

Dun Wang and I discussed near-term and medium-term publishing plans. The top priority is to finish his paper on image differencing. I asked him to work hard on explaining how it is totally different from all other image differencing methods, because it uses the past images to learn regression coefficients, but builds the model of the present (target) image from other pixels in that image itself. That is, it is more general image modeling, really. And that's why it performs so well! Of course it requires a great data set for training.

Leistedt is using Gaussian Processes inside a physical model for galaxy observations given Doppler shifts. This is a completely flexible data-driven model, but constrained to obey the redshift physics implied by special relativity. That makes for a very powerful method for predicting galaxy colors at other redshifts, given an observation at a single (training) redshift. He can make photometric redshift predictions, make k corrections (my favorite), simulate future data, and train the photometric redshifts in survey A from a training set that exists only in survey B. All awesome! We discussed the scope of paper zero.

#AstroHackWeek, day five

Last day of #AstroHackWeek! What a week it has been, for me, anyway. So many projects. In the final wrap-up session, I was just blown away by each hack and accomplishment. You can get some feel for it at the hackpad. That site is subject to change and editing, so it isn't a static result, but you should still get a sense of the awesome.

In the wrap-up, Matt Mechtley (ASU) surprised me by showing that our molecular reconstruction does not depend extremely sensitively on the assumption that the molecules are viewed from an isotropic direction distribution! That is, he generated data with a (very strong) dipole pattern in the direction distribution, not aligned with any symmetry of the molecule, and we still reconstruct well. That bodes very well for the method.

Earlier in the day, I discussed with Dalya Baron (TAU) what directions to move in on the molecular imaging project. She decided to move the model towards the proper diffraction imaging model, with the photons generated not by the molecule directly in real space but by the squared norm of the Fourier Transform of the molecule. Because our methodology depends on having correct, analytic derivatives of a marginalized likelihood, this requires taking our derivatives through a sum of products of Fourier Transformed basis functions and their complex conjugates. Oddly, Baron wants to figure that out and code it up! (Oh, and unit tests, of course.) I'm impressed, and excited about the results. Watch this space.

#AstroHackWeek, day four

Today was the fourth day of #AstroHackWeek 2016. My goodness is everyone tired! A week of 10-to-15-hour days is a bit too much! Fortunately, the crowd was woken up by two awesome talks on optimization, one by Yu Feng (UCB) and one by Grigor Aslanyan (UCB). Feng talked about optimization in general, classifying methods, and helping us to understand the use of gradients and hessians, and trust regions. Aslanyan focused on the use of optimization to solve linear-algebra problems, and explained extremely clearly the conjugate gradient method. These talks, taken together, constituted the best live introduction to optimization I have ever seen. These were followed by Dan Foreman-Mackey (UW), who explained, in his usual extremely clear and engaging way, the basics of sampling, assessing convergence, and presentation of sampling-based results. So excited to get our pedagogical paper on MCMC out there.

I pair-coded with Adrian Price-Whelan (Princeton) our new project to sample the single-line binary star (or expolanet) system exactly using Simple Monte Carlo. Our innovation is to sample explicitly only in the non-linear parameters, and to deal with the linear parameters through exact marginalization at the rejection step, and through exact sampling at the output step. This all seems to work, and we can see the posterior pdf become less multi-moded as the number (or density) of observations increases.

Meanwhile, Dalya Baron (TAU) and Matt Mechtley (ASU) changed the toy microscopy problem to a better test problem (we had too much symmetry in our toy “molecule”). Everything works end-to-end, so in the end-of-day stand-up meeting, Baron stood up and explained our marginalized likelihood and how we are optimizing it, and Mecthley showed the data and the results. These pretty-much shocked the audience: Our images are so bad (a few photons each) and yet such a complex molecule can be inferred. We ended the day by discussing what directions to go next with this; I am a bit confused about the scope for paper one.

#AstroHackWeek, day three

Early in the morning, before #AstroHackWeek 2016 day three, I had a long phone call with Andy Casey (Cambridge) about The Cannon, RAVE, and Gaia DR1. He is trying to produce a catalog of detailed abundances from RAVE matched to Gaia T-GAS before the DR1 deadline, which is in two weeks! There are many hard parts to this project. One is to make a model of the red giants and a model of the main sequence, and somehow understand that these two models are physically consistent. Another is to get a training set of detailed abundances for main-sequence stars. We also talked about DR1 zero-day projects. I am still stumped as to what, exactly, I am going to do!

I also got a great email this morning from Megan Bedell (Chicago), demonstrating that the convection explanation is reasonable for the radial-velocity scatter she sees in her HARPS data. As my loyal reader may recall, we spent the last few months demonstrating that there is no evidence (and damn, did we search) that the HARPS pipeline is leaving radial-velocity information on the table. If it isn't, then the radial-velocity scatter must come from intrinsic stellar noise (or something much worse). What Bedell has shown is that the quantitative amplitude and granularity of stellar surface convection is sufficient to lead to meter-per-second jitter. Duh! Now we have to figure out how to deal with that. I have ideas, related to The Cannon.

The afternoon and evening of #AstroHackWeek was at GitHub headquarters in San Francisco, which is a (dangerously) fun place to hack. Jonathan Whitmore (Silicon Valley Data Science) gave a great presentation about all the crazy things you can do with a Jupyter notebook, which blew me away, and Phil Marshall (SLAC) gave a presentation to the company about how GitHub is integrated into astrophysics research these days (oh, and what features we would like).

From my (limited) perspective, the most important thing about day three was that Adrian Price-Whelan (Princeton) and I had the final realization that the importance sampling we wanted to do for the radial-velocity problem will not work. Sad, but true: We can't get a complex enough linear model that maps cleanly enough onto the Kepler problem. So we decided to switch gears today to Phil Marshall's favorite: Simple Monte Carlo. What we will do is sample the prior extremely densely, and then rejection sample to the posterior using the likelihood function. This is usually impossible! We will make it possible in this case by capitalizing on the linearity of two of our parameters: These two parameters we can analytically marginalize out at every sample in the four-dimensional space of non-linear parameters. That's our job for tomorrow.

Fail fast. That's what we are trying to do.

#AstroHackWeek, day two

Today was the second day of #AstroHackWeek 2016. Josh Bloom (UCB) spoke in the morning about machine learning, with a great set of tutorials based on the Jupyter notebook and lots of experience introducing scientists to supervised machine learning. He emphasized Random Forest, of course!

In the afternoon, Matt Mechtley (ASU) and Dalya Baron (TAU) generalized my toy one-d-angle molecular imaging problem to a three-d angle problem. Mechtley also sped up and improved the generation of the fake data, and Baron also got the inference (including the stochastic gradient descent) to work end-to-end. Exciting!

Adrian Price-Whelan (Princeton) and I capitalized on the finite time span of any radial-velocity data set to make a safe period grid for inference, and then did exact sampling at each period in a mixture-of-sinusoids model. That's awesome! But we realized that if we try to go to second order—that is, put in a sine and cosine at half the period as well as at the period—there is no way to control the amplitudes such that we maintain both linearity (which is crucial to our exact sampling) and interpretability (in terms of the true orbital parameters, which is what we really want to sample in). So we decided to implement the first-order problem only today and make a judgement about whether it works for our purposes tomorrow. We got a long way, and made some awesome plots that clearly show the multimodality of the posterior for sparse data sets. I think this could still be a useful tool, even if it doesn't do everything for all customers.

#AstroHackWeek, day one

They say you shouldn't mess with the timeline! #AstroHackWeek was so busy and full, I ended up not blogging properly during the week, and am writing these blog posts after the fact, based on telegraphic notes taken day-of. This is not uncommon here at Hogg's Research, and, for that, I apologize: Even when I write a post after the fact, I (misleadingly) date it for the day to which it corresponds (and give it a time stamp of one minute before midnight). One of the many reasons that this blog should not be seen as a precise historical document is that these after-the-fact blog posts can certainly be contaminated by present knowledge.

Today was the kick-off day for #AstroHackWeek, our now annual meeting at which participants learn about computational data analysis, and also work on their own computational data analysis projects. This year we had the meeting at the Berkeley Institute for Data Science (and partially supported it with the Moore-Sloan Data Science Environment that spans UCB, UW, and NYU). It was organized (beautifully) by Kyle Barbary (UCB) and Phil Marshall (SLAC).

In the morning today, both I and Jake VanderPlas (UW) spoke, about the basics of probabilistic inference. Then we had what Phil Marshall calls a “stand-up”, at which every participant introduced her or himself, said what it was they wanted to learn, and said what it was they knew well and could help with. They also said what they wanted to do or produce, if there was a well-defined plan.

In the stand-up and early in the hack session, Adrian Price-Whelan (Princeton) talked about joining matched sequential colormaps into diverging colormaps, with one option (or, really, style) emphasizing values near zero, and one emphasizing values far from zero. He had immediate success, and showed some nice results. One amusing thing that might bear fruit later in the week is that the author of the (currently ascendent) Viridis colormap is apparently owner of one of the BIDS desks in our vicinity this week. The conditions on a colormap are many and in tension: There are b/w printer issues, colorblindness issues, there are no-saturate-to-white and black issues, there are small-scale resolution issues, and etc.

I started (perhaps foolishly) two hacks. The first, which I started with Dalya Baron (TAU) and Matt Mechtley (ASU), was to make the demonstration I have of low-photon-rate, direct molecular imaging much more realistic. My demo, which I mention here, is an extreme toy, and there are many directions to make it less toy-like, and improve the internal engineering. I spent time with Baron and Mechtley getting them up to speed on what works and what doesn't, and what needs to change. The easiest change to make first is to go from one-dimensional angle sampling to full three-dimensional sampling in Euler angles (or, equivalently, projection matrices).

My second hack is to somehow, some way, build an MCMC sampler that can successfully and believably sample from all the modes in the multimodal posterior pdfs that we get in standard radial-velocity fitting problems (think: finding exoplanets and binary stars by measuring precise radial velocities). When the observations are sparse, the number of qualitatively different orbital solutions is large, and no sampler that I know of convincingly samples them all. Very late in the day, over coffee, Adrian Price-Whelan, Dan Foreman-Mackey (UW), and I had a very good idea: Sample exactly in a linear problem (mixture of sinusoids) that can be sampled more-or-less analytically, transform those samples into samples at the nearest points in the parameter space of the orbit-fitting problem, and then use importance sampling to get a provably correct (in the limit) sampling from the true posterior pdf. We have a plan for tomorrow!

MCMC tutorial and stellar ages

I spent my research time today working on finishing the instruction manual on MCMC that Foreman-Mackey and I are writing as part of my Data Analysis Recipes series. Our goal is to get this posted to arXiv this week. I enjoy writing but learned (once again) in the re-reading that I am not properly critical of my own prose. It really takes an outsider—or a long break—to see what needs to be fixed. The long-term preservation of science is in the hands of scientists, so the writing we do is important! Anyway, enough philosophy; this is pedagogy, not research: I am trying to make this document the most useful thing it can be. I also read and commented on a big new paper by Anna Y. Q. Ho on red giant masses and ages measured in the LAMOST project. Her paper includes an incredible map of the stellar ages as a function of position on the sky, and the different components of the Galaxy are obvious!

new space!

In a low-research day, I got my first view of the new location of the NYU Center for Data Science, in the newly renovated building at 60 Fifth Ave. The space is a mix of permanent, hoteling, and studio space for faculty, researchers, staff, and students, designed to meet very diverse needs and wants. It is cool! I also discussed briefly with Daniela Huppenkothen (NYU) the scope of her first paper on the states of GRS 1915, the black-hole source with extremely complex x-ray timing characteristics.

spectral signatures of convection; photo-zs without training data

I have spent part of the summer working with Megan Bedell (Chicago) to see if there is any evidence that radial velocity measurements with the HARPS instrument might be being affected by calibration issues or helped by taking some kind of hierarchical approach to calibration. We weren't building that hierarchical model, we were looking to see if there is evidence in the residuals for information that a hierarchical model could latch on to. We found nothing, to my surprise. I think this means that the HARPS pipelines are absolutely awesome. I think they are closed-source, so we can't do much but inspect the output.

Given this, we decided to start looking at stellar diagnostics—if it isn't the instrument calibration, then maybe it is actually the star itself: We need to ask whether we can we see spectral signatures that predict radial velocity. This is a very general causal formulation of the problem: We do not expect that a star's spectrum will vary with the phase of an exoplanet's orbit (unless it is a very hot planet!), so if anything about the spectrum predicts the radial velocity, we have something to latch on to. The idea is that we might see the spectral signature of hot up-welling or cold down-welling at the stellar surface. There is much work in this area, but I am not sure than anyone has done anything truly data driven (in the style, for example, of The Cannon). We discussed first steps towards doing that, with Bedell assigned plotting tasks, and me writing down some methodological ideas.

Over lunch, Boris Leistedt and I caught up on all the various projects we like to discuss. He has had the breakthrough that—if you build a proper generative model for galaxy imaging data—you don't need to have spectroscopic training sets, nor good galaxy spectral models, to get good photometric redshifts. The idea is that once you have multi-band photometry, you can predict the appearance of any observed galaxy as it would appear any other redshift using a flexible, non-parametric SED model that isn't tied to any physical galaxy model. The idea is that we use all of, but only, what we believe about how the redshift works, physically. Most machine-learning methods aren't required to get the redshift physics right, and most template-based models assume lots of auxilliary things about stars and stellar populations and dust. We also realized that, if done correctly, this method could subsume into itself the cross-correlation redshifts that the LSST project is excited about.

the best image differencing ever

I had the pleasure today of reading two draft papers, one by Dun Wang on our alternative to difference imaging based on our data-driven pixel-level model of the Kepler K2 data, and the other by Huanian Zhang (Arizona) on H-alpha emission from the outskirts of distant galaxies. Wang's paper shows (what I believe to be) the most precise image differences ever created. Of course we had amazing data to start with! But his method for image differencing is unusual; it doesn't require any model of either PSF nor the difference between them. It just empirically figures out what linear combinations of pixels in the target image predict each pixel in the target image, using the other images to determine these predictor combinations. It works very well and has been used to find microlensing events in the K2C9 data, but it has the disadvantage that it needs to run on a variability campaign; it can't be run on just two images.

The Zhang paper uses enormous numbers of galaxy-spectrum pairs in the SDSS spectroscopic samples to find H-alpha emission from the outskirts of (or—more precisely—angularly correlated with) nearby galaxies. He detects a signal! And it is 30 times fainter than any previous upper limit. So it is big news, I think, and has implications for the radiation environments of galaxies in the nearby Universe.

halo occupation and assembly bias

My research highlight today was a conversation with MJ Vakili about the paper he wrote this summer about halo occupation and what's known as “assembly bias”. Perhaps the most remarkable thing about contemporary cosmology is that the dark-matter-only simulations do a great job of explaining the large-scale structure in the galaxy distribution, despite the fact that we don't understand galaxy formation! The connection is a “halo occupation function” that puts galaxies into halos. It turns out that incredibly simple prescriptions work.

I have always been suspicious about halo occupation, because galaxy halos are not fundamental objects in gravity or cosmology; they are defined by a prescription, running on the output of a simulation. That is, they are just effective crutches, used for convenience. There was no reason to put any reality onto a halo (or a sub-halo or anything of the sort). Really there is just a density field! However, empirically, the halo description of the Universe has been both easy and useful.

Now that cosmology is seeking ever higher precision, work has started along the lines of asking what halo properties (mass, velocity amplitude, concentration, and so on) are relevant to the galaxies that form within them. The answer from the data seems to be that mass is the main driving factor. The community has expected a bias or occupation that depends on the time of formation of the halo (which itself relates to he halo concentration parameter). Vakili has been testing this, and the main punchline is that if the effect is there, it is a small one! It is a great result and he is nearly ready to submit.

My question is: Can we step out of the halo box and consider all the ways we might put galaxies into the dark-matter field? Could the data tell us what is most relevant?

dynamics of M33

In a day short of research (because: getting ready to teach again!), I spent some time working with the Simons Foundation to prepare for the #GaiaSprint, which is coming up in 8 weeks. After that I had lunch with Ekta Patel (Arizona), who has been working on the dynamics of the Local Group, and especially understanding the orbits of M31 and M33.

probabilistic redshifts

In the morning I had a long and overdue conversation with Alex Malz, who is attempting to determine galaxy one-point statistics given probabilistic photometric redshift information. That is, each galaxy (as in, say, the LSST plan and some SDSS outputs) is given a posterior probability over redshifts rather than a strict redshift determination. How are these responsibly used? It turns out that the answer is not trivial: They have to be incorporated into a hierarchical inference, in which the (often implicit) interim priors used to make the p(z) outputs is replaced by a model for the distribution of galaxies. That requires (a) mathematics of probability, and (b) knowing the interim priors. One big piece of advice or warning we have for current and future surveys is: Don't produce probabilistic redshifts unless you can produce the exact priors too! Some photometric redshift schemes don't even really know what their priors are, and this is death.

In the afternoon, I discussed various projects with John Moustakas (Siena), around Gaia and large galaxies. He mentioned that he is creating a diameter limited catalog and atlas of galaxies. I am very interested in this, but we had to part ways before discussing further.