Nonparametric estimation of an unknown probability distribution using maximum likelihood and Bayesian approaches

2013-08-01Suppose we observe a random sample Y₁,...,Y_N, of independent but not necessarily identically distributed random variables Yᵢ,∈ ℝᵈ, for i = 1,...,N. Assume also that the conditional density of Yᵢ given θᵢ is known and denoted by pᵢ(Yᵢ|θᵢ), where the θᵢ’s are unobserved random parameters that are independent and identically distributed with common but unknown distribution function F. ❧ The objective is to estimate F given the data Y₁,...,Y_N. We used two different approaches to get the estimate of F: Nonparametric Maximum Likelihood (NPML) and Nonparametric Bayesian (NPB). ❧ For the nonparametric maximum likelihood approach, convex analysis shows that the maximum likelihood (ML) estimator Fᴹᴸ of F is a discrete distribution with no more than N support points where N is the number of subjects in the population. The position and weights of the support points are unknown. We developed four different methods of estimation: NPEM 1, NPEM 2, EM1+IP and EM2+IP. These methods are based on the Expectation–Maximization (EM) and Primal-Dual Interior Point (IP) algorithms. The NPEM 1 algorithm uses the EM algorithm to maximize the log-likelihood function over the set of support points and weights. NPEM 2 fixes the large grid of point at first and them uses EM algorithm to maximize the log-likelihood only over the set of weights. The NPEM 2 algorithm is much faster than the NPEM 1 algorithm, but to find the optimal solution a very large grid of points has to be used. Therefore in case of multidimensional parameters, NPEM 2 can be impractical. The EM1+IP algorithm is a hybrid of the NPEM 1 and IP algorithms in the way that support points are updated the same way as in NPEM 1, but the weights are updated using Primal-Dual Interior Point (IP) algorithms. This method is much faster then NPEM 1 and comparable in speed with NPEM 2, however it does not need a large grid of points. EM1+IP requires no more then N support points to start the estimation process and therefore can be used in cases of multidimensional parameters. The EM2+IP algorithm is a hybrid of NPEM 2 and IP. It starts with a large grid of support points and then updates the weights using the IP algorithm. This method gives results instantly if compared to the speeds of the rest of the methods. However the optimality of the results, as in case of the NPEM 2 algorithm, very much depends on the initial grid of points. ❧ Another highlight of this thesis is an ability of building nonparametric maximum likelihood estimates of distributions of the parameters of linear stochastic models, i.e. population models with process noise. The advantages of the models with process noise are that, in addition to the measurements errors, the uncertainties in the model itself are taken into the account, i.e. process noise. For example, in case of the PK problem the errors like dose errors, sample timing errors, model misspecifications. are not included in the deterministic models. The stochastic models on the other hand can overcome this problem. ❧ We used linear Kalman-Bucy filtering [J70] to find the likelihoods of these models and then employed all the above methods to find the maximum likelihood estimates of an unknown distribution F. ❧ In the NPB approach, the distribution F is itself a random variable. The first requirement is to place a prior distribution on F. A convenient way of doing this is with the so-called Dirichlet process (DP) prior, see [F73, WW97, WW98]. A Dirichlet process prior on a random distribution F is defined by a probability distribution G₀ and a scale parameter α₀, where G₀ is our prior guess of F and the parameter α₀ will be the strength of the guess. We write DP(α₀G₀) for the Dirichlet process. ❧ The main reason we consider the NPB approach in the PK Population Problem is because of its ability to calculate Bayesian credibility intervals, no matter what the sample size is. This is not possible with NPML in general, however it is possible to get bootstrapped confidence intervals for NPML which are somewhat similar to Bayesian credibility intervals. In both cases we get approximations of accuracy of the estimates. Even though NPML and NPB give us an estimates of accuracy, the Bayesian nature of the estimated population distribution F is more preferable. ❧ Two algorithms were developed to find NBP estimates of the distribution F, Marginal NPB and Full Conditional NPB. The Marginal NPB approach uses Gibbs sampling to estimate the posterior expectation of the distribution F given the data, but it does not solve the classical NPB problem of finding the posterior distribution of F given the data. The Full NPB approach uses the Stick-Breaking representation of the Dirichlet process and gives a way to sample from the posterior distribution of F given the data. This approach allows us to build credibility intervals on all the desired quantities of the distribution F. ❧ Finally we established a very nice connection between NPML and NPB methods. Using the Full NPB approach we can sample potential estimates of the true distribution F, moreover we can calculate the value of the log-likelihood function for each of those distributions. All these values of likelihoods are smaller than or equal to the maximum likelihood value that was found by NPML methods. The equality is possible with a big enough number of draws from the posterior distribution of F given the data