Demographic studies of cosmic populations must contend with measurement errors and selection effects. We survey some of the key ideas astronomers have developed to deal with these complications, in the context of galaxy surveys and the literature on corrections for Malmquist and Eddington bias. From the perspective of modern statistics, such corrections arise naturally in the context of multilevel models, particularly in Bayesian treatments of such models: hierarchical Bayesian models. We survey some key lessons from hierarchical Bayesian modeling, including shrinkage estimation, which is closely related to traditional corrections devised by astronomers. We describe a framework for hierarchical Bayesian modeling of cosmic populations, tailored to features of astronomical surveys that are not typical of surveys in other disciplines. This thinned latent marked point process framework accounts for the tie between selection (detection) and measurement in astronomical surveys, treating selection and measurement error effects in a self-consistent manner.

The field of astronomy has arrived at a turning point in terms of size and complexity of both datasets and scientific collaboration. Commensurately, algorithms and statistical models have begun to adapt --- e.g., via the onset of artificial intelligence --- which itself presents new challenges and opportunities for growth. This white paper aims to offer guidance and ideas for how we can evolve our technical and collaborative frameworks to promote efficient algorithmic development and take advantage of opportunities for scientific discovery in the petabyte era. We discuss challenges for discovery in large and complex data sets; challenges and requirements for the next stage of development of statistical methodologies and algorithmic tool sets; how we might change our paradigms of collaboration and education; and the ethical implications of scientists' contributions to widely applicable algorithms and computational modeling. We start with six distinct recommendations that are supported by the commentary following them. This white paper is related to a larger corpus of effort that has taken place within and around the Petabytes to Science Workshops (https://petabytestoscience.github.io/).

This Astro2020 State of the Profession Consideration White Paper highlights the growth of astrostatistics and astroinformatics in astronomy, identifies key issues hampering the maturation of these new subfields, and makes recommendations for structural improvements at different levels that, if acted upon, will make significant positive impacts across astronomy.

Through the lens of the LSST Science Collaborations' experience, this paper advocates for new and improved ways to fund large, complex collaborations at the interface of data science and astrophysics as they work in preparation for and on peta-scale, complex surveys, of which LSST is a prime example. We advocate for the establishment of programs to support both research and infrastructure development that enables innovative collaborative research on such scales.

Over the past century, major advances in astronomy and astrophysics have been largely driven by improvements in instrumentation and data collection. With the amassing of high quality data from new telescopes, and especially with the advent of deep and large astronomical surveys, it is becoming clear that future advances will also rely heavily on how those data are analyzed and interpreted. New methodologies derived from advances in statistics, computer science, and machine learning are beginning to be employed in sophisticated investigations that are not only bringing forth new discoveries, but are placing them on a solid footing. Progress in wide-field sky surveys, interferometric imaging, precision cosmology, exoplanet detection and characterization, and many subfields of stellar, Galactic and extragalactic astronomy, has resulted in complex data analysis challenges that must be solved to perform scientific inference. Research in astrostatistics and astroinformatics will be necessary to develop the state-of-the-art methodology needed in astronomy. Overcoming these challenges requires dedicated, interdisciplinary research. We recommend: (1) increasing funding for interdisciplinary projects in astrostatistics and astroinformatics; (2) dedicating space and time at conferences for interdisciplinary research and promotion; (3) developing sustainable funding for long-term astrostatisics appointments; and (4) funding infrastructure development for data archives and archive support, state-of-the-art algorithms, and efficient computing.

We describe work in progress by a collaboration of astronomers and statisticians developing a suite of Bayesian data analysis tools for extrasolar planet (exoplanet) detection, planetary orbit estimation, and adaptive scheduling of observations. Our work addresses analysis of stellar reflex motion data, where a planet is detected by observing the "wobble" of its host star as it responds to the gravitational tug of the orbiting planet. Newtonian mechanics specifies an analytical model for the resulting time series, but it is strongly nonlinear, yielding complex, multimodal likelihood functions; it is even more complex when multiple planets are present. The parameter spaces range in size from few-dimensional to dozens of dimensions, depending on the number of planets in the system, and the type of motion measured (line-of-sight velocity, or position on the sky). Since orbits are periodic, Bayesian generalizations of periodogram methods facilitate the analysis. This relies on the model being linearly separable, enabling partial analytical marginalization, reducing the dimension of the parameter space. Subsequent analysis uses adaptive Markov chain Monte Carlo methods and adaptive importance sampling to perform the integrals required for both inference (planet detection and orbit measurement), and information-maximizing sequential design (for adaptive scheduling of observations). We present an overview of our current techniques and highlight directions being explored by ongoing research.

