Conditional_Probability_Computation.m

if (~memoryless)
  if (rounds>1)
      disp('Conditional Probability Computation');
      % P(E) sum of every probability of every combination that exibit E
      dimensioni=length(SuppliersNew);
      P_co=coalition(length(SuppliersNew));
      [ri,co]=size(P_co);
      P_post=[];
      for i=1:1:ri
         aux_pos=P_co(i,:).*probability_old;
         aux_neg=(1-P_co(i,:)).*(1-probability_old);
         P_post=[P_post;(aux_pos+aux_neg)];
      end   
      P_Tutti=prod(probability_old,2);
      P_Nessuno=prod(1-probability_old,2);
      % global scenario.
      if (GlobalContext==1)
          P_Stable=0.8;
          P_Up=0.1;
          P_Down=0.1;
      else
        if (GlobalContext==2)
            P_Stable=0.3;
            P_Up=0.2;
            P_Down=0.5;
        else
            P_Stable=0.3;
            P_Up=0.5;
            P_Down=0.2;
        end
      end  
      %sum(prod(P_post,2).*P_Up) for market down
      %P_E=sum(prod(P_post,2).*P_Down)+(P_Tutti*P_Up)+(P_Nessuno*P_Down);
      %P_num=prod((((decision<probability_old).*probability_old)+((1-(decision<probability_old)).*(1-probability_old))),2).*P_Up;
      actors=[];
      actors.den=[];
      actors.num=[];
      actors.att=[];
      actors.Pa_e=[];
      for attori=1:1:length(SuppliersNew)
          % different probability table for different attacks
          if get(handles.radiobutton3,'Value')
         %    attfield='ProdCap congiunto
          Matrice=[MatProcCapCostMercStable;MatProcCapCostMercDown;MatProcCapCostMercUp];
          else
              if get(handles.radiobutton4,'Value')
                  %scenario= 'Lead Time congiunto
              else
                  if get(handles.radiobutton5,'Value')
                    %UnitProdCost
                    Matrice=[MatCostMercStable;MatCostMercDown;MatCostMercUp];
                  else
                     if get(handles.radiobutton6,'Value')
                       %ProdCap
                       Matrice=[MatProcCapMercStable;MatProcCapMercDown;MatProcCapMercUp];
                     else
                         if get(handles.radiobutton7,'Value')
                            %LeadTime
                         end          
                     end
                  end
              end
          end    
           %2 note each actor's demand (evidence)
           %3 find each branch for the evidence of interest without knowing
           % the market status. 
           %4 select branch
           if ((SuppliersNew(attori).Demand) < x_prec(attori))
               %if less then positive
               index_row=find(Matrice(:,dimensioni+attori)>0);
           else if ((SuppliersNew(attori).Demand) < x_prec(attori))
                   %if greatest then negative
                    index_row=find(Matrice(:,dimensioni+attori)<0);
               else
                   index_row=find(Matrice(:,dimensioni+attori)==0);
                   % equal
               end
           end    
          %5 Compute denominator 
          %6 numerator only one each time
          [rig,colo]=size(Matrice);
          rig=rig/3;
          P_den=sum([(prod(P_post(index_row<=rig,:),2).*P_Stable)',(prod(P_post((index_row((index_row>rig)&(index_row<=(rig*2)))-rig),:),2).*P_Down)',(prod(P_post((index_row(index_row>(rig*2))-(rig*2)),:),2).*P_Up)']);
          actors(attori).den=P_den;
          actors(attori).num=[(prod(P_post(index_row<=rig,:),2).*P_Stable)',(prod(P_post((index_row((index_row>rig)&(index_row<=(rig*2)))-rig),:),2).*P_Down)',(prod(P_post((index_row(index_row>(rig*2))-(rig*2)),:),2).*P_Up)'];
          actors(attori).att=[sum(Matrice(index_row,1:dimensioni)')];
          actors(attori).Pa_e=[actors(attori).num./actors(attori).den];
          %
          %P_num= prod(P_post,2).*P_Up;
          %P_a=P_num./P_E;
      
          % Sum up for number of attacks
          % P(k=0|E) none attack
          % P(k=1|E) only one
          % .... 
      
          Pk_e=[];
          %auxsum=(sum(P_co')');
          for k=1:1:(length(SuppliersNew)-1)
              Pk_e=[Pk_e;sum(actors(attori).Pa_e(actors(attori).att==k))];
          end  
          actors(attori).Pk_e=Pk_e;
                 
          % With the reactivity for every actor T compute P(T|k=?) for
          % every k
          load('react_default.mat')
          t=0:1/length(SuppliersNew):1;
          t2=t.^2;
          t3=t2.*t;
          xreact=react.ax*t3+react.bx*t2+react.cx*t+react.x0;
          yreact=react.ay*t3+react.by*t2+react.cy*t+react.y0;
          
          % P(T|E)= P(T|k=0)P(k=0|E) +
          % P(T|k=1)P(k=1|E)..... etc
          % La P(T|E) is P_Cond
          P_Cond(attori)=sum((yreact(2:dimensioni).*actors(attori).Pk_e'));
          probability(attori)=1-((1-P_Cond(attori))*(1-probability(attori)));
        end
  else
      P_Cond=1;
      probability=P_Cond.*probability;
  end
else
  P_Cond=1;    
end