用matlab怎么解带有分布函数的方程

momo67799 发表于 2012-4-11 09:41

对于5和6，可以考虑变量替换，如x3=x1+x2，这样的话就可以正常使用mvncdf函数了，但它又带来了一个新问题，新的多元变量为[x1 x3]，此时协方差矩阵要重新给出。积分是个十分费力的事。

mj2012 发表于 2012-4-11 11:09
对于5和6，可以考虑变量替换，如x3=x1+x2，这样的话就可以正常使用mvncdf函数了，但它又带来了一个新问题 ...

我用二分查找法解x1和x2，涉及到积分的式子我试过，是正确的，但是放在循环语句里，就运行错误，麻烦您帮我看一下是哪里出现了错误，谢谢！
ZERO=10^-4;
ZERO1=10^-8;
a=150;
mu1=400;sigma1=112;
mu2=300;sigma2=105;
x2_min=0;x2_max=800;
Error=1;
     while(Error&&(x2_max-x2_min)>ZERO)
     x2=(x2_min+x2_max)/2;
     x1_min=0;x1_max=800;
     Error1=1;
        while(Error1&&(x1_max-x1_min)>ZERO1)
        x1=(x1_min+x1_max)/2;
        syms x y
        P1=vpa(int(int(mvnpdf([x y],[mu1 mu2],[sigma1^2 0;0 sigma2^2]),y,x1+x2+a-x,+inf),x,x1,x1+a));
        Phi1=-40*normcdf(x1,mu1,sigma1)+20*( normcdf(x1+a,mu1,sigma1)-normcdf(x1,mu1,sigma1))+45*(1-normcdf(x1+a,mu1,sigma1))+20*P1;  
            if Phi1<-ZERO1
                 x1_max=x1;
                 elseif Phi1>ZERO1
                 x1_min=x1;
                 else
                 Error1=0;
            end
        end
        x1
        syms x y
        P2=vpa(int(int(mvnpdf([x y],[mu1 mu2],[sigma1^2 0;0 sigma2^2]),y,x2,+inf),x,x1+a,+inf));
        P3=vpa(int(int(mvnpdf([x y],[mu1 mu2],[sigma1^2 0;0 sigma2^2]),y,x2+a,+inf),x,0,x1));
        P4=vpa(int(int(mvnpdf([x y],[mu1 mu2],[sigma1^2 0;0 sigma2^2]),y,x2,x1+x2+a-x),x,x1,x1+a));
        P5=vpa(int(int(mvnpdf([x y],[mu1 mu2],[sigma1^2 0;0 sigma2^2]),y,x2,x2+a),x,0,x1));
        Phi2=30*normcdf(x2,mu2,sigma2)+90*P2+90*P1+70*P4+90*P3+70*P5-55;
        if Phi2<-ZERO
        x2_max=x2;
        elseif Phi2>ZERO
        x2_min=x2;
        else
        Error=0;
        end
     end
x1
x2