Integrals involving symbolic probability distributions

I have the following script:

a1 = 49.5;
a2 = 49.5;
b1 = 153.4787/a1;
b2 = 40.1/a2;

syms x p
assume(p<=0)
fph_n = symfun(int(gampdf(x,a1,b1)*gampdf(-p+x,a2,b2),x,0,Inf), p);

However, I keep getting MuPad errors and not sure why. Can anyone explain?

Answer 1

[Walter Roberson]

Unfortunately gampdf is only defined for numeric x; it is not part of the Symbolic toolbox.

The Symbolic Toolbox does have a symbolic equivalent, but it does not have a nice interface to it. See https://www.mathworks.com/help/symbolic/mupad_ref/stats-gammapdf.html and notice this is MuPAD, so to get at it from MATLAB you would need to use feval(symengine) or evalin(symengine)

[Walter Roberson]

My tests suggest that if you change the 153.4787 by +-/5 * 10^-5 (that is, you assume 153.4787 might be a rounded figure rather than 1534787/10000 exactly) then the answer can change in the third significant figure.

If you call the potential change in the 153.4787 as delta_b1 then the formula comes out as

(1817928720888968425241057364445160534589625558323953650711543669959064474728153250319263125155461001289698034825787330596840738400622283462809672291356658793406778632642350981205012154144160933954124098219184541017148077030935754728920144227491758223360000000000000000000000000000000000000/20671005336345936511218405359554688582335560687949793)*115^(1/2) * (-p)^(99/2) * exp((495/401)*p) * MeijerG([[0], []], [[97/2, -99/2], []], -(99/802) * (22770*delta_b1+4407787) * p / (349471+2277*delta_b1)) / (((22770*delta_b1+4407787)^97)^(1/2) * Pi^(3/2) * (349471+2277*delta_b1))

https://www.mathworks.com/matlabcentral/answers/296400-meijerg-function-using-matlab-s-symbolic-toolbox

[Aaron J. Hendrickson]

Can u please explain how you got that solution. The parameters I used were for example. I would want to be able to solve for arbitrary a's and b's.

[Walter Roberson]

The way I got the solution was that I went into Maple, and I assigned the specific values to the variables you indicated, and then did roughly

with(Statistics):
dist1 := u->PDF(RandomVariable(GammaDistribution(b1, a1)), u):
dist2 := u->PDF(RandomVariable(GammaDistribution(b2, a2)), u):
simplify(int(dist1(x)*dist2(x-p), x = 0 .. infinity)) assuming p<0, a1>=0, a2>=0, b1>=0, b2>=0;
simplify(convert(%,MeijerG)) assuming p<0, a1>=0, a2>=0, b1>=0, b2>=0;

With specific numeric values not assigned, the result corresponds to

-3 * sin(pi * a2) * exp(-p / b1) * (-gamma(a2 + a1 + 3) ^ 2 * gamma(1 - a2) * (b1 + b2) ^ (-a2 - a1 - 1) * p ^ 2 * (sin(pi * a1) * cos(pi * a2) + cos(pi * a1) * sin(pi * a2)) * (b1 ^ (a2 - 1) * b2 ^ (2 + a1) + 3 * b1 ^ (a2 + 1) * b2 ^ a1 + b1 ^ (2 + a2) * b2 ^ (a1 - 1) + 3 * b2 ^ (1 + a1) * b1 ^ a2) * hypergeom([1 - a1], [2 - a1 - a2], (b1 + b2) * p / b1 / b2) / 3 + pi * (a2 + a1 + 2) * (a2 * ((2 ./ 3 * (b2 ^ (-a2) * b1 ^ (-a1) + b1 ^ (1 - a1) * b2 ^ (-1 - a2) / 2 + b1 ^ (-a1 - 1) * b2 ^ (1 - a2) / 2) * (a1 + a2) * p ^ (2 + a1 / 2) + p ^ (3 + a1 / 2) * (b2 ^ (-a2) * b1 ^ (-a1 - 1) + b1 ^ (1 - a1) * b2 ^ (-2 - a2) / 3 + b1 ^ (-a1 - 2) * b2 ^ (1 - a2) / 3 + b2 ^ (-1 - a2) * b1 ^ (-a1))) * (a1 + a2) * (-1 ./ (p)) ^ (-a2) + ((-(-p) ^ (2 + a2) * (a1 + a2) * b2 ^ (-1 - a2) / 3 + b2 ^ (-2 - a2) * (-p) ^ (a2 + 3) / 3) * b1 ^ (1 - a1) + (-(-p) ^ (2 + a2) * (a1 + a2) * b1 ^ (-a1 - 1) / 3 + b1 ^ (-a1 - 2) * (-p) ^ (a2 + 3) / 3) * b2 ^ (1 - a2) + b2 ^ (-a2) * b1 ^ (-a1 - 1) * (-p) ^ (a2 + 3) - 2 ./ 3 * b1 ^ (-a1) * (-3 ./ 2 * b2 ^ (-1 - a2) * (-p) ^ (a2 + 3) + b2 ^ (-a2) * (-p) ^ (2 + a2) * (a1 + a2))) * p ^ (a1 / 2)) * hypergeom([a2 + 1], [a2 + a1 + 3], (b1 + b2) * p / b1 / b2) - 4 ./ 3 * (-(2 * a2 * (b2 ^ (-a2) * b1 ^ (-a1) + b1 ^ (1 - a1) * b2 ^ (-1 - a2) / 2 + b1 ^ (-a1 - 1) * b2 ^ (1 - a2) / 2) * p ^ (2 + a1 / 2) + p ^ (1 + a1 / 2) * (a2 + a1 + 2) * (a1 - 1 + a2) * (b2 ^ (-a2) * b1 ^ (1 - a1) + b2 ^ (1 - a2) * b1 ^ (-a1))) * (a2 + a1 + 1) * (a1 + a2) * (-1 ./ (p)) ^ (-a2) / 4 + p ^ (a1 / 2) * (-p) ^ (2 + a2) * a2 * (b2 ^ (-a2) * b1 ^ (-a1) + b1 ^ (1 - a1) * b2 ^ (-1 - a2) / 2 + b1 ^ (-a1 - 1) * b2 ^ (1 - a2) / 2)) * hypergeom([a2], [a2 + a1 + 2], (b1 + b2) * p / b1 / b2)) * (1 ./ (p)) ^ (-a1 / 2) * (a2 + a1 + 1) * gamma(a1) * (a1 + a2)) / p ^ 2 / gamma(a2 + a1 + 3) / sin(pi * (a1 + a2)) / pi / (b1 + b2) / (a1 - 1 + a2) / (a1 + a2) / (a2 + a1 + 1) / (a2 + a1 + 2) / gamma(a1)

[Aaron J. Hendrickson]

Thank you.

[Walter Roberson]

symLIST = @(varargin)feval(symengine,'DOM_LIST',varargin{:});
symRANGE = @(a,b) feval(symengine, '_range', a, b);

syms a1 b1 a2 b2 x p
dist1 = feval(symengine, 'stats::gammaPDF', a1, b1);
dist2 = feval(symengine, 'stats::gammaPDF', a2, b2);
int(feval(symengine, dist1, x) .* feval(symengine, dist2, -p+x), x, 0, inf)

Unfortunately it leaves it unevaluated unless you give specific numeric a1, b1, a2, b2, and p.