# Fitting distribution to data with heteroscedastic errors

I want to fit a distribution to a set of univariate observations that were obtained using two different measurement techniques, which have with different measurement errors. One thought would be to fit a mixture model to the data, in which I first fit independent distributions to the values obtained from each of the two measurement techniques, then mix them according to the relative measurement errors. However, this seems like a crude solution since the distributions are constructed independently.

Another thought is using a kernel density estimate with the kernel width dependent on the measurement technique, but for communication and future modeling purposes a parametric distribution is definitely preferable. Is there a relatively straightforward way to formulate this with a Bayesian approach (in PyMC3 or Stan)?

I can provide an analytical framework for the problem, but I'm not familiar enough with PyMC3 or Stan to translate that framework into either of those languages. (I write my own MCMC samplers in Mathematica!)

Let $y = (y_1, \ldots, y_n)$ denote an unobserved variable where $$p(y|\mu,\tau) = \prod_{i=1}^n \textsf{N}(y_i|\mu,\tau)$$ and $\tau$ is a variance. The actual observations are given by $z = (z_1, \ldots, z_n)$, where $$p(z|y,\lambda_1,\lambda_2) = \prod_{i=1}^{n_1} \textsf{N}(z_i|y_i,\lambda_1) \prod_{i=n_1+1}^n \textsf{N}(z_i|y_i,\lambda_2)$$ and $\lambda_j$ is a variance. The model is completed by providing a prior for the parameters $\phi = (\mu,\tau,\lambda_1,\lambda_2)$. This prior will likely impose some independence assumptions such as $$p(\phi) = p(\mu,\tau)\,p(\lambda_1)\,p(\lambda_2) .$$ The joint posterior for the unknowns given the observations $z$ can be expressed as $$p(y,\phi|z) \propto p(z|y,\lambda_1,\lambda_2)\,p(y|\mu,\tau)\,p(\phi) .$$