The radial velocity method is one of the most successful techniques for detecting exoplanets. It works by detecting the velocity of a host star induced by the gravitational effect of an orbiting planet, specifically the velocity along our line of sight, which is called the radial velocity of the star. As astronomical instrumentation has improved, radial velocity surveys have become sensitive to low-mass planets that cause their host star to move with radial velocities of 1 m/s or less. While analysis of a time series of stellar spectra can in theory reveal such small radial velocities, in practice intrinsic stellar variability (e.g., star spots, convective motion, pulsations) affects the spectra and often mimics a radial velocity signal. This signal contamination makes it difficult to reliably detect low mass planets and planets orbiting magnetically active stars. A principled approach to recovering planet radial velocity signals in the presence of stellar activity was proposed by Rajpaul et al. (2015) and involves the use of a multivariate Gaussian process model to jointly capture time series of the apparent radial velocity and multiple indicators of stellar activity. We build on this work in two ways: (i) we propose using dimension reduction techniques to construct more informative stellar activity indicators that make use of a larger portion of the stellar spectrum; (ii) we extend the Rajpaul et al. (2015) model to a larger class of models and use a model comparison procedure to select the best model for the particular stellar activity indicators at hand. By combining our high-information stellar activity indicators, Gaussian process models, and model selection procedure, we achieve substantially improved planet detection power compared to previous state-of-the-art approaches.

We present {\em Spitzer} secondary-eclipse observations of the hot Jupiter HAT-P-13 b in the 3.6 {\micron} and 4.5 {\micron} bands. HAT-P-13 b inhabits a two-planet system with a configuration that enables constraints on the planet's second Love number, \math{k\sb{2}}, from precise eccentricity measurements, which in turn constrains models of the planet's interior structure. We exploit the direct measurements of \math{e \cos \omega} from our secondary-eclipse data and combine them with previously published radial velocity data to generate a refined model of the planet's orbit and thus an improved estimate on the possible interval for \math{k\sb{2}}. We report eclipse phases of \math{0.49154 \pm 0.00080} and \math{0.49711 \pm 0.00083} and corresponding \math{e \cos \omega} estimates of \math{-0.0136 \pm 0.0013} and \math{-0.0048 \pm 0.0013}. Under the assumptions of previous work, our estimate of \math{k\sb{2}} of 0.81 {\pm} 0.10 is consistent with the lower extremes of possible core masses found by previous models, including models with no solid core. This anomalous result challenges both interior models and the dynamical assumptions that enable them, including the essential assumption of apsidal alignment. We also report eclipse depths of 0.081\% {\pm} 0.008\% in the 3.6 {\micron} channel and 0.088 \% {\pm} 0.028 \% in the 4.5 {\micron} channel. These photometric results are non-uniquely consistent with solar-abundance composition without any thermal inversion.

Observational astronomy in the time-domain era faces several new challenges. One of them is the efficient use of observations obtained at multiple epochs. The work presented here addresses faint object detection with multi-epoch data, and describes an incremental strategy for separating real objects from artifacts in ongoing surveys, in situations where the single-epoch data are summaries of the full image data, such as single-epoch catalogs of flux and direction estimates for candidate sources. The basic idea is to produce low-threshold single-epoch catalogs, and use a probabilistic approach to accumulate catalog information across epochs; this is in contrast to more conventional strategies based on co-added or stacked image data across all epochs. We adopt a Bayesian approach, addressing object detection by calculating the marginal likelihoods for hypotheses asserting there is no object, or one object, in a small image patch containing at most one cataloged source at each epoch. The object-present hypothesis interprets the sources in a patch at different epochs as arising from a genuine object; the no-object (noise) hypothesis interprets candidate sources as spurious, arising from noise peaks. We study the detection probability for constant-flux objects in a simplified Gaussian noise setting, comparing results based on single exposures and stacked exposures to results based on a series of single-epoch catalog summaries. Computing the detection probability based on catalog data amounts to generalized cross-matching: it is the product of a factor accounting for matching of the estimated fluxes of candidate sources, and a factor accounting for matching of their estimated directions. We find that probabilistic fusion of multi-epoch catalog information can detect sources with only modest sacrifice in sensitivity and selectivity compared to stacking.

Modern astronomy increasingly relies upon systematic surveys, whose dedicated telescopes continuously observe the sky across varied wavelength ranges of the electromagnetic spectrum; some surveys also observe non-electromagnetic "messengers," such as high-energy particles or gravitational waves. Stars and galaxies look different through the eyes of different instruments, and their independent measurements have to be carefully combined to provide a complete, sound picture of the multicolor and eventful universe. The association of an object's independent detections is, however, a difficult problem scientifically, computationally, and statistically, raising varied challenges across diverse astronomical applications. The fundamental problem is finding records in survey databases with directions that match to within the direction uncertainties. Such astronomical versions of the record linkage problem are known by various terms in astronomy: cross-matching, cross-identification, and directional, positional, or spatio-temporal coincidence assessment. Astronomers have developed several statistical approaches for such problems, largely independently of related developments in other disciplines. Here we review emerging approaches that compute (Bayesian) probabilities for the hypotheses of interest: possible associations, or demographic properties of a cosmic population that depend on identifying associations. Many cross-identification tasks can be formulated within a hierarchical Bayesian partition model framework, with components that explicitly account for astrophysical effects (e.g., source brightness vs. wavelength, source motion, or source extent), selection effects, and measurement error. We survey recent developments, and highlight important open areas for future research.

