prob.m

用于计算贝叶斯频谱概率

function [sigma,stloge]=prob(G,d,ino,ifun);
% new test for eigenvalues/eigenfunctions
zloge=10e6;%
%number of functions used
for k=1:ifun
    for ell=1:ifun
        g(k,ell)=0;
        for i=1:ino
            %This is matrix called gij 
            %by Bretthorst.
            %It is NOT diagonal. Equation 3.4
            g(k,ell)=g(k,ell)+G(k,i)*G(ell,i);
        end
    end
end
%correct to here
[X,D]=eig(g);
%maybe to here
for i=1:ino
     for k=1:ifun
         H(k,i)=0;
     end
 end
for i=1:ino
    for j=1:ifun;   %j is the eigenvalue index
        for k=1:ifun
            H(j,i)=H(j,i)+(1/sqrt(D(j,j)))*G(k,i)*X(k,j);
        end
    end
end
H2=0;
for j=1:ifun
    H1=0;
    for ell=1:ino
        H1=H1+d(ell)*H(j,ell);
        %eqn 3.13
    end
   HI(j)=H1;           % may be excessive use.. never uses again******WHY??
    H2=H2+H1*H1;
end
H2BAR=H2/ifun;%eqn 3.15
%our d2 is Bretthorst's Y2
d2=0;
for i=1:ino
    d2=d2+d(i)*d(i);
end
d2=d2/ino;
qq=1-ifun*H2BAR/(ino*d2);
%eqn 3.17 except for power
%stloge is student t distribution
stloge=log(qq)*((ifun-ino)/2);
%the above is 3.17 in log form
%stloge is student t distribution's logarithm
ahold=stloge-zloge;
st=0;
if(abs(zloge)~=0) st=exp(ahold);
end
    sigma=sqrt((ino/(ino-ifun-2))*(d2-ifun*H2BAR/ino));
    %eqn 3.22
phat=ifun*H2BAR*st
;     %eqn 4.15 on P57
%end 
%cases
%Bretthorst used prob as a subroutine, going back to it from 
%where parameter was set