# Problem plotting expression involving Generalized hypergeometric functions $_2F_2 \left(.,.,. \right)$

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. I've got the same plotting issue with Matlab.

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. Perhaps there is another way to derive the above equation with other functions or there might be a way to simplify it and get more accurate results.

• In the last case, simply changing 12.85 to 11.85 causes a massive change in the expectation and produces a large negative value. I suspect a numerical issue. Using NExpectation gives a better result: With[{κ = 12.85, θ = 0.336551, a = 3, b = 1/5}, G = GammaDistribution[κ, θ]; a * NExpectation[q[Sqrt[b] y], y \[Distributed] G] ] (* 0.127737 *) – flinty Jul 28 '20 at 13:26
• It's impolite of you to expect an answer without providing the code to reproduce. – Mariusz Iwaniuk Jul 28 '20 at 14:11

## 1 Answer

It can be solved by increasing the precision with SetPrecision. See the excerpt below. This question has also been solved here: https://mathoverflow.net/questions/366754/inaccurate-results-for-the-analytical-expression-of-mathbbe-left-a-mathcal/366791#366798

prec = 100;
m = 16;
a = 4 (1 - (1/Sqrt[m]));
b = 3/(m - 1);

pdb = -12;
p = 10^(pdb/10);

betag = (1.5)^2;

betah = (2.3)^2;

bits = 8;

Q = 2^bits;

n = 8;

ef = SetPrecision[
n p (betag  betah) (1 + ((n - 1)/16) (Q^2) (Sin[Pi/Q]^2)), prec];

ef2 = SetPrecision[
n (p^2) (1/
256) (((betag^2) (betah^2))) (512 (n +
1) + (32 (n - 1) (Q^2))/(Pi^2) + ((n -
1) (Q^2) (-32 Cos[(4 Pi)/Q] +
Pi (Sin[
Pi/Q]^2) (16 (Pi + 4  n  Pi) + (n -
2) Q ((n - 3) Pi Q (Sin[Pi/Q]^2) +
16 Sin[(2 Pi)/Q]))))/(Pi^2)), prec];

a1 = SetPrecision[(ef2 - (ef^2)), prec];
a2 = SetPrecision[(ef2 - 5 (ef^2)), prec];
a3 = SetPrecision[-6 (ef^2), prec];

κ = SetPrecision[(-a2 + Sqrt[(a2^2) - 4 a1 a3])/(2 a1), prec];

θ = SetPrecision[Sqrt[ef/(κ (κ + 1))], prec];

f1 = SetPrecision[
2^(-1 - κ/2) a b^(-κ/2) θ^-κ, prec];
f2 = SetPrecision[(κ HypergeometricPFQ[{1/2 + κ/2,
1 + κ/2}, {3/2, 3/2 + κ/2}, 1/(2 b θ^2)])/(
Sqrt Sqrt[b] θ Gamma[(3 + κ)/2]), prec];
f3 = SetPrecision[
HypergeometricPFQ[{1/2 + κ/2, κ/2}, {1/2,
1 + κ/2}, 1/(2 b θ^2)]/Gamma[1 + κ/2], prec];

SetPrecision[f1 (f3 - f2), prec]