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. 
 A: 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.
