Spectral analysis in high-energy astrophysics typically requires sampling from difficult posterior distributions, e.g., multimodal distributions, in a highly structured model with complex data collection mechanisms. Markov chain Monte Carlo (MCMC) methods have been widely explored in such tasks. However, traditional MCMC methods suffer from slow convergence incurred by the local trap problem. In this thesis, we propose a general population-based MCMC strategy, which is able to speed up convergence significantly. Specifically, we initialize multiple chains from dispersed starting values and perform a two-step jump with our proposed between-chain jump. The novel between-chain jumping proposal tries to move to the neighborhood of the iterate in another chain, which encourages full exploration of the parameter space. As for the numerical illustration, we apply the proposed method to fit thermal models to Capella data. The results indicate that our method is effective in two aspects. First, it can be applied in Bayesian model selection to decide which mixture model is more appropriate; secondly, it can be served as an exploratory step to identify modes. The strength shown from the second aspect can help determine the importance of those temperature components and decide upon the relative proportion among those modes. In addition, we show that our method is also useful in other applications including Bayes factor computation for Bayesian model selection, variance component estimation in mixed effect models, and sensor network localization.

Missing data occur frequently in surveys, clinical trials as well as other real data studies. In the analysis of incomplete data, one needs to correctly identify the missing mechanism and then adopt appropriate statistical procedures. Recently, the analysis of missing data has gained more and more attention. People start to investigate the missing data analysis in several different areas. This dissertation concerns two projects. First, we propose a Bayesian solution to data analysis with non-ignorable missingness. The other one is the non-parametric test of missing mechanism for incomplete multivariate data.

First, Bayesian methods are proposed to detect non-ignorable missing and eliminate potential bias in estimators when non-ignorable missing presents. Two hierarchical linear models, pattern mixture model and selection model, are applied to a real data example: the National Assessment of Education Progress (NAEP) education survey data. The results show that the Bayesian methods can correctly recognize the missingness mechanism and provide model-based estimators which can eliminate the possible bias due to non-ignorable missing. We also evaluate the goodness-of-fit of these two proposed models using two methods: the comparison of the real data with the predictive posterior distribution and the residual analysis by cross validation. A simulation study compares the performance of the two proposed Bayesian methods with the traditional design-based methods under different missing mechanisms and show the good properties of the Bayesian methods. Further, we discuss the three commonly used model selection criteria: the Bayes factor, the deviance information criterion and the minimum posterior predictive loss approach. Due to the complicated calculation of the Bayes factor and the uncertainty of the DIC, we conduct the last approach, which fails to correctly detect the real model structure for the hierarchical linear model.

Second, as an alternative to the fully specified model-based Bayesian method, a novel non-parametric test is proposed to detect the missing mechanism for multivariate missing data. The proposed test does not need any distributional assumptions and is proven to be consistent. A simulation study demonstrates that it has well controlled type I error and satisfactory power against a variety of alternative hypotheses.

We were motivated by the two major limitations of the current research approaches on the North Atlantic Oscillation (NAO) based on empirical orthogonal functions (EOF) analysis: (i) long-term stationary assumptions; (ii) lack of measures of uncertainty, and proposed and developed a time-varying low-dimensional representation for spatio-temporal data in this thesis.

The low-dimensional representation is based on a structured spatial covariance matrix using a certain number of structured basis functions with certain parametric forms. Initially, we developed the Parametric Basis Function (PBF) spatial covariance model in a stationary scenario and provided the statistical inference in both maximum likelihood and Bayesian analysis frameworks.

We further extended the model by introducing time-varying parameters to develop the time-varying parametric basis function (TV-PBF) model in the state space model framework. The Bayesian approach with MCMC techniques was used to make inference for the TV-PBF model. The model is able to provide smoothly changing patterns of the 1st EOFs NAO over time which can serve as an alternative representation for the spatio-temporal NAO data.

Systematic instrumental uncertainties in astronomical analyses have been generally ignored in data analysis due to the lack of robust principled methods, though the importance of incorporating instrumental calibration uncertainty is widely recognized by both users and instrument builders. Ignoring calibration uncertainty can cause bias in the estimation of source model parameters and can lead to underestimation of the variance of these estimates. Lee et al. (2011) introduced a so-called pragmatic Bayesian method to address this problem. The method is "pragmatic" in that it introduces an ad hoc technique that simplifies computation by assuming that the current data is not useful in narrowing the uncertainty for the calibration product, i.e., that the prior and posterior distributions for the calibration products are the same.

