Problems sampling from a $pdf$ over $SO\left(3\right)$

I have a probability density function over $SO\left(3\right)$, which I am trying to sample from. The $pdf$ is given as a generalized fourier series:

$$f\left(\omega,\theta,\phi\right)=\sum s_{\lambda}^{n}Z_{\lambda}^{n}\left(\omega,\theta,\phi\right)$$

where $s_{\lambda}^{n}$ are the coefficients and $Z_{\lambda}^{n}\left(\omega,\theta,\phi\right)$ are basis functions (symmetrized hyperspherical harmonics), and $SO\left(3\right)$ is parameterized by the variables $\left(\omega,\theta,\phi\right)$, which are the rotation angle, and spherical coordinates of the rotation axis, respectively. I want to sample values of $\left(\omega,\theta,\phi\right)$ (i.e. rotations) from this distribution, but I'm having trouble doing so.

So far I have essentially tried two methods: rejection sampling, a discrete method.

Both methods have given me something that is qualitatively similar to what I would expect, but they seem to have an erroneous uniform distribution superimposed. So I have two questions:

(1) Any ideas why I might be getting this uniform noise?

(2) Any suggestions of how to fix it or a better way to sample from this kind of distribution?

Further Details: To check my sampling method I used some software to generate a known distribution and sample rotations. I then computed $s_{\lambda}^{n}$ from these "correct samples". Then I tried to generate samples myself from the spectral form of $f\left(\omega,\theta,\phi\right)$ given above. Next I calculated $s_{\lambda}^{n}$ from my samples and compared the two. My samples led to the low order terms being generally too small in magnitude and the higher order terms being too large in magnitude.

The discrete method that I used consisted of generating a "grid" of points over $SO\left(3\right)$, and using these as the centers of bins. The bins were sized proportional to the probability density at the bin center. Then uniform samples were generated and the number that fell in each bin was then proportional to the associated probability of the respective bins.

Again, both methods produced a sort of background noise which looks like a uniform distribution on top of the correct distribution.

As, a side note, it seems like I ought to be able to exploit the form of this expression for efficient sampling, but I haven't made use of the spectral decomposition at all in my attempts.

-
as another side note, I tried setting the zero order term to 0, just for kicks, since the zero order term is the uniform portion, and I STILL got the uniform noise that I mentioned. – okj May 22 '13 at 19:21
In case you're not aware of it, a great method in general is inverse transform sampling. Presumably you don't have the analytic form for the inverse of the CDF corresponding to your PDF, but perhaps you can compute the CDF by numerically integrating your PDF, and then compute the inverse numerically as well. I have done this before to great effect. You might need to map your variables to one-dimension; I'm not sure how a multi-dimensional CDF is defined. – Douglas B. Staple May 31 '13 at 2:26
@DouglasB.Staple: I had thought of inverse transform sampling, it would definitely be the ideal method if it is possible. The problem is that in my case the domain of the function is the unit hypersphere $S^3$, so the random variables $\omega$, $\theta$, and $\phi$ are not independent. This means that I need to express the pdf as: $$f_{\Omega,\Theta,\Phi}\left(\omega,\theta,\phi\right)=f_{\Omega}\left(\omega \right)f_{\Theta\mid\Omega}\left(\theta\mid\omega\right)f_{\Phi\mid\Omega,\Theta‌​}\left(\phi\mid\Omega,\Theta\right)$$ How to get these numerically? – okj May 31 '13 at 18:58