Probability distribution of two lines intersection in 3D I am given two lines in 3D spaces. Each line is attached to its own coordinate systems which have distinct origins.
Each line always goes through its coordinate system origin but its direction is uncertain. This uncertainty is modeled by a gaussian multivariate distribution with respect to the azimut and elevation (polar coordinates) defining the line direction.
Thus, the uncertainty for each line can be represented by a conic whose axe is defined by the azimut and elevation expectation. This is what I tried to draw on following picture (my apologizes, it is not perfect).
line intersection
As gaussian distributions have an infinite support, the conic shape and delimitation is somewhat arbitrary, the true shape would be e.g. to compute
the volume $x^T\Sigma x \le 1$, $\Sigma$ being the covariance matrix of the PDF of the given line.
My main question is that I would like to characterize the probability of lines to intersect. That is to say I would like at least to compute the expectation of the intersection and its matrix covariance.
One of the difficulty comes from the fact that line PDFs are expressed in distincts coordinate systems (distincts origins at least).
Let's use another cartesian coordinate system, i.e. whose origin $O$ is distinct from $O_1$ and $O_2$, the origin of respective line coordinate systems.
The probability $p({\mathbf x})$ that lines intersect around position ${\mathbf x}$ is the product of both lines multivariate PDFs :
$$p({\mathbf x}) = p_1(Q_1({\mathbf x}))p_2(Q_2({\mathbf x}))$$
where $Q_1$ and $Q_2$ transform the cartesian coordinate of ${\mathbf x}$ into its coordinates in each line system, except that we only need the azimut and elevation (range ignored).
So, I would like to compute :
$$E(\mathbf{x}) = \int \mathbf{x}p(\mathbf{x})\,d\mathbf{x}$$
which turns out to be rather complex :
$$E({\mathbf x}) = \frac{1}{(2\pi)^2\sqrt{|\Sigma_1||\Sigma_2|}}\int\mathbf{x}e^{-\frac{1}{2}\biggl( (Q_1(\mathbf{x})-\mu_1)^T\Sigma_1^{-1}(Q_1(\mathbf{x})-\mu_1) + (Q_2(\mathbf{x})-\mu_2)^T\Sigma_2^{-1}(Q_2(\mathbf{x})-\mu_2)\biggr)}\,d\mathbf{x}$$
Note that $\mathbf{x}$ is a 3D vector, but $Q_i(\mathbf{x})$ is a 2D vector composed of the azimu and elevation towards the given point in the $i^{th}$ line coordinate system. So I use the bivariate normal distribution.
I would like to avoid computing numerically this integral. With univariate normal distribution the maximum probability for $X$ is also its expectation $E(X)~,~~X\sim N(\mu, \sigma^2)$.
Do you think computing the maximum of $p(\mathbf{x})$ by deriving it would give me the expectation ? Then, how to compute the covariance matrix ?
Thanks for your kind attention
Regards
 A: Let's suppose that the coordinate systems are all identical, so the $Q_1$ and $Q_2$ are the same, and that we're working with a point where the elevation is zero, and the azimuth is far from $0$ or $2\pi$, i.e., far from the singularity in the parameterization. 
The thing you've written down

The probability $p({\mathbf x})$ that lines intersect around position ${\mathbf x}$ is the product of both lines multivariate PDFs :$$p({\mathbf x}) = p_1(Q_1({\mathbf x}))p_2(Q_2({\mathbf x}))$$

