Given the Likelihood $X \sim{\cal N}(\mu (\theta ),\Sigma )$,and $\theta$ uniform rv,how do i compute the joint distribution $f({X})$ analytically? I know the likelihood function to be $X \sim{\cal N}(\mu (\theta ),\Sigma )$ with $X$ a vector of two random normal variables $X = [{X_1},{X_2}]$, with known var-cov matrix ${\Sigma }$ which is diagonal and ins't dependent on $\theta$ and the vector of means which are dependent on $\theta$: 
$$\mu (\theta ) = [{\mu _1}(\theta ),{\mu _2}(\theta )]$$
with $$\begin{array}{l}
{\mu _1}(\theta ) = d\sin \theta \\
{\mu _2}(\theta ) = d\cos \theta 
\end{array}$$
with $d$ a constant, and $\theta$ uniformly distributed: $\theta  \sim U\left( {\frac{{ - \pi }}{2},\frac{{  \pi }}{2}} \right)$.
I want to compute the joint distribution $f({x_1},{x_2})$ (which is the function in the denominator of posterior distribution (e.g. in Bayesian estimator)
$$f({x_1},{x_2}) = \int\limits_{ - \pi /2}^{\pi /2} f ({x_1},{x_2}|\theta )f(\theta )d\theta$$
How do I calculate the integral analytically?
(PS I need to find this function not to find the Bayesian estimator)
Thank you in advance
 A: This particular problem leads to an integral that has no analytic solution, and so you need to use numerical methods to evaluate it.  I will consider a slightly generalised version of your problem where we initially allow the variance matrix to be non-diagonal (i.e., we allow correlation between the variables), but I'll also show what happens when you assume the variance is diagonal.  To present this generalisation I will let $\sigma_1$ and $\sigma_2$ denote the standard deviations of the observed variables and I will let $-1<\rho<1$ be the correlation between these values.
Allowing initially for the above generalisation, the squared Mahalanobis distance in this problem is:
$$\begin{aligned}
D^2(\theta) 
&= (\mathbf{x} - \boldsymbol{\mu}(\theta))^\text{T} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}(\theta)) \\[6pt]
&= d^2 
\begin{bmatrix} \sin \theta & \cos \theta \end{bmatrix} 
\begin{bmatrix} \sigma_1^2 & \rho \sigma_1 \sigma_2 \\ \rho \sigma_1 \sigma_2 & \sigma_2^2 \end{bmatrix}^{-1}
\begin{bmatrix} \sin \theta \\ \cos \theta \end{bmatrix} \\[6pt]
&= \frac{d^2}{(1-\rho)^2 \sigma_1^2 \sigma_2^2} \cdot 
\begin{bmatrix} \sin \theta & \cos \theta \end{bmatrix} 
\begin{bmatrix} \sigma_2^2 & -\rho \sigma_1 \sigma_2 \\ -\rho \sigma_1 \sigma_2 & \sigma_1^2 \end{bmatrix}
\begin{bmatrix} \sin \theta \\ \cos \theta \end{bmatrix} \\[6pt]
&= \frac{d^2}{(1-\rho)^2 \sigma_1^2 \sigma_2^2} \cdot 
\Big( \sigma_2^2 \sin^2 \theta - 2 \rho \sigma_1 \sigma_2 \sin \theta \cos \theta + \sigma_1^2 \cos^2 \theta \Big). \\[6pt]
\end{aligned}$$
Thus, the integral you want to compute is of the form:
$$\int \exp \Big( A \cdot \sin^2 \theta + B \cdot \sin \theta \cos \theta + C \cdot \cos^2 \theta \Big) \ d \theta,$$
which does not have a general analytic solution.  In the case where $\rho = 0$ the variance matrix is diagonal and the above integral simplifies by taking $B=0$.  The only time that the integral has an analytic solution is when $\rho = 0$ and $\sigma_1 = \sigma_2 = \sigma$, in which case the squared Mahalanobis distance simplifies to:
$$D^2(\theta) = \frac{d^2}{\sigma^2} \cdot 
\Big( \sin^2 \theta + \cos^2 \theta \Big) = \frac{d^2}{\sigma^2}.$$
As you can see, the integral you want to compute is the exponential of a "quadratic" trigonometric function.  This will require numerical methods.
