Likelihood function for a hyperplane I have a hyperplane defined by $y = b_1x_1 + b_2x_2 + b_3x_3$. 
I have a list of $n$ data points $(y,x_1,x_2,x_3)_i, i=1\ldots n$. 
These data points have been measured with errors drawn from a normal distribution known variances $(\sigma^2_y, \sigma^2_{x_1}, \sigma^2_{x_2}, \sigma^2_{x_3})_i$. Each set of data points has a different set of known variances. 
I want to calculate the likelihood of a set of coefficients $(b_1, b_2, b_3)$, i.e., the probability of the data given those coefficients. What is the appropriate likelihood function? 
 A: The tricky part is that there are many implicit parameters.  This becomes clear when the problem is carefully formulated, so let's distinguish the reality from the data. 
Conventionally, we use Greek letters for the (unknown) reality, so express the hyperplane as the set of points
$$\mathcal H(\beta) = \{(x_1,x_2,x_3,y)\mid y - \beta_1 x_1 - \beta_2 x_2 - \beta_3 x_3= 0\}.$$
Your model supposes the probability density of an observation $(x_1,x_2,x_3,y)$ being made at some point with unknown coordinates $(\xi_1,\xi_2,\xi_3,\eta)$ is proportional to
$$f(x_1,\ldots,y; \xi_1,\ldots,\eta; \sigma_{x_1},\ldots,\sigma_y) \propto \frac{1}{\sigma_{x_1}\cdots\sigma_y}\exp\left(-\frac{(x_1-\xi_1)^2}{2\sigma_{x_1}^2} - \cdots --\frac{(y-\eta)^2}{2\sigma_{y}^2}\right).$$
Substituting $\beta_1\xi_1+\beta_2\xi_2+\beta_3\xi_3$ for $\eta$ expresses this as a function of the $x_i,$ $y,$ the $\xi_i,$ the $\sigma_{*},$ and the $\beta_i.$  The $x_i$ and $y$ are the data; the $\sigma_{*}$ are known values (it might therefore be notationally better to use Latin letters and call them $s_{x_1}$ etc); and the rest--the $\xi_i$ and $\beta_i$--are the parameters. Assuming your $n$ measurements are independent, the likelihood will be the product of these probability density values.  Since the data and the $\sigma_{*}$ are known, the likelihood will be a function of $3n$ $\xi_i$ values and the three $\beta$ values: $3(n+1)$ parameters in toto.

For any given $\beta=(\beta_1,\beta_2,\beta_3),$ however, the likelihood will be maximized by projecting each $(x_1,x_2,x_3,y)$ orthogonally onto $\mathcal{H}(\beta)$ with respect to the Euclidean metric determined by the corresponding $\sigma_{*},$ thereby estimating all the $(\xi_1,\xi_2,\xi_3)$ values conditional on $\beta.$  Consequently, the problem of maximizing the likelihood reduces to a three-dimensional optimization.
The orthogonal projection of $(x_1,\ldots,y)$ onto $\mathcal{H}(\beta)$ with respect to the $(\sigma_{x_1}, \ldots, \sigma_y)$ metric is
$$(\xi_1,\ldots,\eta) = (x_1,\ldots,y) - \frac{\beta_1x_1 + \cdots + \beta_3x_3 - y}{\sigma_{x_1}\beta_1^2+\cdots+\sigma_{x_3}\beta_3^2 + \sigma_y}(\sigma_{x_1}\beta_1,\ldots,\sigma_{x_3}\beta_3,-\sigma_y).$$
Upon substituting the right hand side for the left hand variables in the likelihood you will obtain an explicit formula for it in terms of $\beta$ (given the data and the $\sigma_{*}$).