Time-correlated noise is a significant source of uncertainty when modeling exoplanet light-curve data. A correct assessment of correlated noise is fundamental to determine the true statistical significance of our findings. Here we review three of the most widely used correlated-noise estimators in the exoplanet field, the time-averaging, residual-permutation, and wavelet-likelihood methods. We argue that the residual-permutation method is unsound in estimating the uncertainty of parameter estimates. We thus recommend to refrain from this method altogether. We characterize the behavior of the time averaging's rms-vs.-bin-size curves at bin sizes similar to the total observation duration, which may lead to underestimated uncertainties. For the wavelet-likelihood method, we note errors in the published equations and provide a list of corrections. We further assess the performance of these techniques by injecting and retrieving eclipse signals into synthetic and real Spitzer light curves, analyzing the results in terms of the relative-accuracy and coverage-fraction statistics. Both the time-averaging and wavelet-likelihood methods significantly improve the estimate of the eclipse depth over a white-noise analysis (a Markov-chain Monte Carlo exploration assuming uncorrelated noise). However, the corrections are not perfect, when retrieving the eclipse depth from Spitzer datasets, these methods covered the true (injected) depth within the 68\% credible region in only $\sim$45--65\% of the trials. Lastly, we present our open-source model-fitting tool, Multi-Core Markov-Chain Monte Carlo ({MC$^3$}). This package uses Bayesian statistics to estimate the best-fitting values and the credible regions for the parameters for a (user-provided) model. {MC$^3$} is a Python/C code, available at https://github.com/pcubillos/MCcubed.

We present a method to simultaneously model the dust far-infrared spectral energy distribution (SED) and the total infrared $-$ carbon monoxide (CO) integrated intensity $(S_{\rm IR}-I_{\rm CO})$ relationship. The modelling employs a hierarchical Bayesian (HB) technique to estimate the dust surface density, temperature ($T_{\rm eff}$), and spectral index at each pixel from the observed far-infrared (FIR) maps. Additionally, given the corresponding CO map, the method simultaneously estimates the slope and intercept between the FIR and CO intensities, which are global properties of the observed source. The model accounts for correlated and uncorrelated uncertainties, such as those present in Herschel observations. Using synthetic datasets, we demonstrate the accuracy of the HB method, and contrast the results with common non-hierarchical fitting methods. As an initial application, we model the dust and gas on 100 pc scales in the Magellanic Clouds from Herschel FIR and NANTEN CO observations. The slopes of the $\log S_{\rm IR}-\log I_{\rm CO}$ relationship are similar in both galaxies, falling in the range 1.1$-$1.7. However, in the SMC the intercept is nearly 3 times higher, which can be explained by its lower metallicity than the LMC, resulting in a larger $S_{\rm IR}$ per unit $I_{\rm CO}$. The HB modelling evidences an increase in $T_{\rm eff}$ in regions with the highest $I_{\rm CO}$ in the LMC. This may be due to enhanced dust heating in the densest molecular regions from young stars. Such simultaneous dust and gas modelling may reveal variations in the properties of the ISM and its association with other galactic characteristics, such as star formation rates and/or metallicities.

The Second Workshop on Extreme Precision Radial Velocities defined circa 2015 the state of the art Doppler precision and identified the critical path challenges for reaching 10 cm/s measurement precision. The presentations and discussion of key issues for instrumentation and data analysis and the workshop recommendations for achieving this precision are summarized here. Beginning with the HARPS spectrograph, technological advances for precision radial velocity measurements have focused on building extremely stable instruments. To reach still higher precision, future spectrometers will need to produce even higher fidelity spectra. This should be possible with improved environmental control, greater stability in the illumination of the spectrometer optics, better detectors, more precise wavelength calibration, and broader bandwidth spectra. Key data analysis challenges for the precision radial velocity community include distinguishing center of mass Keplerian motion from photospheric velocities, and the proper treatment of telluric contamination. Success here is coupled to the instrument design, but also requires the implementation of robust statistical and modeling techniques. Center of mass velocities produce Doppler shifts that affect every line identically, while photospheric velocities produce line profile asymmetries with wavelength and temporal dependencies that are different from Keplerian signals. Exoplanets are an important subfield of astronomy and there has been an impressive rate of discovery over the past two decades. Higher precision radial velocity measurements are required to serve as a discovery technique for potentially habitable worlds and to characterize detections from transit missions. The future of exoplanet science has very different trajectories depending on the precision that can ultimately be achieved with Doppler measurements.

Several planetary satellites apparently have subsurface seas that are of great interest for, among other reasons, their possible habitability. The geologically diverse Saturnian satellite Enceladus vigorously vents liquid water and vapor from fractures within a south polar depression and thus must have a liquid reservoir or active melting. However, the extent and location of any subsurface liquid region is not directly observable. We use measurements of control points across the surface of Enceladus accumulated over seven years of spacecraft observations to determine the satellite's precise rotation state, finding a forced physical libration of 0.120 $\pm$ 0.014{\deg} (2{\sigma}). This value is too large to be consistent with Enceladus's core being rigidly connected to its surface, and thus implies the presence of a global ocean rather than a localized polar sea. The maintenance of a global ocean within Enceladus is problematic according to many thermal models and so may constrain satellite properties or require a surprisingly dissipative Saturn.

