high-resolution dust maps, data-driven SED templates

Maria Kapala (Cape Town) showed up today, to discuss analyses of multi-wavelength imaging. Our idea, which was born years ago in work with Lang, is to build a dust map that constrains dust density and temperature and emission properties using all the Herschel bands, but works at the angular resolution of the best band. The usual practice is to smooth everything to the worst band!

Also had a long conversation with Boris Leistedt (NYU) about learning the templates simultaneously with the redshifts in a template-based photometric-redshift system. This is the right thing to do: It captures the causal model that is inherent in the template-based systems, but also captures the data-driven advantages of a machine-learning method. I am interested to know how accurately and at what resolution we could recover templates in realistic fake-data tests. We scoped a first paper on the subject.

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.

ABC, IMF, and habitability

While Foreman-Mackey and Barclay set off on a tangent to measure (or limit) exoplanet masses using n-body models of exoplanet systems observed by Kepler, I had a great phone call with Schaefer (CMU), Cisewski (CMU), Weller (CMU), and Lang about using Approximate Bayesian Computation (ABC) to ask questions about the universality of the high-mass initial mass function (IMF) in stellar clusters observed in the PHAT survey. The idea behind ABC is to do a kind of rejection sampling from the prior to make an approximation to posterior sampling in problems where it is possible to generate data sets from the model (and priors) but impractical or impossible to write down a likelihood function.

The reason we got this conversation started is that way back when we were writing Weisz et al on IMF inference, we realized that some of the ideas about how high-mass stars might form in molecular clouds (and thereby affect the formation of other less-massive stars) could be written down as a data-generating process but not as a computable likelihood function. That is, we had a perfect example for ABC. We didn't do anything about it from there, but maybe a project will start up on this. I think there might be quite a few places in astrophysics where we can generate data with a mechanistic model (a simulation or a semi-analytic model) but we don't have an explicit likelihood anywhere.

At the end of the day, Sarah Ballard (UW) gave a great Physics Colloquium on habitable exoplanets and asteroseismology, and how these two fields are related. They are related because you only know the properties of the exoplanet as well as you can understand the properties of the star, and asteroseismology rocks the latter. She mentioned anthropics momentarily, which reminded me that we should be thinking about this: The anthropic argument in exoplanet research is easier to formulate and think about than it is in cosmology, but figuring it out on the easier problem might help with the harder one.

WFC3, exoplanet searching

At Computer-vision-meets-astronomy today, Fadely showed us all some example HST WFC3 images, some models of the PSF, and some comparisons between model and observed stars. I had never put two-and-two together, but the PHAT project (on the periphery of which I lie) has taken some absolutely awesome WFC3 images for the purposes of array calibration: The PHAT images (being of M31) are absolutely teeming with stars. Indeed, it is impressive that the PHAT team can photometer them at all. We discussed strategies for flat-field determination given that we have a good but not perfect PSF model and a lot of heterogeneous data.

After that but before lunch, we more-or-less decided that while Foreman-Mackey works on a Kepler light-curve likelihood function paper, Angus (Oxford) should start work on a Kepler light-curve exoplanet search paper, making use of the same machinery. This is a great division of labor (I hope) and might eventually bring us close to the goal of everything we have been doing with Kepler, to wit, finding Earth-like planets on year-ish orbits around Sun-like stars. Pleased.

data analysis consulting, mixture models

Julianne Dalcanton (UW) gave a great Galaxy Coffee talk about mapping the dust in M31, by a very clever mixture-model fitting approach, fitting the extinctions towards red-giant stars in the PHAT data. She shows amazing angular resolution and incredible relationships with emission measures from infrared and millimeter. And that after a great Galaxy Coffee talk from Aaron Dutton (MPIA) summarizing the meeting "The Physical Link between Galaxies and their Halos". He did the very clever / sensible thing of choosing the three things about the meeting that most impressed him, and only talking about those.

I got back into my "data analysis guru" mode today, with long conversations with a group (including Smolcic and Groves and many others) that is trying to detect very faint sources in very deep JVLA imaging of blank fields: How do you know that the sources you are seeing are real, and how do you measure their properties? Much of the non-triviality comes from the fact that the raw data are interferometry visibilities and the maps are made with (relatively speaking) black boxes. I was a very strong advocate of jackknife (or full likelihood modeling, which is a great plan but outrageously hard given where everyone is right now).

I also spoke a bit with Watkins (MPIA) and van de Ven (MPIA) about modeling dynamical systems in the presence of a large-amplitude background, such that the model must include not just the dynamical system of interest (a stellar cluster, in this case) but also the foreground or background (the non-cluster stars of the galaxy or Galaxy, in this case). I worked through my understanding of the mixture models that the Loyal Reader (tm) knows so much about. It gets confusing when the mixture amplitudes are conditioned on observables or data that are not explicitly being modeled in the likelihood; for example in the Watkins case, the velocity distribution is a mixture of components, but the mixture amplitudes depend on position. I originally recommended modeling position and velocity simultaneously, but given the crazy selection of stars they face, it is better to model velocity conditioned on position. This makes the mixture less trivial.

