dependency between regression coefficients and probability distribution

Question

let's consider regression problem. Given set of training data $\{(x_i,y_i)\}_{i=1}^N$, $x_i \in \mathbb{R}^n$ and $y_i \in \mathbb{R}$, find prediction function $y = f(x)$, e.g. in RBF regression case $y = \sum_{i=0}^m a_i g(x_i,c_i)$, where $g(x_i,c_i)$ is radial function with center $c_i$. It's well known that solution can be find by least-squares method from training data as $a = (G^\mathrm{T} G)^{-1} G^\mathrm{T} y$.

Here is another question, if we don't have training sample, but we know the joint distribution density $p(x,y)$, e.g. $p(x,y) = \sum_{i=1}^r b_i k_i(x,y)$. How to estimate regression coefficients $a_i$ from distribution $p(x,y)$ (coefficients $b_i$)?

Answer 1

You can consider the probability distribution as a continuum of sample points and use the same formula, $a=(G^\mathrm{T}G)^{-1}G^\mathrm{T}y$, except that matrix multiplication now means integrating over all $x$ and $y$, weighted by the probability distribution, instead of summing over a finite number of points.

Edit in response to TheBug's comment:

Sorry, I was too vague and general in my original answer.

I don't know what "radial function with center $c_i$" is in your problem statement; I'll write it down as I understand it and I think you'll be able to adapt it to your notation/case.

We want to find $a_i$ to make $y=\sum_i a_i g_i (x)$ a best fit.

Discrete case with sample points:

$\sum_i \left(y_i - \sum_j a_j g_j(x_i)\right)^2 \rightarrow \mathrm{min}$

$\sum_i \left(y_i - \sum_j a_j g_j(x_i)\right)g_k(x_i) = 0$ for all $k$

In matrix form: $G^\mathrm{T}y - G^\mathrm{T}Ga=0$ with $G_{ik} = g_k(x_i)$ and $(G^\mathrm{T}G)_{kj} = \sum_i G_{ik}G_{ij}$, and thus $a=(G^\mathrm{T}G)^{-1}G^\mathrm{T}y$.

Continuous case with probability distribution:

\[\int \mathrm{d}^nx \;\mathrm{d}y\; p (x,y) \left(y - \sum_j a_j g_j(x))\right)^2 \rightarrow \mathrm{min}\]

\[\int \mathrm{d}^nx \;\mathrm{d}y\; p (x,y) \left(y - \sum_j a_j g_j(x))\right)g_k(x)= 0\]

Now the equivalent of $G$ is an object with one discrete index (j/k) and one continuous index (x and y). The equivalent of the matrix multiplication $G^\mathrm{T}G$ is the integration

\[\int \mathrm{d}^nx \;\mathrm{d}y\; p (x,y) g_j(x) g_k(x) =: A_{kj}\;,\]

which integrates over the continuous index (corresponding to the sum over $i$ in the discrete case) and leaves a "normal" matrix $A$ with the two discrete indices $j$ and $k$ which plays the role of $G^\mathrm{T}G$. The equivalent of the matrix multiplication $G^\mathrm{T}y$ is the integration

\[\int \mathrm{d}^nx \;\mathrm{d}y\; p (x,y) y g_k(x) =: b_k\;,\]

which integrates over the continuous index and leaves a "normal" vector $b$ with the discrete index $k$. The result is again a "normal" matrix equation with discrete indices: $Aa=b$, and thus $a=A^{-1}b$.

I hope that was more helpful; sorry for being a bit terse the first time around; feel free to ask more questions if I haven't made it clear enough.

Edit in response to the additional question about convergence with increasing training set size:

I think phrasing the question as we both did (i.e. the probability in every neighbourhood of the "correct" values tending to 1), we can prove it using the continuity of the inverse.

If we were trying to calculate the expected values of the fitted parameters, we'd have to deal with the non-linearity and potential singularity of the matrix inversion. The average of $G^\mathrm{T}G$ over the probability distribution is already the matrix $A$, independent of the size of the training set, but the same is not true for the inverses of these matrices. In fact the inverse of $G^\mathrm{T}G$ doesn't even exist until the training set contains at least as many points as we're trying to fit, and even then, there's always a subspace (of measure 0) of training sets (of measure 0) that lead to a singular matrix, and in integrating over the inverse we'd be integrating over all these singularties.

But if we're only interested in the probability in neighbourhoods tending to 1, I think we don't have to worry about all that. The expected values of $G^\mathrm{T}G$ and of $G^\mathrm{T}y$ are the "correct" values, and the variances of their entries decrease as $n^{-1}$ with the size $n$ of the training set. Thus, the probability of $G^\mathrm{T}G$ and $G^\mathrm{T}y$ being within a given neighbourhood of their "correct" values tends to 1. If we choose the neighbourhood small enough, it will not contain any singular matrices. By the continuity of matrix inversion (restricted to such a well-behaved neighbourhood), for every given neighbourhood of the "correct" parameters we can choose a sufficiently small neighbourhood of the "correct" values of $G^\mathrm{T}G$ and $G^\mathrm{T}y$ such that $a=(G^\mathrm{T}G)^{-1}G^\mathrm{T}y$ will lie in the given neighbourhood, and hence the probability of the fitted parameter values lying within the given neighbourhood will also tend to 1.