Ultra-high energy cosmic rays (UHECRs) are atomic nuclei with energies over ten million times energies accessible to human-made particle accelerators. Evidence suggests that they originate from relatively nearby extragalactic sources, but the nature of the sources is unknown. We develop a multilevel Bayesian framework for assessing association of UHECRs and candidate source populations, and Markov chain Monte Carlo algorithms for estimating model parameters and comparing models by computing, via Chib's method, marginal likelihoods and Bayes factors. We demonstrate the framework by analyzing measurements of 69 UHECRs observed by the Pierre Auger Observatory (PAO) from 2004-2009, using a volume-complete catalog of 17 local active galactic nuclei (AGN) out to 15 megaparsecs as candidate sources. An early portion of the data ("period 1," with 14 events) was used by PAO to set an energy cut maximizing the anisotropy in period 1; the 69 measurements include this "tuned" subset, and subsequent "untuned" events with energies above the same cutoff. Also, measurement errors are approximately summarized. These factors are problematic for independent analyses of PAO data. Within the context of "standard candle" source models (i.e., with a common isotropic emission rate), and considering only the 55 untuned events, there is no significant evidence favoring association of UHECRs with local AGN vs. an isotropic background. The highest-probability associations are with the two nearest, adjacent AGN, Centaurus A and NGC 4945. If the association model is adopted, the fraction of UHECRs that may be associated is likely nonzero but is well below 50%. Our framework enables estimation of the angular scale for deflection of cosmic rays by cosmic magnetic fields; relatively modest scales of $\approx\!3^{\circ}$ to $30^{\circ}$ are favored. Models that assign a large fraction of UHECRs to a single nearby source (e.g., Centaurus A) are ruled out unless very large deflection scales are specified a priori, and even then they are disfavored. However, including the period 1 data alters the conclusions significantly, and a simulation study supports the idea that the period 1 data are anomalous, presumably due to the tuning. Accurate and optimal analysis of future data will likely require more complete disclosure of the data.

Exoplanet WASP-14b is a highly irradiated, transiting hot Jupiter. Joshi et al. calculate an equilibrium temperature Teq of 1866 K for zero albedo and reemission from the entire planet, a mass of 7.3 +/- 0.5 Jupiter masses and a radius of 1.28 +/- 0.08 Jupiter radii. Its mean density of 4.6 g/cm3 is one of the highest known for planets with periods less than 3 days. We obtained three secondary eclipse light curves with the Spitzer Space Telescope. The eclipse depths from the best jointly fit model are $0.224\%$ +/- $0.018\%$ at 4.5 {\mu}m and $0.181\%$ +/- $0.022\%$ at 8.0 {\mu}m. The corresponding brightness temperatures are 2212 +/- 94 K and 1590 +/- 116 K. A slight ambiguity between systematic models suggests a conservative 3.6 {\mu}m eclipse depth of $0.19\%$ +/- $0.01\%$ and brightness temperature of 2242 +/- 55 K. Although extremely irradiated, WASP-14b does not show any distinct evidence of a thermal inversion. In addition, the present data nominally favor models with day night energy redistribution less than $~30\%$. The current data are generally consistent with oxygen-rich as well as carbon-rich compositions, although an oxygen-rich composition provides a marginally better fit. We confirm a significant eccentricity of e = 0.087 +/- 0.002 and refine other orbital parameters.

A preliminary study of a set of well-isolated pulses in GRB light curves indicates that simple pulse models, with smooth and monotonic pulse rise and decay regions, are inadequate. Examining the residuals of fits of pulses to such models suggests the following patterns of departure from the smooth pulse model of Norris et al. (2005): A Precursor Shelf occurs prior to or concurrent with the exponential Rapid Rise. The pulse reaches maximum intensity at the Peak Plateau, then undergoes a Rapid Decay. The decay changes into an Extended Tail. Pulses are almost universally characterized by hard-to-soft evolution, arguing that the new pulse features reflect a single physical phenomenon, rather than artifacts of pulse overlap.

This perspective chapter briefly surveys: (1) past growth in the use of Bayesian methods in astrophysics; (2) current misconceptions about both frequentist and Bayesian statistical inference that hinder wider adoption of Bayesian methods by astronomers; and (3) multilevel (hierarchical) Bayesian modeling as a major future direction for research in Bayesian astrostatistics, exemplified in part by presentations at the first ISI invited session on astrostatistics, commemorated in this volume. It closes with an intentionally provocative recommendation for astronomical survey data reporting, motivated by the multilevel Bayesian perspective on modeling cosmic populations: that astronomers cease producing catalogs of estimated fluxes and other source properties from surveys. Instead, summaries of likelihood functions (or marginal likelihood functions) for source properties should be reported (not posterior probability density functions), including nontrivial summaries (not simply upper limits) for candidate objects that do not pass traditional detection thresholds.

