Sampling 6D Relative Poses from a Normal Distribution I'm fairly new to statistical mathematics, and needed to solve a problem that may rely heavily on it. 
Context: I'm working on 7D robot poses - i.e. 3 translation dimensions and 4 rotations dimensions (quaternion) constitute a single vector. A kernel $k(T_p^q{'}, T_p^q)$  accepts as input a reference pose and a target pose, and returns a single scalar value indicating the "similarity" of the target pose to the set of reference poses. Here, $T_p^q$ is the target pose whereas $T_p^q{'}$ is a reference pose from the set of reference poses. This is expected to be a multi-modal Normal Distribution. Note that $T_p^q$ is the pose of object p with reference to object q. There may be 2 or more objects present in the scene. I want to visualize the output of the kernel against a number of parameters and inputs.
Task: For every pair of objects in the scene, generate 7D relative poses that correspond to a given normal distribution. These generated poses will become the reference set of poses.
Approach: For the sake of ease, I intend to sample 6D poses and convert them to 7D ones (Euler to Quaternion transformations). It seems simple enough to directly generate 6 values sampled from a desired univariate normal distribution and mash them together into a single vector to generate 6D poses in the world reference frame. However, the kernel $k$ works on relative poses. Therefore, before processing through $k$, the poses are all converted into relative poses.
Problem: I have no clue if and how the transformation of poses from one reference frame to another affects the underlying distribution of the poses. An example-
Given the set of objects $$\mathcal{O}=\{o_1,...,o_M\}$$
The set of reference poses was originally a set of object poses in the world reference frame $W$-$$R=\{T_{o_1}^W, ..., T_{o_M}^W; R_{o_i}^W \sim \mathcal{N}_6(\mu, \sigma^2) \forall o_i \in \mathcal{O}\}$$ which were generated as per the method defined above. Now, if I use the standard pose transformation mathematics to generate all the reference poses,
$$R'=\{R_{o_1},...,R_{o_M}\},$$
$$R_{o_1}=\{T_{o_1}^{o_2},...,T_{o_1}^{o_M}\},\ and\ so\ on$$
will these reference poses still follow a normal distribution? If so, what are the parameters of this normal distribution? If not, how do I ensure that all poses in a given set of relative reference poses, say, $R_{o_1}$, follow one normal distribution, and so do the poses in $R_{o_2}$, and so on and so forth? Or even have all relative poses follow the same normal distribution, as was originally the case in $R$?
I have perused through a bunch of literature on sampling from normal distributions, as well as quaternion transforms and relative pose computation, but this unholy combination of all three seems to be missing. And maybe I'm not skilled enough yet to see how to combine it all myself.
Thank you!
 A: It seems like you should be sampling from a vector space, that way you can properly define your normal distribution.  This is very common in defining noise and uncertainty about state estimates in kalman Filtering.
What I would do is use the lie algebra associated with the SE(3) lie group (homogenous transforms).  Define the mean transform, $T$, then sample "errors" to that mean as 6-vectors normally about (0,0,0,0,0,0).  To convert an "error" sample into a transform, use the exponential map which converts the algebra to the lie group:
$$
T_i = T \exp\left(\begin{bmatrix} (\omega)^\times & v \\ 0 &0 \end{bmatrix}\right)
$$
 Where your sample in $\mathbb{R}^6$ is partitioned as 
$$
\xi = \begin{bmatrix} \omega \\ v \end{bmatrix}
$$
And $(\cdot)^\times$ is the skew-symmetric matrix
While I believe this very naturally solves the relative nature of your problem, I believe this is also much more principled than using Euler angles to generate your samples, because Euler angles don't form a vector space, and therefore, "normally distributed Euler angles" doesn't make any sense.
Consider, if you have a set of three Euler angles and you tweak pitch alone, you are actually changing the coordinate frame that defined roll.  The intermediate frames used in defining euler angles can cause a lot of headache when trying to consider correlations between samples.
Using the lie algebra is better because $\xi$ actually is a vector.  You can add and subtract $\xi$'s, consider their magnitude of deviation from the mean, calculate covariance, etc... which makes them a perfect basis for a multivariate normal.
In addition, the inverse of the above is a principled way to find the difference between two transforms.