The goal might be to compute the posterior predictive distribution for the latent variable: $$p(y_{n+1}|z) = \iint \textsf{N}(y_{n+1}|\mu,\tau)\,p(\mu,\tau|z)\,d\mu\,d\tau ,$$ where $$p(\mu,\tau|z) = \iiint p(y,\phi|z)\,dy\,d\lambda_1\,d\lambda_2 .$$ Let $\{(y^{(r)},\phi^{(r)}\}_{r=1}^R$ denote draws from the joint posterior distribution (computed via MCMC). The predictive distribution is approximated by $$p(y_{n+1}|z) \approx \frac{1}{R} \sum_{r=1}^R \textsf{N}(y_{n+1}|\mu^{(r)},\tau^{(r)}) .$$

Whether this predictive distribution is approximately normal or not depends strongly on the prior for $\mu$ and $\tau$. For example, if the prior is parametric, such as $$p(\mu,\tau) = p(\mu)\,p(\tau) = \textsf{N}(\mu|m_0,s_0^2)\,\textsf{Inv-Gamma}(\tau|a_y/2,b_y/2) ,$$ then the predictive distribution will likely be approximately normal. By contrast, if one adopts a nonparametric prior for $(\mu,\tau)$ such as the Dirichlet Process (DP) prior, then the posterior predictive distribution will reflect the underlying distribution for $y$, which could be highly non-normal.

The OP asked for more detail on the DP-related prior. As a preliminary, it is convenient to re-express the prior for $y$ as follows: \begin{align}\label{like_mix} y_i|\theta_y &\stackrel{\text{iid}}{\sim} F(\theta_y) \\ \label{parametric_H} \theta_y &\sim H , \end{align} where $\theta_y = (\mu,\tau)$ and the densities for the distributions $F(\theta_y)$ and $H$ are given by \begin{align} f(y_i|\theta_y) &= \textsf{N}(y_i|\mu,\tau) \\ h(\theta_y) &= \textsf{N}(\mu|m_0,s_0^2)\, \textsf{Inv-Gamma}(\tau|a_y/2,b_y/2) . \end{align} As we acquire data, we learn about $\theta_y$. We now turn to a more general prior that has the property that --- instead of learning about $\theta_y$ --- we learn about the distribution for $\theta_y$.

We generalize the parametric prior by providing a hierarchical prior for $\theta_y$ that involves a Dirichlet Process (DP): \begin{align} \theta_y|G &\sim G \\ G &\sim \textsf{DP}(\xi,H) , \end{align} where $\textsf{DP}(\xi,H)$ represents a Dirichlet Process. In this prior, the prior for $\theta_y$ is the probability measure $G$ which itself is random and depends on $\xi$ and $H$. The mean of $G$ is $H$ (the base distribution): $E[G] = H$ and the variation of $G$ around $H$ is controlled by the $\xi$ (the concentration parameter).

There is an explicit representation for the DP prior that involves an infinite-order mixture. In preparation, let $\psi = (w,\theta)$ denote the parameters of the mixture where $w = (w_1, w_2, \ldots )$ is an infinite collection of nonnegative mixture weights that sum to one, $\theta = (\theta_1, \theta_2, \ldots)$ is a corresponding collection of mixture parameters, and $\theta_c = (\mu_c,\tau_c)$. In addition, let $\delta_{\theta_c}$ denote a point mass located at $\theta_c$. The prior can be expressed as \begin{align} \theta_y|\psi &\sim G = \sum_{c=1}^\infty w_c\,\delta_{\theta_c} \\ w|\xi &\sim \textsf{Stick}(\xi) \\ \theta_c &\stackrel{\text{iid}}{\sim} H . \end{align} In this representation, the randomness of $G$ follows from the randomness of the mixture weights $w$ (which depend on $\xi$) and of the locations $\theta$ (which depend on $H$). As we acquire data, we learn about $\psi = (w,\theta)$ and consequently we learn about the distribution for $\theta_y$.

The form of the prior that will prove most useful is obtained by integrating out $\theta_y$, which has the effect of transfering the infinite mixture from $\theta_y$ to $y$: \begin{align}\label{x_mix} y_i|\psi &\stackrel{\text{iid}}{\sim} \sum_{c=1}^\infty w_c\,F(\theta_c) \\ w|\xi &\sim \textsf{Stick}(\xi) \\ \theta_c &\stackrel{\text{iid}}{\sim} H . \end{align} All that remains is to describe $\textsf{Stick}(\xi)$, the distribution for the mixture weights. (In addition, the DPM can be extended by providing a prior for $\xi$ as well.)

The standard prior for $w$ is given by $$w|\xi \sim \textsf{Stick}(\xi),$$ where $\textsf{Stick}(\xi)$ denotes the stick-breaking distribution given by $$\label{gem1} w_c = v_c \prod_{\ell=1}^{c-1} (1-v_\ell) \qquad\text{where v_c \stackrel{\text{iid}}{\sim} \textsf{Beta}(1, \xi)}.$$ The concentration parameter $\xi$ controls the rate at which the weights decline on average. In particular, the weights decline geometrically in expectation: $$E[w_c|\xi] = \xi^{c-1}\, (1+\xi)^{-c}.$$ Note $E[w_1|\xi] = 1/(1+\xi)$ and $E\!\left[\sum_{c=m+1}^\infty w_c|\xi\right] = \big(\xi/(1+\xi)\big)^m$.

• Great, thank you! Would you be able to explain how a DP prior can be applied here? My (very basic) understanding is that it would be placing a distribution of distributions on $(\mu, \tau)$ so that each data point is modeled as belonging to a different distribution or cluster - is this somewhat similar to a kernel density estimate with a different kernel for each data point? How is this set of distributions specified, e.g., would I have to specify the type of parametric distribution(s) that make up the set of distributions? – Aorus Sep 15 '17 at 15:34
• @Aorus I'll extend my answer to explain how to setup the DP prior. (Your understanding seems roughly correct.) – mef Sep 15 '17 at 16:38