In general we would all be better off if we understood mixture models much better. They obviate hard classification and capture a lot of our ideas about how our data are generated.

MCMC and model grids, stacking

For some reason these days, I keep getting asked about running MCMC in situations where the model is only defined on a discrete grid. Answer: No problem! You can either run MCMC also on that grid (with discrete proposal distributions) or else you can run MCMC in a continuous space, but snap-to-grid for the likelihood calculation (and then snap back off when you are done). Things got a bit hairier when the PHAT team (Weisz and Fouesneau are in town for a code sprint, Gordon was on the phone) were asking about the same but with non-trivial priors and exceedingly non-uniform model grids. So I decided to write down the full answer. I didn't finish by the end of the day.

It being spring break, Price-Whelan also spent a spa day down at NYU, to re-start our project on co-adding (or, really on not co-adding) imaging data. We are showing that photometric modeling (or measurement) in unstacked data beats the same in stacked data, even for sources too faint to see at any epoch. That is, you might need to stack the data in order to see the source, but you don't need to in order to detect or measure it. Worse than don't need to: You get more precision by not stacking. Duh!

MCMC initialization and convergence

Fouesneau (UW) and I discussed and adjusted the initialization for his ensemble sampling (with emcee) fits of King models to young stellar clusters in the PHAT data. Our pretty consistent experience is that you should initialize the ensemble in a pretty small ball in parameter space and then burn it in to a fair sampling. We also looked at the autocorrelation times, which are not stably measured in short chains, and only stably measured when the chains are long enough that you are properly converged. All of the experience we have developed in MCMC sampling for inference in typical astronomy problems ought to be passed on to the community somewhere! After Foreman-Mackey finishes his current exoplanet paper, we might take a couple weeks and write the how to do MCMC document.

fun day

On my way to work I ran into Fergus, taking photos in preparation for a SIGGRAPH submission. I helped him out for an hour. Fouesneau (UW) made plots to evaluate the quality of his PHAT young stellar cluster fits. They look great; I think he has nailed this fitting. Next up: Do the fits in multiple bands and check that the color distribution is tighter than it used to be. At lunch Foreman-Mackey figured out that the dry-erase glass coffee table I am building (yes, building) could include a back-projected external monitor. That lost us some time with online shopping for compact, short-throw projectors. At MCMC meeting (Goodman, Hou, Foreman-Mackey, Fadely, and I) discussed combinatoric degeneracies and their ubiquity and hardness. The fact that you can reorder a set of planets or stars or whatever in your model and leave the likelihood unchanged is a real problem: You either fix this massive, perfect degeneracy in the prior (by enforcing order) or else live with it. Either way, it hurts performance in sampling with almost all known methods. Late in the day, as I was describing my dream of a data-driven modeling approach to APOGEE chemical abundances, Fadely had a brainstorm: We could consider random subsets of the spectrum to avoid unknown or uncharacterized data or noise issues. That made Foreman-Mackey immediately say random forest, a method I never thought I would find myself using. Over the next hour I became more and more convinced that random forest could do almost exactly what I need. There were many other breakthroughs in that conversation, including desiderata for the data-driven model and ways to test it or start out.

Fouesneau, Geha

Morgan Fouesneau (UW) arrived for a few days of paper-finishing intensity. He is trying to measure the unresolved (confused) light in young star clusters observed in the PHAT data on M31. We discussed the model and I asked for some diagnostics. They look good; the models seem to be working better than Fouesneau himself thought! Marla Geha (Yale) also showed up for some spectroscopic calibration consulting. We discussed the projects I had been doing with Roweis (and Bolton a bit) at the time of Roweis's death. She may want to give them a shot! We started by specifying some visualizations to make from her current calibrations. The idea is: If we can find invariants of the calibration, we can strongly regularize the calibration fits. Roweis and I were doing this in full generality, but the first thing to do is just look at the data.

more galaxies and more dust

In between talks by KG Lee (MPIA), Bovy, and Farrar (NYU), Mykytyn and I continued to work on the Atlas measurements, and Kapala and I looked some more into the relationship between Herschel dust maps and PHAT extinction measurements on luminous stars. The latter is quite a mess and we can't figure out what parts of the mess are problems with the star fitting or real variance in the extinction, and what parts of the mess are problems with the star fitting or real correlations between star and dust properties. We discussed with Groves.

PHAT Camp, day five

I worked on wording, equations, and pseudo-code for completeness inclusion in the various PHAT projects underway. Gordon and Weisz successfully ran the PDMF fitting code from end-to-end on both simulated and real clusters. Also Maria Kapala (MPIA) came by and we looked at the pixel histogram in various PHAT UV HST images of the M31 disk. She is looking to see if it is possible to measure the light from the unresolved and undetected stars in the UV as a shift in the pixel histogram in the UV images. The histograms were nearly incomprehensible, making me wonder if MultiDrizzle is messing up the data. Of course I have never been fond of Drizzle-like algorithms, because they treat the data as overlapping square bins of photons and not samples of a pixel-convolved-psf-convolved intensity field. More on that here next week (for Brewer's sake if no-one else's). We may try to chase down the reasons for the pixel histogram issues next week.

