Integrals involving symbolic probability distributions

Question

Aaron J. Hendrickson 2017-7-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/349135-integrals-involving-symbolic-probability-distributions

评论： Walter Roberson 2017-7-18

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2017-7-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/349135-integrals-involving-symbolic-probability-distributions#answer_274497

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)

5 个评论
显示 3更早的评论隐藏 3更早的评论

Walter Roberson 2017-7-18

在 MATLAB Online 中打开

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 2017-7-18

Thank you.

Walter Roberson 2017-7-18

在 MATLAB Online 中打开

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)