# Maximum Likelihood estimation of a parametric density of a univariate response, given multidimensional data?

I have some data represented by vectors $$\mathbf{x}_1,\mathbf{x}_2,\ldots,\mathbf{x}_n \in \mathbb{R}^m$$, that is, each sample is a collection of $$m$$ features. In addition, I have a response corresponding to each data point: $$y_1,\ldots,y_n \in \mathbb{R}$$.

I would like to train a parametric model on this data, in order to be able to predict the probability density of the response $$y(\mathbf{x})$$ corresponding to a new sample $$\mathbf{x} \in \mathbb{R}^m$$. (I really need to have a probability density for $$y(\mathbf{x})$$, not just a punctual estimation.)

So my idea would be to select a family of parametric densities defined through a parameter $$\boldsymbol{\theta}\in\mathbb{R}^k$$, $$f(y;\boldsymbol{\theta})$$. Then, I would like to run a maximum-likelihood procedure to find the best possible function $$\theta$$ from a class of functions $$\Theta$$ mapping $$\mathbb{R}^m\mapsto\mathbb{R}^k$$: $$\max_{\theta \in \Theta} \sum_{i=1}^n \log f(y_i,\theta(\mathbf{x}_i)).$$

Question 1 Is this a well-studied problem, and if yes, how is it called in the literature? Or is there an alternative approach that is known to work better?

Question 2 My first, naive idea, is to restrict myself to a class $$\Theta$$ of linear functions, so the above problem reduces to determining $$km$$ coefficients. In this case, any reference to some algorithm performing this task would be appreciated.

Question 3 If we consider a family of gaussians whose parameters $$\mu$$ and $$\sigma$$ are linear functions of $$\mathbf{x}$$, i.e. $$y \sim \mathcal{N}(\boldsymbol{\mu}^T\mathbf{x}, (\boldsymbol{\sigma}^T\mathbf{x})^2)$$ for some vectors $$\boldsymbol{\mu},\boldsymbol{\sigma}\in \mathbb{R}^m$$, then the optimization problem to solve is $$\min_{\boldsymbol{\mu},\boldsymbol{\sigma}\in \mathbb{R}^m}\quad \sum_{i=1}^n \frac{1}{2}\left(\frac{y_i-\boldsymbol{\mu}^T\mathbf{x}_i}{ \boldsymbol{\sigma}^T\mathbf{x}_i} \right)^2 + \log(\boldsymbol{\sigma}^T\mathbf{x}_i).$$ What is the state-of-the-art to solve this problem ?

Or maybe, because we want small predicted variances, we could add a penalization term: $$+ \lambda \sum_i (\boldsymbol{\sigma}^T\mathbf{x}_i )^2$$ in the function to minimize?

Question 4 If we assume that $$f$$ is a mixture of laws $$f_1,\ldots,f_N$$ defined by a vector of probabilities $$\boldsymbol{\pi}$$, is there a way to construct a function $$\pi$$ that maps $$\mathbf{x}$$ to the probability vector $$\boldsymbol{\pi}=\pi(\mathbf{x})$$, so that the density of the response $$y(\mathbf{x})$$ is given by $$\sum_{j=1}^N \pi_j(\mathbf{x}) f_j(y).$$ I have the feeling that using a neural network might be appropriate to do construct the function $$\pi$$, but again, I would be grateful if someone could point out a good reference for this specific problem.

Question 5

Here is a tentative approach to answer Question 4. How appropriate do you think this is? Is this a standard approach (or a stupid one)? Do you know some references?

We could use the EM algorithm to model the augmented vectors $$(\mathbf{x}_i,y_i)\in\mathbb{R}^{m+1}$$ as coming from a mixture of multivariate gaussian variables: $$\left(\left[\begin{array}{c}\mathbf{x}_i\\y_i\end{array}\right] | Z=k\right) \sim \mathcal{N}\left( \left[\begin{array}{c}\boldsymbol{\mu}_k\\ m_k\end{array}\right], \left[\begin{array}{cc}\Sigma_k & \boldsymbol{v}_k\\ \boldsymbol{v}_k^T & \sigma_k^2\end{array}\right] \right), \quad\textrm{with}\ P(Z=k) = w_k.$$

