how to solve a problem about marginal

1 次查看（过去 30 天）

Skydriver 2019-4-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/457789-how-to-solve-a-problem-about-marginal

评论： Walter Roberson 2019-9-15

I have a coding :

h_data = [0; 0.1568;   0.2112;    0.264;  0.3403;   0.469;  0.58;	0.7];
% deltaPGA = h_data(2:8) - h_data(1:7);
Mwgrid = 4.35:0.3:9.15;
amaxgrid = linspace(0,3,36); %0:0.01:3;
grid = meshgrid(Mwgrid,amaxgrid);
PGAgrid = 0.19:0.02:0.39;
for i=1:length(h_data)-1
    %calculating Conditional Probability(Mw given PGA>h) in equation (5)    
    pdfMwPGAh = marginal(:,i);
    
    %calculating Pr(amax,PGA>h) in equation (9)   
    [PGA_temp,~] = hist(PGA,PGAgrid);
    pdfPGA = PGA_temp/sum(PGA_temp);
    
    %calculating joint probability between amax and PGA in equation (9)
    %p(amax,PGA)={1/amaxSigln(amax)(2phi)^1/2
    %exp[-{ln(amax)-Mln(amax)^2/2Sig^2ln(amax)]Deltaamax}p(PGA)
    pdfamaxPGA=zeros(length(amaxgrid),length(PGAgrid));
    for j=1:length(amaxgrid)
        for k=1:length(pdfPGA)
            if PGAgrid(k)>h_data(i+1)
                pdfamaxPGA(j,k)=lognpdf(amaxgrid(j),-0.15+(-0.13+1)*log(PGAgrid(k)),0.52).*0.01.*pdfPGA(k);
            end
        end
    end
    %Sum Joint Probability between amax and PGA in equation (8)
    pdfamaxPGAh = sum(pdfamaxPGA');  
    
    %determine conditional probability in equation (7)
    if sum(pdfamaxPGAh)~=0
        pdfamaxPGAh = pdfamaxPGAh/sum(pdfamaxPGAh);
    end
    
    %calculating Pr(amax,Mw,PGA>h) = Pr(amax,PGA>h)*Pr(Mw given
    %PGA>h)formulation (5)
    for j=1:length(amaxgrid)
        for k=1:length(pdfMwPGAh)
            pdfamaxMwPGAh(j,k,i)= pdfamaxPGAh(j)*pdfamaxPGAh(j) * pdfMwPGAh(k);
        end
    end
end