I have a random variable $X$ which follows the following probability density function,
$$ p(x) = \frac{1}{4\pi} \Big[ \operatorname{erf}\Big(\frac{k\mu-x+2\pi}{\sqrt{2}k\sigma}\Big) - \operatorname{erf}\Big(\frac{k\mu-x}{\sqrt{2}k\sigma}\Big) \Big]$$
I wanted to find the probability density function for the function of the random variable $Y = \cos(X)$.
I followed the CDF method to find this.
$$ F_Y(y) = p(Y \leq y) = p(\cos(X) \leq y) = p(X \leq \cos^{-1}(y))$$
$$ F(y) = \int_{-\infty}^{\cos^{-1}(y)} \frac{1}{4\pi} \Big[ \operatorname{erf}\Big(\frac{k\mu-x+2\pi}{\sqrt{2}k\sigma}\Big) - \operatorname{erf}\Big(\frac{k\mu-x}{\sqrt{2}k\sigma}\Big) \Big] dx $$
The terms $\mu$, $\sigma$ and $k$ are constants.
$$ k \quad \epsilon \quad Z; k \geq 0$$ $$ \sigma \quad \epsilon \quad R; \sigma > 0 $$
Mathematica Output: $$ F_Y(y) = \frac{1}{8 \pi ^{3/2}} \Bigg[ -\sqrt{\pi } \left(2 i \log \left(\sqrt{1-y^2}+i y\right)-2 k \mu+\pi \right) \text{erf}\left(\frac{-2 i \log \left(\sqrt{1-y^2}+i y\right)+2 k \mu-\pi }{2 \sqrt{2} k \sigma}\right) \\ +\sqrt{\pi } \left(2 i \log \left(\sqrt{1-y^2}+i y\right) -2 k \mu-3 \pi \right) \text{erf}\left(\frac{-2 i \log \left(\sqrt{1-y^2}+i y\right) +2 k \mu+3 \pi }{2 \sqrt{2} k \sigma}\right)- \\ 2 \sqrt{2} k \sigma \left(\sqrt{1-y^2}+i y\right)^{\frac{i (2 k \mu+3 \pi )}{2 k^2 \sigma^2}} \exp \left(-\frac{(2 k \mu+3 \pi )^2-4 \log ^2\left(\sqrt{1-y^2}+i y\right)}{8 k^2 \sigma^2}\right)+ \\ 2 \sqrt{2} k \sigma \left(\sqrt{1-y^2}+i y\right)^{-\frac{i (\pi -2 k \mu)}{2 k^2 \sigma^2}} \exp \left(-\frac{(\pi -2 k \mu)^2-4 \log ^2\left(\sqrt{1-y^2}+i y\right)}{8 k^2 \sigma^2}\right)+4 \pi ^{3/2} \Bigg] $$
When I do this integral, it gives me a complex valued function on mathematica. How can this be possible that the CDF is complex valued? If it is so, what should I consider as the CDF? Should I consider the real part only?
============== EDIT ==========
How I found $p(x)$ ?
It comes from another problem of finding the PDF of a random variable $ X = k U + \beta $, where k is an integer. $U$ comes from a Gaussian distribution and $\beta$ comes from a uniform distribution.
$$ U \sim \mathcal{N}(\mu, \sigma) $$ and $$ \beta \sim \mathcal{U}(0, 2\pi) $$
As it is an addition, I found the PDF of $X$ as the convolution of the other 2 PDFs.
$$ p(x) = \int_{0}^{2\pi} \frac{1}{2\pi \sqrt{2 \pi k^2 \sigma^2}} e^{-(kx - k\mu - \tau)/(2k^2\sigma^2)} d\tau $$
$$ p(x) = \frac{1}{4\pi} \Big[ \operatorname{erf}\Big(\frac{k\mu-x+2\pi}{\sqrt{2}k\sigma}\Big) - \operatorname{erf}\Big(\frac{k\mu-x}{\sqrt{2}k\sigma}\Big) \Big] $$
I thought this was correct, because when I integrated $p(x)$ from $-\infty$ to $+\infty$, I got $1$.