Conditional joint probability of Gaussian random variables Update: I figured out a way to solve this (still some questions remain)
The problem basically in all its generalization is:  


*

*We are given a matrix $W \in \mathbb{R}^{m \times k}$ whose entries $w_{ij}$ are independent Gaussian random variables distributed as per $w_{ij} \sim \mathcal{N}(0,\sigma^2)$.

*We are given that the sum of each row $R_i > 0$. i.e., $\sum_{j=1}^{k}w_{ij} = R_i  > 0$ for $1 \le i \le m$.

*We have to find the probability that the sum of each column $\sum_{i=1}^{m}w_{ij} = C_j > 0$.


My method: (after referring Density of Gaussian Random variable conditioned on sum)  
We need to calculate (conditioning on row sums is implicit) 
$$P(C_1 > 0)P(C_2 > 0 | C_1> 0)...P(C_k > 0 | C_1, C_2...C_{k-1} > 0)$$
For $P(C_1 > 0)$ we use the linked answer to get $\mu_{C_1} = m\sigma\sqrt{\frac{2}{\pi k}}$ and $\sigma_{C_1}^2 = m\sigma^2\left(1 - \frac{2}{\pi k}\right)$.
Here I have approximated the skew normal distribution with a normal distribution.
Now to calculate $P(C_r > 0 | C_1, C_2...C_{r-1} > 0)$ we can first compute the posterior distribution of $w_{ij}$ conditioned on $C_{j} > 0$ as
$$w_{ij}|_{C_j > 0} \sim \hat{\mathcal{N}}\left(\sigma\sqrt{\frac{2}{\pi m}}, \sigma^2\left(1 - \frac{2}{\pi m}\right)\right)$$
Here I have put a hat showing that I am approximating it by Gaussian. Now the posterior means and variances can be obtained using a truncated normal distribution for $R_i$.
However when I try this out and simulate it, the above method is underestimating the actual probability. Whereas if I simply consider the $C_j$'s independent and multiply their individual probabilities then I am little over-estimating but it is coming much closer to the actual answer.
Is it possible that the approximation of skewed distribution to normal distribution might be interfering? If so is there a fix to this problem?
Thanks
I am leaving the old version of the question as it is:
I have 6 random variables $X, Y, Z, U, V, W$ all iid distributed as $\mathcal{N}(0, \sigma^2)$.
Also I am given that $X + Y + Z > 0$ and  $U + V + W > 0$. I have to find the probability that $X + U > 0$ and $Y + V > 0$ and $Z + W > 0$. 
Now if these three random variables $X + U = A, Y + V = B, Z + W = C$ were independent then I could have simply calculated the posterior probabilities of each of these variables as described in the linked answer and then multiplied them together.
However I think they are dependent so I have to do something like
$$P(A > 0, B > 0, C > 0) = P(A > 0)P(B>0 | A > 0)P(C>0|A>0,B>0)$$
Now conditional probability of $A > 0$ can be got from the linked answer. Then to calculate $P(B>0 | A > 0)$ I shall have to redo the analysis after considering the posterior probabilities of $X$ and $U$.
Am I correct when I say that the random variables are dependent? If so how should I approach this problem?
I shall be grateful if someone can guide me or give me some hints or pointers regarding this.
Thanks
 A: (Incomplete) First you should know that 
$$
\begin{bmatrix}
A \\ B \\ C \\ D 
\end{bmatrix} = 
\begin{bmatrix}
X + U \\ Y + V \\ Z + W \\ X + Y + Z 
\end{bmatrix}
= \begin{bmatrix}
1 & 0 & 0 & 1 & 0 & 0 \\
0 & 1 & 0 & 0 & 1 & 0 \\
0 & 0 & 1 & 0 & 0 & 1 \\
1 & 1 & 1 & 0 & 0 & 0 \\
\end{bmatrix}
\begin{bmatrix}
X \\ Y \\ Z \\ U \\ V \\ W
\end{bmatrix}$$
So $(A, B, C, D)$ is just an affine transform of $(X, Y, Z, U, V, W)$ and thus it also jointly follows a multivariate normal. Obviously the mean vector is still $\mathbf{0}$ and the covariance matrix is given by
$$ \begin{bmatrix}
2\sigma^2 & 0 & 0 & \sigma^2  \\
0 & 2\sigma^2 & 0 & \sigma^2  \\
0 & 0 & 2\sigma^2 & \sigma^2  \\
\sigma^2 & \sigma^2 & \sigma^2 & 3\sigma^2 \\
\end{bmatrix}$$ 
We omitted $U + V + W$ here as it is linearly dependent on $A, B, C, D$:
$$ U + V + W = A + B + C - D$$
Now you want to compute
$$ \begin{align} 
&~ \Pr\{A > 0, B > 0, C > 0|D > 0, A + B + C - D > 0\} \\
= &~ \frac { \Pr\{A > 0, B > 0, C > 0, D > 0, A + B + C - D > 0\} } {\Pr\{X + Y + Z > 0, U + V + W > 0\}} \\
= &~ 4\Pr\{A > 0, B > 0, C > 0, D > 0, A + B + C - D > 0\} 
\end{align}$$
By using Cholesky decomposition on the covariance matrix (actually it is easy to observe too), we know that
$$ \begin{bmatrix}
A \\ B \\ C \\ D 
\end{bmatrix} 
\stackrel {d} {=}
\sqrt{2}\sigma \begin{bmatrix}
T_1 \\  T_2  \\  T_3  \\ 
\displaystyle \frac {1} {2} (T_1 + T_2 + T_3) + \frac {\sqrt{3}} {2} T_4
\end{bmatrix}$$
where $(T_1, T_2, T_3, T_4)$ are i.i.d. standard normal. So the required probability becomes
$$ 4\Pr\{T_1 > 0, T_2 > 0, T_3 > 0, T_1 + T_2 + T_3 + \sqrt{3}T_4 > 0,
T_1 + T_2 + T_3 - \sqrt{3}T_4 > 0\}$$