This paper summarizes a presentation for a panel discussion on "The Future of Astrostatistics" held at the Statistical Challenges in Modern Astronomy V conference at Pennsylvania State University in June 2011. I argue that the emerging needs of astrostatistics may both motivate and benefit from fundamental developments in statistics. I highlight some recent work within statistics on fundamental topics relevant to astrostatistical practice, including the Bayesian/frequentist debate (and ideas for a synthesis), multilevel models, and multiple testing. As an important direction for future work in statistics, I emphasize that astronomers need a statistical framework that explicitly supports unfolding chains of discovery, with acquisition, cataloging, and modeling of data not seen as isolated tasks, but rather as parts of an ongoing, integrated sequence of analyses, with information and uncertainty propagating forward and backward through the chain. A prototypical example is surveying of astronomical populations, where source detection, demographic modeling, and the design of survey instruments and strategies all interact.

This paper is an invited commentary on Tamas Budavari's presentation, "On statistical cross-identification in astronomy," for the Statistical Challenges in Modern Astronomy V conference held at Pennsylvania State University in June 2011. I begin with a brief review of previous work on probabilistic (Bayesian) assessment of directional and spatio-temporal coincidences in astronomy (e.g., cross-matching or cross-identification of objects across multiple catalogs). Then I discuss an open issue in the recent innovative work of Budavari and his colleagues on large-scale probabilistic cross-identification: how to assign prior probabilities that play an important role in the analysis. With a simple toy problem, I show how Bayesian multilevel modeling (hierarchical Bayes) provides a principled framework that justifies and generalizes pragmatic rules of thumb that have been successfully used by Budavari's team to assign priors.

The Earth is continuously showered by charged cosmic ray particles, naturally produced atomic nuclei moving with velocity close to the speed of light. Among these are ultra high energy cosmic ray particles with energy exceeding 5x10^19 eV, which is ten million times more energetic than the most energetic particles produced at the Large Hadron Collider. Astrophysical questions include: what phenomenon accelerates particles to such high energies, and what sort of nuclei are energized? Also, the magnetic deflection of the trajectories of the cosmic rays makes them potential probes of galactic and intergalactic magnetic fields. We develop a Bayesian hierarchical model that can be used to compare different association models between the cosmic rays and source population, using Bayes factors. A measurement model with directional uncertainties and accounting for non-uniform sky exposure is incoporated into the model. The methodology allows us to learn about astrophysical parameters, such as those governing the source luminosity function and the cosmic magnetic field.

The dayside of HD 149026b is near the edge of detectability by the Spitzer Space Telescope. We report on eleven secondary-eclipse events at 3.6, 4.5, 3 x 5.8, 4 x 8.0, and 2 x 16 microns plus three primary-transit events at 8.0 microns. The eclipse depths from jointly-fit models at each wavelength are 0.040 +/- 0.003% at 3.6 microns, 0.034 +/- 0.006% at 4.5 microns, 0.044 +/- 0.010% at 5.8 microns, 0.052 +/- 0.006% at 8.0 microns, and 0.085 +/- 0.032% at 16 microns. Multiple observations at the longer wavelengths improved eclipse-depth signal-to-noise ratios by up to a factor of two and improved estimates of the planet-to-star radius ratio (Rp/Rs = 0.0518 +/- 0.0006). We also identify no significant deviations from a circular orbit and, using this model, report an improved period of 2.8758916 +/- 0.0000014 days. Chemical-equilibrium models find no indication of a temperature inversion in the dayside atmosphere of HD 149026b. Our best-fit model favors large amounts of CO and CO2, moderate heat redistribution (f=0.5), and a strongly enhanced metallicity. These analyses use BiLinearly-Interpolated Subpixel Sensitivity (BLISS) mapping, a new technique to model two position-dependent systematics (intrapixel variability and pixelation) by mapping the pixel surface at high resolution. BLISS mapping outperforms previous methods in both speed and goodness of fit. We also present an orthogonalization technique for linearly-correlated parameters that accelerates the convergence of Markov chains that employ the Metropolis random walk sampler. The electronic supplement contains light-curve files and supplementary figures.

I describe ongoing work developing Bayesian methods for flexible modeling of arrival time series data without binning, aiming to improve detection and measurement of X-ray and gamma-ray pulsars, and of pulses in gamma-ray bursts. The methods use parametric and semiparametric Poisson point process models for the event rate, and by design have close connections to conventional frequentist methods currently used in time-domain astronomy.

The transiting exoplanet WASP-18b was discovered in 2008 by the Wide Angle Search for Planets (WASP) project. The Spitzer Exoplanet Target of Opportunity Program observed secondary eclipses of WASP-18b using Spitzer's Infrared Array Camera (IRAC) in the 3.6 micron and 5.8 micron bands on 2008 December 20, and in the 4.5 micron and 8.0 micron bands on 2008 December 24. We report eclipse depths of 0.30 +/- 0.02%, 0.39 +/- 0.02%, 0.37 +/- 0.03%, 0.41 +/- 0.02%, and brightness temperatures of 3100 +/- 90, 3310 +/- 130, 3080 +/- 140 and 3120 +/- 110 K in order of increasing wavelength. WASP-18b is one of the hottest planets yet discovered - as hot as an M-class star. The planet's pressure-temperature profile most likely features a thermal inversion. The observations also require WASP-18b to have near-zero albedo and almost no redistribution of energy from the day-side to the night side of the planet.

