%TE-MAP均衡
%copyright:flyingxiang
function  [Lec,saidaequa]=Lblocoequaiter(y,var,C1equa,C0equa)
%得到均衡的輸出和外部信息

S=4;   %狀態數 
P=zeros(4,4,length(y));
oute=[-1.629 -.815 0 0;0 0 .001 .815;-.815 -.001 0 0;0 0 .815 1.629]; %對應網格圖的輸出矩陣

A1=[0 1 0 0;0 0 0 1;0 1 0 0;0 0 0 1];
A0=[1 0 0 0;0 0 1 0;1 0 0 0;0 0 1 0];
pos=[1 1 0 0;0 0 1 1;1 1 0 0;0 0 1 1];
gamma=zeros(S,S,length(y));
B1=zeros(S,S,length(y));
B0=zeros(S,S,length(y));
num=zeros(1,length(y));
den=zeros(1,length(y));
priori_equa=zeros(S,S,length(y));

%生成對應的先驗信息
for k=1:length(y)
    priori_equa(:,:,k)=[C0equa(k),C1equa(k),0,0;0,0,C0equa(k),C1equa(k);C0equa(k),C1equa(k),0,0;0,0,C0equa(k),C1equa(k)];
    
end


%計算gamma
for k=1:length(y)   
   for i=1:4
      for j=1:4
         arg=-(y(k)-oute(i,j))^2;
         P(i,j,k)=priori_equa(i,j,k)*exp(arg/(2*var))*(1/sqrt(2*pi*var));
       
           end
       end
   end
gamma=P;

N=length(y)+1;

%初始化前向遞規向量f
f=zeros(S,1,N);
f(:,1,1)=[1;0;0;0] ;
c=zeros(1,N);
c(1)=1;
%遞規計算前向向量f
for i=1:(N-1)
   f(:,1,i+1)=gamma(:,:,i)'*f(:,1,i);
   c(i+1)=sum(f(:,1,i+1));
   f(:,1,i+1)=f(:,1,i+1)./c(i+1);
end

%初始化后向遞規向量b
b=zeros(S,1,N);
b(:,1,N)=[1;1;1;1];
%遞規計算后向向量b
for i=(N-1):-1:1
   b(:,1,i)=gamma(:,:,i)*b(:,1,i+1);
   cb(i)=sum(b(:,1,i));
   b(:,1,i)=b(:,1,i)./cb(i);

end
b;


for n=1:length(y)
B1(:,:,n)=A1.*gamma(:,:,n);
B0(:,:,n)=A0.*gamma(:,:,n);
end

 
 
 L=zeros(1,N-1);
 saida=zeros(1,N-1);

for k=1:(N-1)
    
   num(k)=f(:,:,k)'*B1(:,:,k)*b(:,:,k+1);
    
   den(k)=f(:,:,k)'*B0(:,:,k)*b(:,:,k+1);
      Lec(k)=log((num(k))/(den(k)));
end
 
for i=1:length(y)
   if Lec(i)>=0
      saidaequa(i)=1;
   elseif Lec(i)<0
      saidaequa(i)=-1;
   end
end