Then, the distribution of the response $$y$$ given the value of a new sample $$\mathbf{x}$$ and the value of the latent variable $$Z$$ would be:

$$(y|\mathbf{x},Z=k) \sim \mathcal{N}\left( m_k + \boldsymbol{v}_k^T \Sigma_k^{-1} (\mathbf{x}-\boldsymbol{\mu}_k),\ \sigma_k^2-\boldsymbol{v}_k^T \Sigma_k^{-1}\boldsymbol{v}_k \right).$$

Finally, we obtain the distribution of $$y|\mathbf{x}$$ as a mixture of gaussians, with weights $$P(Z=k|\mathbf{x})$$ that can be obtained from the Bayes formula.

Basically, the idea is to assume that the responses $y_1,\ldots,y_n$ come from a multivariate normal distribution $\mathcal{N}(\boldsymbol{\mu}, \Sigma)$, whose mean and variance depend from the data $X=(\mathbf{x}_1,\ldots,\mathbf{x}_n)$. A typical approach is to assume a particular form for the functions $\mu()$ and $\Sigma()$, for example the $\mu_i$'s can be equal to some affine function applied to $\mathbf{x}_i$, and $\Sigma_{ij}$ can be a combination of some squared exponential kernel $K(\mathbf{x}_i,\mathbf{x}_j)$ and a diagonal noise $\sigma_{\rm{noise}}^2 \delta_{ij}$. These functions depend from a set of hyperparameters $\boldsymbol{\theta}$, which can be optimized by maximum-likelihood. Then, if we try to predict the answer $y^*$ of a new sample $\mathbf{x}^*$, we can use the formula that gives the posterior $y^* | (y_1,\ldots,y_n)$: $$y^* | (y_1,\ldots,y_n) \sim \mathcal{N}\left(\mu(\mathbf{x}^*) + \mathbf{v_*}^T \Sigma(X)^{-1} (\boldsymbol{y}-\boldsymbol{\mu}(X)),\quad \sigma_*^2 - \mathbf{v_*}^T \Sigma(X)^{-1} \mathbf{v_*} \right),$$ where $\sigma_*^2:=\Sigma(\mathbf{x}^*,\mathbf{x}^*)$ is the prior variance at $\mathbf{x}^*$ and $\mathbf{v_*}$ is the vector of prior cross-covariances between $\mathbf{x}^*$ and the data: $(\mathbf{v_*})_i = \Sigma(\mathbf{x}^*,\mathbf{x}_i)$.
Question 3 For my particular application, I can indeed solve this maximum likelihood problem to find the vectors $\boldsymbol{\mu}$ and $\boldsymbol{\sigma}$. Then, I can use these optimal vectors of coefficients to construct a starting point for the (much harder) problem of optimizing the hyper-parameters of a gaussian process, cf. my answer to Question 1. I found that the following strategy works well: solve alternately the optimization problem with respect to $\boldsymbol{\mu}$ and $\boldsymbol{\sigma}$. The optimization in $\boldsymbol{\mu}$ is a standard least square problem for a fixed $\boldsymbol{\sigma}$, and the optimization in $\boldsymbol{\sigma}$ can be done by any gradient-based algorithm, e.g. the BFGS algo.
• What I proposed in Question 5 seems to be close to the solution provided by Gaussian processes (cf. my answer to question 1), but uses a mixture of Gaussian rather than a Kernel, and works in the space of extended vectors $(\mathbf{x}_i,y_i)$. I wonder if the link between these two approaches has already been studied ? Oct 5 '15 at 11:02