In the thesis, we focus on incorporating calibration uncertainty into a principled Bayesian X-ray spectral analysis, specifically we account for uncertainty in the so-called effective area curve and the photon redistribution matrix. X-ray spectral analysis models the distribution of the energies of X-ray photons emitted from an astronomical source. The effective area curve of an X-ray detector describes its sensitive as a function of the energy of incoming photons, and the photon redistribution matrix describes the probability distribution of the recorded (discrete) energy of a photon as a function of the true (discretized) energy. Starting with the effective area curve, we follow Lee et al. (2011) and use a principle component analysis (PCA) to efficiently represent the uncertainty. Here, however, we leverage this representation to enable a principled, fully Bayesian method to account for calibration uncertainty in high-energy spectral analysis. For the photon redistribution matrix, we first model each conditional distribution as a normal distribution and then apply PCA to the parameters describing the normal models. This results in an efficient low-dimensional summary of the uncertainty in the redistribution matrix. Our methods for both calibration products are compared with standard analysis techniques and the pragmatic Bayesian method of Lee et al. (2011). The advantage of the fully Bayesian method is that it allows the data to provide information not only for estimation of the source parameters but also for the calibration product; we demonstrate this for the effective area curve. In this way, our fully Bayesian approach can yield more accurate and efficient estimates of the source parameters, and valid estimates of their uncertainty. Moreover, the fully Bayesian approach is the only method that allows us to make a valid inference about the effective area curve itself, quantifying which possible curves are most consistent with the data.

This dissertation focuses on developing statistical models and methods to address data-analytic challenges in astrostatistics---a growing interdisciplinary field fostering collaborations between statisticians and astrophysicists. The astrostatistics projects we tackle can be divided into two main categories: modeling solar activity and Bayesian analysis of stellar evolution. These categories from Part I and Part II of this dissertation, respectively.

The first line of research we pursue involves classification and modeling of evolving solar features. Advances in space-based observatories are increasing both the quality and quantity of solar data, primarily in the form of high-resolution images. To analyze massive streams of solar image data, we develop a science-driven dimension reduction methodology to extract scientifically meaningful features from images. This methodology utilizes mathematical morphology to produce a concise numerical summary of the magnetic flux distribution in solar ``active regions'' that (i) is far easier to work with than the source images, (ii) encapsulates scientifically relevant information in a more informative manner than existing schemes (i.e., manual classification schemes), and (iii) is amenable to sophisticated statistical analyses.

In a related line of research, we perform a Bayesian analysis of the solar cycle using multiple proxy variables, such as sunspot numbers. We take advantage of patterns and correlations among the proxy variables to model solar activity using data from proxies that have become available more recently, while also taking advantage of the long history of observations of sunspot numbers. This model is an extension of the Yu et al. (2012) Bayesian hierarchical model for the solar cycle that used the sunspot numbers alone. Since proxies have different temporal coverage, we devise a multiple imputation scheme to account for missing data. We find that incorporating multiple proxies reveals important features of the solar cycle that are missed when the model is fit using only the sunspot numbers.

In Part II of this dissertation we focus on two related lines of research involving Bayesian analysis of stellar evolution. We first focus on modeling multiple stellar populations in star clusters. It has long been assumed that all star clusters are comprised of single stellar populations---stars that formed at roughly the same time from a common molecular cloud. However, recent studies have produced evidence that some clusters host multiple populations, which has far-reaching scientific implications. We develop a Bayesian hierarchical model for multiple-population star clusters, extending earlier statistical models of stellar evolution (e.g., van Dyk et al. 2009, Stein et al. 2013). We also devise an adaptive Markov chain Monte Carlo algorithm to explore the complex posterior distribution. We use numerical studies to demonstrate that our method can recover parameters of multiple-population clusters, and also show how model misspecification can be diagnosed. Our model and computational tools are incorporated into an open-source software suite known as BASE-9. We also explore statistical properties of the estimators and determine that the influence of the prior distribution does not diminish with larger sample sizes, leading to non-standard asymptotics.

In a final line of research, we present the first-ever attempt to estimate the carbon fraction of white dwarfs. This quantity has important implications for both astrophysics and fundamental nuclear physics, but is currently unknown. We use a numerical study to demonstrate that assuming an incorrect value for the carbon fraction leads to incorrect white-dwarf ages of star clusters. Finally, we present our attempt to estimate the carbon fraction of the white dwarfs in the well-studied star cluster 47 Tucanae.