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

Question

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.

Answer 1

So, after some weeks I have come up with some answers to some of my above questions.

Question 1 It turns out that this problem has been studied a lot within the framework of Gaussian processes. A good reference seems to be "Machine Learning: A Probabilistic Perspective", from Kevin P. Murphy.

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.