I would like to take double integral over standard normal distribution but it does not work.

1 次查看（过去 30 天）
Donghwan Hyun
Donghwan Hyun 2020-9-24
评论： Donghwan Hyun 2020-9-24
Hi everyone, I'm trying to take double integral over the pdf to find out the probablity,
However using the code below, it just took the time and memory(over 9GB) without any error message, but it did not give any solution.
What I'm trying: where .
syms z_1 z_2 rho
phi(z_1,z_2,rho) = (exp(-(((z_1)^2)+((z_2)^2)-(2*rho*(z_1)*(z_2)))/(2*(1-rho^2)))/(2*pi*sqrt(1-rho^2)))
phi(z_1, z_2, rho) =
exp((z_1^2 - 2*rho*z_1*z_2 + z_2^2)/(2*rho^2 - 2))/(2*pi*(1 - rho^2)^(1/2))
syms mu_1 mu_2 sigma lambda mu_M
Pr1 = int(int(phi,z_2,[-inf ((mu_M-mu_2)/sigma)]),z_1,[-inf ((lambda-mu_1)/sigma)])
I would appreciate your help.
  6 个评论
显示 4更早的评论
Bjorn Gustavsson
Bjorn Gustavsson 2020-9-24
@Bruno, For the uncorrelated case things behave nicely:
>> syms x z xlim zlim real
>> TwoDtrivial = int(int(exp(-x^2-z^2),'x',-inf,xlim),'z',-inf,zlim)
TwoDtrivial =
(pi*(erf(xlim) + 1)*(erf(zlim) + 1))/4
For the case I sketched below (in its simplest version where the transformation is just a 45-degree rotation) I end up with:
>> int(int(exp(-x^2-z^2),'x',z-zlim,zlim-z),'z',-inf,zlim)
ans =
int(-pi^(1/2)*erf(z - zlim)*exp(-z^2), z == -Inf..zlim)
Not a solution that doesn't require an integration - but at least something like a solution...

请先登录，再进行评论。

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

采纳的回答

Bjorn Gustavsson
Bjorn Gustavsson 2020-9-24
Your problem might be that the symbolic toolbox have no way of knowing that your correlation is between -1 and 1. If you try the trivial calculation for the case where there is no correlation it is unproblematically fast. One way forward might be to transform your coordinates to one where there is no correlation. This would have to be combined with a change in the integration-limits, but by a simple sketch this seems "reasonably" manageable to do - the rectangular region in your original variables ought to become a triangular region in the transformed variables.
HTH

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Number Theory 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by