I describe ongoing work on development of Bayesian methods for exploring periodically varying phenomena in astronomy, addressing two classes of sources: pulsars, and extrasolar planets (exoplanets). For pulsars, the methods aim to detect and measure periodically varying signals in data consisting of photon arrival times, modeled as non-homogeneous Poisson point processes. For exoplanets, the methods address detection and estimation of planetary orbits using observations of the reflex motion "wobble" of a host star, including adaptive scheduling of observations to optimize inferences.

We observed two secondary eclipses of the exoplanet WASP-12b using the Infrared Array Camera on the Spitzer Space Telescope. The close proximity of WASP-12b to its G-type star results in extreme tidal forces capable of inducing apsidal precession with a period as short as a few decades. This precession would be measurable if the orbit had a significant eccentricity, leading to an estimate of the tidal Love number and an assessment of the degree of central concentration in the planetary interior. An initial ground-based secondary eclipse phase reported by Lopez-Morales et al. (0.510 +/- 0.002) implied eccentricity at the 4.5 sigma level. The spectroscopic orbit of Hebb et al. has eccentricity 0.049 +/- 0.015, a 3 sigma result, implying an eclipse phase of 0.509 +/- 0.007. However, there is a well documented tendency of spectroscopic data to overestimate small eccentricities. Our eclipse phases are 0.5010 +/- 0.0006 (3.6 and 5.8 microns) and 0.5006 +/- 0.0007 (4.5 and 8.0 microns). An unlikely orbital precession scenario invoking an alignment of the orbit during the Spitzer observations could have explained this apparent discrepancy, but the final eclipse phase of Lopez-Morales et al. (0.510 -0.006 / +0.007) is consistent with a circular orbit at better than 2 sigma. An orbit fit to all the available transit, eclipse, and radial-velocity data indicates precession at <1 sigma; a non-precessing solution fits better. We also comment on analysis and reporting for Spitzer exoplanet data in light of recent re-analyses.

I describe a framework for adaptive scientific exploration based on iterating an Observation--Inference--Design cycle that allows adjustment of hypotheses and observing protocols in response to the results of observation on-the-fly, as data are gathered. The framework uses a unified Bayesian methodology for the inference and design stages: Bayesian inference to quantify what we have learned from the available data and predict future data, and Bayesian decision theory to identify which new observations would teach us the most. When the goal of the experiment is simply to make inferences, the framework identifies a computationally efficient iterative ``maximum entropy sampling'' strategy as the optimal strategy in settings where the noise statistics are independent of signal properties. Results of applying the method to two ``toy'' problems with simulated data--measuring the orbit of an extrasolar planet, and locating a hidden one-dimensional object--show the approach can significantly improve observational efficiency in settings that have well-defined nonlinear models. I conclude with a list of open issues that must be addressed to make Bayesian adaptive exploration a practical and reliable tool for optimizing scientific exploration.

I discuss an issue arising in analyzing data from astronomical surveys: accounting for measurement uncertainties in the properties of individual sources detected in a survey when making inferences about the entire population of sources. Source uncertainties require the analyst to introduce unknown ``incidental'' parameters for each source. The number of parameters thus grows with the size of the sample, and standard theorems guaranteeing asymptotic convergence of maximum likelihood estimates fail in such settings. From the Bayesian point of view, the missing ingredient in such analyses is accounting for the volume in the incidental parameter space via marginalization. I use simple simulations, motivated by modeling the distribution of trans-Neptunian objects surveyed in the outer solar system, to study the effects of source uncertainties on inferences. The simulations show that current non-Bayesian methods for handling source uncertainties (ignoring them, or using an ad hoc incidental parameter integration) produce incorrect inferences, with errors that grow more severe with increasing sample size. In contrast, accounting for source uncertainty via marginalization leads to sound inferences for any sample size.

We analyze the Time-Tagged Event (TTE) data from observations of gamma ray bursts (GRBs) and soft gamma repeaters (SGRs) by the Burst and Transient Source Experiment (BATSE). These data provide the best available time resolution for GRBs and SGRs. We have performed an extensive search for weak periodic signals in the frequency range 400 Hz to 2500 Hz using the burst records for 2203 GRBs and 152 SGR flares. The study employs the Rayleigh power as a test statistic to evaluate the evidence for periodic emissions. We find no evidence of periodic emissions from these events at these frequencies. In all but a very few cases the maximum power values obtained are consistent with what would be expected by chance from a non-periodic signal. In those few instances where there is marginal evidence for periodicity there are problems with the data that cast doubt on the reality of the signal. For classical GRBs, the largest Rayleigh power occurs in bursts whose TTE data appear to be corrupted. For SGRs, our largest Rayleigh power, with a significance of 1%, occurs in one record for SGR 1900+14 (at 2497 Hz), and in no other outbursts associated with this source; we thus consider it unlikely to represent detection of a real periodicity. From simulations, we deduce that the Rayleigh test would have detected significant oscillations with relative amplitude ~10% about half the time. Thus, we conclude that high frequency oscillations, if present, must have small relative amplitudes.