is not correct. For if $U$ is a small box around $y$ (i.e., a small region of space corresponding to a small range of azimuth/elevation pairs), then $$\int_U p_1(Q_1(x))~dx
$$
represents the probability that the first line has azimuth and elevation nearly equal to those of $y$...which includes all points in the cone from the origin through all points of $U$, many of which (like, say $37y$) are likely to be far from $y$. A similar statement is true for $p_2(Q_2(...))$. Note that even if each of these probabilities is high, the lines in question generally will not intersect! So in taking a limit as the size of $U$ goes to zero, you're still not computing anything about the probability of two lines intersecting. 
It's possible that despite this, you still actually want to do the integral above. If so, go for it. 
I'd also like to point out that in problems like this, working in azimuth/elevation coordinates is usually a weird choice (unless you're considering, say, a pair of telescopes with azimiuth/elevation drives, and you're uncertain about the control circuitry or something). One problem is that the values you might compute will change completely when the orientation of the coordinate systems (but not the actual point $x$) is changed --- that's usually a bad sign: you'd really like your answers to be a function of the geometry rather than of the coordinate system you've chosen for parameterizing your geometry. Without knowing your particular problem, I can't be sure this is relevant, but I thought I'd mention it. 
A: Thanks for your answeer.
In response, I would like to give you a bit more of the typical context of this question : imagine that you manage several sensors, each one tracking several objects. Sensors are usually placed at distinct locations and many of them will only give you the azimut and elevation of the detected objects in the coordinate system originating on the sensor.
Then, you are asked whether spots detected by sensors correspond to the very same object or not.
This is an issue known as data fusion, there are already some good algorithms to address it. I focus on a particular issue here, I cannot tell you more.
Spots from each sensor comes with some uncertainty due to the sensor measure error itself, the noise received by the sensor and data processing within the sensor itself.
So, if the criteria for deciding to associate some spots from each sensor was that their corresponding lines cross, I would never make any association.
This is why I try to assess the PDF of line intersection in 3D space, that is to say the probability that lines will cross the same space volume, that is to say lines are close enough to consider they point towards the same object.
Trouble is that the density function is not a PDF, because not normalized. However one can still compute an expectation and a covariance matrix.
Regards
A: I am trying in the sequel to develop as much rigorously as possible how I would like to solve my first question.
3D space is represented by a cartesian coordinate system  $(O, x, y, z)$  originating in point $O$.
We define two lines in 3D space, $L_1$ and $L_2$, characterized each by one point $O_i$,
together with their azimut ($a_i$) and elevation ($e_i$) in the coordinate system
attached to point $O_i \scriptstyle{(\ne O)}$ with same axes as $(O, x, y, z)$.
The lines direction is subject to some uncertainties which are represented by a bivariate
gaussian distribution with respect to azimut and elevation angles :
\begin{equation}
p_i(\mathbf{x_i}) = \frac{1}{2\pi\sqrt{|\Sigma_i|}}e^{-\frac{1}{2}[(\mathbf{x_i} - \mu_i)^T\Sigma_i^{-1}(\mathbf{x_i} - \mu_i)]}
\end{equation}
with
\begin{equation}
  \mathbf{x_i} = \begin{bmatrix}
    a_i \\ e_i
  \end{bmatrix}~,~~
  \mathbf{\mu_i} = \begin{bmatrix}
    \mathrm{E}(A_i) \\  \mathrm{E}(E_i)
  \end{bmatrix}~,~~
  \Sigma_i = \begin{bmatrix}
    (\sigma_i^a)^2 & \rho_i \\
    \rho_i & (\sigma_i^e)^2 \\
  \end{bmatrix}~,~~\rho_i = \mathrm{Cov}(A_i, E_i)
\end{equation}
The probability of lines to meet in one point is equal to zero as that of any punctual
event in a continous space. However, one can speak about the probability of each line to
cross a given volume around one point.
Say $M$ is a point with $(x, y, z)^T$ coordinates in the cartesian system. Its
coordinates in each line coordinate system are :
\begin{equation}
  \begin{bmatrix}
    r \\ a_i \\ e_i
  \end{bmatrix}
  =\begin{bmatrix}
  \dots   \\ \arctan(\frac{y -y_i}{x -x_i}) \\ \arctan\biggl(\frac{z -z_i}{\sqrt{(x -x_i)^2 + (y - y_i)^2}}\biggr)
\end{bmatrix}=
\begin{bmatrix}
    \dots \\ \mathbf{x_i}
  \end{bmatrix}
\end{equation}
where $(x_i, y_i, z_i)^T$ are the coordinates of point $O_i$ in the same cartesian
system as $M$.
The probability for each line to cross a small volume $V$ around that point is then :
\begin{equation}
P[L_i \cap V(M)\ne \emptyset ] = p_i(\mathbf{x_i})\,da_i\,de_i
\end{equation}
Assuming PDF of both lines are independent, the probability of both lines to cross the
same volume around point $M$ is :
\begin{equation}
p_1(\mathbf{x_1})p_2(\mathbf{x_2})\,da_1\,de_1\,da_2\,de_2
\end{equation}
I would like then to integrate this formula, but it appears more complex as I first thought, because in latest formula, $a_1$, $a_2$, $e_1$ and $e_2$ are linked by point M.
How could I change the integration variables into these of point M (x, y, z) ?
