I’m trying to plot a graph for the following expectation

$$mathbb{E}left( a mathcal{Q} left( sqrt{b } gamma right) right)=a 2^{-frac{kappa }{2}-1} b^{-frac{kappa }{2}} theta ^{-kappa } left(frac{, _2F_2left(frac{kappa }{2}+frac{1}{2},frac{kappa }{2};frac{1}{2},frac{kappa }{2}+1;frac{1}{2 b theta ^2}right)}{Gamma left(frac{kappa }{2}+1right)}-frac{kappa , _2F_2left(frac{kappa }{2}+frac{1}{2},frac{kappa }{2}+1;frac{3}{2},frac{kappa }{2}+frac{3}{2};frac{1}{2 b theta ^2}right)}{sqrt{2} sqrt{b} theta Gamma left(frac{kappa +3}{2}right)}right)$$

where $a$ and $b$ are constant values, $mathcal{Q}$ is the Gaussian Q-function, which is defined as $mathcal{Q}(x) = frac{1}{sqrt{2 pi}}int_{x}^{infty} e^{-u^2/2}du$ and $gamma$ is a random variable with Gamma distribition, i.e., $f_{gamma}(y) sim frac{1}{Gamma(kappa)theta^{kappa}} y^{kappa-1} e^{-y/theta} $ with $kappa > 0$ and $theta > 0$.

This equation was also found with Mathematica, so it seems to be correct.

Follows some examples, where I have checked the analytical results against the simulated ones.

When $kappa = 12.85$, $theta = 0.533397$, $a=3$ and $b = 1/5$ it returns the correct value $0.0218116$.

When $kappa = 12.85$, $theta = 0.475391$, $a=3$ and $b = 1/5$ it returns the correct value $0.0408816$.

When $kappa = 12.85$, $theta = 0.423692$, $a=3$ and $b = 1/5$ it returns the value $-1.49831$, which is negative. However, the correct result should be a value around $0.0585$.

When $kappa = 12.85$, $theta = 0.336551$, $a=3$ and $b = 1/5$ it returns the value $630902$. However, the correct result should be a value around $0.1277$.

Therefore, the issue happens as $theta$ decreases. For values of $theta > 0.423692$ the analytical matches the simulated results. The issue only happens when $theta <= 0.423692$.

I’d like to know if that is an accuracy issue or if I’m missing something here and if there is a way to correctly plot a graph that matches the simulation.