We present a Bayesian analysis of the energies and arrival times of the neutrinos from supernova SN 1987A detected by the Kamiokande II, IMB, and Baksan detectors, and find strong evidence for two components in the neutrino signal: a long time scale component from thermal Kelvin-Helmholtz cooling of the nascent neutron star, and a brief (~< 1 s), softer component similar to that expected from emission by accreting material in the delayed supernova scenario. In the context of this model, we show that the data constrain the electron antineutrino rest mass to be less than 5.7 eV with 95% probability. Our analysis takes advantage of significant advances that have occured in the years since the detections in both our understanding of the supernova mechanism and our ability to analyze sparse data. As a result there are substantial differences between our inferences and those found in earlier studies. We find that two-component models for the neutrino signal make the data >100 times more probable than single-component models. In addition, the radius and binding energy of the nascent neutron star implied by single-component models deviates significantly from the values predicted by current neutron star models, whereas those implied by models with an accretion component are in complete agreement with the predictions. As a result, two-component models are hundreds to thousands of times more probable than single-component models. The neutrino data thus provide the first direct observational evidence in favor of the delayed supernova scenario over the prompt scenario. (Abridged abstract)

We explore the possible role of evolution in the analysis of data on SNe Ia at cosmological distances. First, using a variety of simple sleuthing techniques, we find evidence that the properties of the high and low redshift SNe Ia observed so far differ from one another. Next, we examine the effects of including simple phenomenological models for evolution in the analysis. The result is that cosmological models and evolution are highly degenerate with one another, so that the incorporation of even very simple models for evolution makes it virtually impossible to pin down the values of $\Omega_M$ and $\Omega_\Lambda$, the density parameters for nonrelativistic matter and for the cosmological constant, respectively. Moreover, we show that if SNe Ia evolve with time, but evolution is neglected in analyzing data, then, given enough SNe Ia, the analysis hones in on values of $\Omega_M$ and $\Omega_\Lambda$ which are incorrect. Using Bayesian methods, we show that the probability that the cosmological constant is nonzero (rather than zero) is unchanged by the SNe Ia data when one accounts for the possibility of evolution, provided that we do not discriminate among open, closed and flat cosmologies a priori. The case for nonzero cosmological constant is stronger if the Universe is presumed to be flat, but still depends sensitively on the degree to which the peak luminosities of SNe Ia evolve as a function of redshift. The estimated value of $H_0$, however, is only negligibly affected by accounting for possible evolution.

We demonstrate that models of resonant cyclotron radiation transfer in a strong field (i.e. cyclotron scattering) can account for spectral lines seen at two epochs, denoted S1 and S2, in the Ginga data for GRB870303. Using a generalized version of the Monte Carlo code of Wang et al. (1988,1989b), we model line formation by injecting continuum photons into a static plane-parallel slab of electrons threaded by a strong neutron star magnetic field (~ 10^12 G) which may be oriented at an arbitrary angle relative to the slab normal. We examine two source geometries, which we denote "1-0" and "1-1," with the numbers representing the relative electron column densities above and below the continuum photon source plane. We compare azimuthally symmetric models, i.e. models in which the magnetic field is parallel to the slab normal, with models having more general magnetic field orientations. If the bursting source has a simple dipole field, these two model classes represent line formation at the magnetic pole, or elsewhere on the stellar surface. We find that the data of S1 and S2, considered individually, are consistent with both geometries, and with all magnetic field orientations, with the exception that the S1 data clearly favor line formation away from a polar cap in the 1-1 geometry, with the best-fit model placing the line-forming region at the magnetic equator. Within both geometries, fits to the combined (S1+S2) data marginally favor models which feature equatorial line formation, and in which the observer's orientation with respect to the slab changes between the two epochs. We interpret this change as being due to neutron star rotation, and we place limits on the rotation period.

The Ginga data for the gamma-ray burst GRB870303 exhibit low-energy dips in two temporally distinct spectra, denoted S1 and S2. S1, spanning 4 s, exhibits a single line candidate at ~ 20 keV, while S2, spanning 9 s, exhibits apparently harmonically spaced line candidates at ~ 20 and 40 keV. We evaluate the statistical evidence for these lines, using phenomenological continuum and line models which in their details are independent of the distance scale to gamma-ray bursts. We employ the methodologies based on both frequentist and Bayesian statistical inference that we develop in Freeman et al. (1999b). These methodologies utilize the information present in the data to select the simplest model that adequately describes the data from among a wide range of continuum and continuum-plus-line(s) models. This ensures that the chosen model does not include free parameters that the data deem unnecessary and that would act to reduce the frequentist significance and Bayesian odds of the continuum-plus-line(s) model. We calculate the significance of the continuum-plus-line(s) models using the Chi-Square Maximum Likelihood Ratio test. We describe a parametrization of the exponentiated Gaussian absorption line shape that makes the probability surface in parameter space better-behaved, allowing us to estimate analytically the Bayesian odds. The significance of the continuum-plus-line models requested by the S1 and S2 data are 3.6 x 10^-5 and 1.7 x 10^-4 respectively, with the odds favoring them being 114:1 and 7:1. We also apply our methodology to the combined (S1+S2) data. The significance of the continuum-plus-lines model requested by the combined data is 4.2 x 10^-8, with the odds favoring it being 40,300:1.

