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_2\left(\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}+1\right)}-\frac{\kappa \, _2F_2\left(\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.

mightbe a bug in Mathematica kernel (I have already found one such issue with the $_2F_1$ function — Wolfram corrected it within a few months after I contacted them), but I suspect simple rounding error here: you seem to be subtractinghugevalues here, the hypergeometric function is of the order $10^11$ to $10^16$ for the given problematic parameters. $\endgroup$ – Mateusz Kwaśnicki Jul 28 at 7:33`SetPrecision[12.85, 50]`

etc. for all numeric constants. $\endgroup$ – Mateusz Kwaśnicki Jul 28 at 12:40