measuring the undetectable

Brewer, Foreman-Mackey, and I have been working hard on sampling all week, so we took a break today to discuss the astrophysics of the problem. We want to measure the number–flux relation (or luminosity function) of stars in a cluster, below where confusion becomes a problem for identifying or photometering stars. There is an old literature from radio astronomy in the sixties and seventies about inferring properties of faint, overlapping sources from looking at statistics of the resulting confusion noise. But nowadays with awesome probabilistic techniques like Brewer's, we might be able to directly generatively model the confused background, and thereby produce probabilistic information about the astrophysical systems that generate it. Could be huge, in these days of confused Herschel data, the PHAT data on M31, the Galactic Center, and so on.

One of the things we discussed is how to make the problem as challenging as possible. We don't want to do any cheating, where the brighter (resolved) stars are effectively giving us the information we seek about the fainter (unresolved) stars.

uncertainties are parameters too

The day started with a call with Rory Holmes about our nearly-finished paper on self-calibration of imaging surveys. One of the things we discussed was journal choice. Not an easy one. We also agreed to move to git and cloud-based code hosting.

Later in the day, Fouesneau (UW), Weisz (UW), and I worked on issues of completeness, cluster membership posterior probabilities, and radial-profile fitting. On the latter point, Fouesneau pointed out that the errors in the measurements of the radial profiles (which are photometric) are probably under-estimated because they are generated by shot noise not in the number of photons, but in the number of stars, which have a huge dynamic range in brightness. He doesn't trust the uncertainties that are reported to him. Not having the ability to re-do the error analysis, we discussed the various things that the uncertainties could depend on, among the data outputs and model inputs we do have. Once we had that written down, we realized that we could just parameterize the dependence of the uncertainties on the measured and model quantities, and fit for them. We pair-coded that in Fouesneau's sandbox and it worked! So we might be treating uncertainties in the way they deserve. It reminded me of conversations I have had in the past with Brendon Brewer (UCSB).

young clusters

I continued working on various aspects of the young clusters in M31 with Fouesneau (UW), Gordon (STScI), and Weisz (UW). We discussed at length how the mixture model works, with a mixing that is position-dependent. We figured out, conceptually, how the cluster properties in the mixture model are constrained strongly by stars that fit the cluster well, and not by stars that don't. We agreed that we need more pedagogical stuff on mixture models, for the paper, for talks, and for the world. I worked with Fouesneau to debug some code that made it seem like emcee was not working; in the end emcee was not the problem. We had a great Korean lunch on 32nd Street before sending Gordon off at Penn Station, and then spent the walk home re-discussing completeness. We decided, tentatively, to do a naive completeness estimation and just note the limitations. Doing the Right Thing (tm) is beyond our current scope, and probably won't change our results very much (for reasons we can quantitatively argue).

completeness

Lang arrived for the day, and after a re-cap of what we are trying to do with the PHAT data, the issues of completeness came up. We had a lively discussion with Fouesneau (UW), Gordon (STScI), and Weisz (UW) of how we should measure the completeness and how we should use those measurements. We realized that there are so many subtleties, there is a paper that could be written on this alone. It is hard to measure and easy to use wrongly! Coding continued, and in the evening, Lang and I re-discussed our first papers from the Tractor. I promised to email Lang minutes of that chat.

young, PHAT clusters

A chunk of the PHAT team—Dalcanton (UW), Fouesneau (UW), Gordon (STScI), Weisz (UW)—arrived today to talk about stellar SED fitting and propagation of uncertainties therein to quantities and studies of interest. (For those of you who do not read and remember absolutely everything I write here, the PHAT project is a Dalcanton-PI six-band imaging survey over a large fraction of the M31 disk to create a catalog of tens of millions of stars and do a lot of the science you can do with that.) Today we made plans for the week-long sprint, which appear to be to use the outputs of SED fitting for every star on a big parameter grid and inputs from a model of the stellar population in a cluster, plus some crazy integration, to build a marginalized likelihood for the cluster parameters. Now everyone is coding like crazy. I think we surprised ourselves when we decided we would work in IDL and not Python. Crazy, but pragmatic.

electromagnetism and massive stars

Inspired in part by our meetings yesterday about Fergus's modeling of imaging data in a coronograph, I worked on a physically motivated re-factor of my physically motivated code to model electromagnetic fields (phase and amplitude) in astronomical telescopes and cameras. I am just a few dozen lines of code away from having a full model (highly approximate) of a simple coronograph.

In the afternoon, Selma de Mink (STScI) gave a nice seminar about extremely massive star evolution. Among many other things, she noted that there is a possibility that low-metallicity, rapidly rotating, massive stars could evolve to very hot temperatures and very high luminosities where no other kinds of stars can be. I think we can find these things in PHAT data on Andromeda; I need to email the team.