Perhaps you can use a model of a piecewise Gaussian variable. Suppose the variable has a mean [itex]a[/itex] and different standard deviation for [itex]x > a[/itex] and [itex]x < a[/itex], i.e. its distribution is:
[tex]
\varphi(x) = \left\{\begin{array}{ll}<br />
A_{1} \exp\left(-\frac{(x - a)^{2}}{2 \sigma^{2}_{1}}\right)&, x > a \\<br />
<br />
A_{2} \exp\left(-\frac{(x - a)^{2}}{2 \sigma^{2}_{2}}\right)&, x < a<br />
\end{array}\right.[/tex]
You have to adjust [itex]A_{1}[/itex] and [itex]A_{2}[/itex] so that:
[tex]
E(X) - a = \int_{-\infty}^{\infty}{(x - a) \varphi(x) \, dx} = 0 \Rightarrow A_{1} \int_{0}^{\infty}{t e^{-\frac{t^{2}}{2 \sigma^{2}_{1}}} \, dt} = A_{2} \int_{0}^{\infty}{t e^{-\frac{t^{2}}{2 \sigma^{2}_{2}}} \, dt} \Rightarrow A_{1} \, \sigma^{2}_{1} = A_{2} \, \sigma^{2}_{2}[/tex]
Of course, the probability density must be normalized:
[tex]
\int_{-\infty}^{\infty}{\varphi(x) \, dx} = 1 \Rightarrow A_{1} \, \int^{\infty}_{0}{e^{-\frac{t^{2}}{2\sigma^{2}_{1}} \, dt} + A_{2} \, \int^{\infty}_{0}{e^{-\frac{t^{2}}{2\sigma^{2}_{2}} \, dt} = 1 \Rightarrow \sqrt{\frac{\pi}{2}} \left(A_{1} \, \sigma_{1} + A_{2} \, \sigma_{2} \right) = 1[/tex]
These two equations allow you to express [itex]A_{1/2}[/itex] in terms of [itex]\sigma_{1/2}[/itex]. Try to find the variance of the variable.
Next, consider the variable:
[tex]
\varepsilon_{i} = a \, X_{i} + b \, Y_{i} + c, \; a^{2} + b^{2} = 1, \; i = 1, \ldots, N[/tex]
If [itex]X_{i}[/itex] and [itex]Y_{i}[/itex] have the above distribution, what is the expectation value and variance for [itex]\varepsilon_{i}[/itex]?
Approximate these variables as having an approximately Normal distribution with the above expectaion values and variances and use the maximum likelihood method, which would reduce to a least-squares method to estimate the parameters of the general linear dependence:
[tex]
a \, x + b \, y + c = 0, \; a^{2} + b^{2} = 1[/tex]