We have conducted pencil-beam searches for outer solar system objects to a limiting magnitude of R ~ 26. Five new trans-neptunian objects were detected in these searches. Our combined data set provides an estimate of ~90 trans-neptunian objects per square degree brighter than ~ 25.9. This estimate is a factor of 3 above the expected number of objects based on an extrapolation of previous surveys with brighter limits, and appears consistent with the hypothesis of a single power-law luminosity function for the entire trans-neptunian region. Maximum likelihood fits to all self-consistent published surveys with published efficiency functions predicts a cumulative sky density Sigma(<R) obeying log10(Sigma) = 0.76(R-23.4) objects per square degree brighter than a given magnitude R.

A recent study of the rotation of the plane of polarization of light from 160 cosmological sources claims to find significant evidence for cosmological anisotropy. We point out methodological weaknesses of that study, and reanalyze the same data using Bayesian methods that overcome these problems. We find that the data always favor isotropic models for the distribution of observed polarizations over counterparts that have a cosmological anisotropy of the type advocated in the earlier study. Although anisotropic models are not completely ruled out, the data put strong lower limits on the length scale $\lambda$ (in units of the Hubble length) associated with the anisotropy; the lower limits of 95% credible regions for $\lambda$ lie between 0.43 and 0.62 in all anisotropic models we studied, values several times larger than the best-fit value of $\lambda \approx 0.1$ found in the earlier study. The length scale is not constrained from above. The vast majority of sources in the data are at distances closer than 0.4 Hubble lengths (corresponding to a redshift of $\approx$0.8); the results are thus consistent with there being no significant anisotropy on the length scale probed by these data.

We use Bayesian methods to analyze the distribution of gamma ray burst intensities reported in the Third BATSE Catalog (3B catalog) of gamma ray bursts, presuming the distribution of burst sources (``bursters'') is isotropic. We study both phenomenological and cosmological source distribution models, using Bayes's theorem both to infer unknown parameters in the models, and to compare rival models. We analyze the distribution of the time-averaged peak photon number flux, F, measured on both 64 ms and 1024 ms time scales, performing the analysis of data based on each time scale independently. Several of our findings differ from those of previous analyses that modeled burst detection less completely. In particular, we find that the width of the intrinsic luminosity function for bursters is unconstrained, and the luminosity function of the actually observed bursts can be extremely broad, in contrast to the findings of all previous studies. Useful constraints probably require observation of bursts significantly fainter than those visible to BATSE. We also find that the 3B peak flux data do not usefully constrain the redshifts of burst sources; useful constraints require the analysis of data beyond that in the 3B catalog (such as burst time histories), or data from brighter bursts than have been seen by BATSE (such as those observed by the Pioneer Venus Orbiter). In addition, we find that an accurate understanding of the peak flux distributions reported in the 3B almost certainly requires consideration of data on the temporal and spectral properties of bursts beyond that reported in the 3B catalog, and more sophisticated modeling than has so far been attempted.

We use Bayesian methods to study anisotropic models for the distribution of gamma ray burst intensities and directions reported in the Third BATSE Catalog (3B catalog) of gamma ray bursts. We analyze data obtained using both the 64 ms and 1024 ms measuring timescales. We study both purely local models in which burst sources (``bursters'') are presumed to be distributed in extended halos about the Galaxy and M31, and mixed models consisting of a cosmological population of standard candle bursters and a local population distributed throughout a standard Bahcall-Soneira dark matter halo with a 2 kpc core. We find that the purely local models we have studied can account for the 3B data as successfully as cosmological models, provided one considers halos with core sizes significantly larger than those used to model the distribution of dark matter. We infer core sizes for the halo distribution that are smaller than one might expect based on popular semiquantitative arguments, and show why such arguments can lead to unwarranted conclusions. We also find that the 3B data do not constrain the width of power-law luminosity functions for burst sources. Our analysis of mixed models finds two families of models that can successfully account for the data: models with up to 20% of observed bursts in a bright local population visible to ~ 50 kpc; and models with up to 50% of observed bursts in a dim local population visible only nearby (to less than a disk scale height). These models fit as well or better than purely cosmological models. They indicate that a surprisingly large local, anisotropic component could be present whose size is comparable to the sizes of hypothetical classes of bursts inferred from analyses of temporal and spectral characteristics.