function minusloglike = BF_Inspector
% BF_INSPECTOR runs through the steps of Bayesian Filtering for your inspection

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        PARAMETERS TO CHANGE                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initial Assumed Probability Mass function for State
% Sstart determined so that Ostart+Cstart+Sstart = 1
Ostart = 1;      % Assumed probability that State starts in Open
Cstart = 0;      % Assumed probability that State starts in Closed

% Transition Probabilities
% 1 = Open, 2 = Closed, 3 = Stuck,
% p11 is Probability of staying Open, 
% p12 is probablity of transitioning from Open to Closed
% Etc.
% p13, p23, p33 are not specified but are determined so that
%               p11+p12+p13 = p21+p22+p23 = p31+p32+p33 = 1 
%
% FROM OPEN
p11 = 0.95;            % Dwell Prob State 1
p12 = 0.05;            % Transition prob 1 --> 2
% p11 + p12 MUST BE <= 1
%
% FROM CLOSED
p21 = .1;            % Transition prob 2 --> 1
p22 = .85;            % Dwell Prob State 2
% p21 + p22 MUST BE <= 1
%
% FROM STUCK
p31 = 0;            % Transition prob 3 --> 1
p32 = 0.003;         % Transition prob 3 --> 2
% p31 + p32 MUST BE <= 1

% Expected Current Measured in each state + NoiseLevel
IOpen = 1;           % Current with Channel is Open
IClosed = 0;         % Current with Channel Closed
IStuck = 0;          % Current with Channel Stuck/Inactive
INoiseLevel = 0.01;  % Measurement Noise Level

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Create needed Matricies
statePMF = [Ostart Cstart (1-Ostart-Cstart)];
trans_matrix = [p11, p12, (1 - p11 - p12);
                p21, p22, (1 - p21 - p22);
                p31, p32, (1 - p31 - p32)];
OutputCurrents = [IOpen IClosed IStuck];

% Initialize Log-Likelihood
loglike = 0;

% print initial PMF
disp('Initial PMF');
printPMF(statePMF);

k = 0;
while(1)
  k = k + 1; 
  disp(['%%%%%%%%%%%% Data Point ',num2str(k),' %%%%%%%%%%%%']);
  disp('-------------Step 1-------------');
  % Prediction
  statePMF = statePMF*trans_matrix;
  disp(['Prediction PMF for step ',num2str(k),...
      ' Conditioned on Previous ', num2str(k-1),' Measurements']);
  printPMF(statePMF);

  % Update
  disp('-------------Step 2-------------');
  data = input('Tell me the data point: ');
  Gaussian = exp(-((data - OutputCurrents)/(2*INoiseLevel)).^2) ...
      / (INoiseLevel * 2 * pi);
  disp('Probability Density of Data given State (Evaluate Bell Curve Funct)');
  printCOND(Gaussian);

  unnormalizedPMF = Gaussian .* statePMF;
  if any(unnormalizedPMF < 0)
    keyboard
  end
  disp('-------------Step 3-------------');
  disp('Conditional Probability Weighted by the Prior (Prediction) PMF')
  printWEIG(unnormalizedPMF)

  marginallike = sum(unnormalizedPMF);
  disp('-------------Step 4-------------');
  disp('The Marginal Likelihood is the sum of the Weighted Conditional Probs');
  disp(['Marginal likelihood ',num2str(marginallike)]);
  % Log-Likelihood 
  loglike = loglike + log(marginallike);
  disp(['Running sum of logs of marginal likelihoods ',num2str(loglike)]);

  statePMF = unnormalizedPMF / marginallike; 
  disp('-------------Step 5-------------');
  disp('The updated state PMF is the Weighted Conditional Probs normalized');
  disp('   by the marginal measurement likelihood (Bayes Rule)');
  disp(['Posterior PMF conditioned on all prior measurements']);
  printPMF(statePMF);

end

%%%%%%%%%%%%%%%%%%%%
% FUNCTIONS NEEDED %
%%%%%%%%%%%%%%%%%%%%

function printPMF(pmf)
% PRINT PMF
disp(['P(Open  ) = ',num2str(pmf(1))]);
disp(['P(Closed) = ',num2str(pmf(2))]);
disp(['P(Stuck ) = ',num2str(pmf(3))]);

function printCOND(c)
% PRINT PMF
disp(['P(Data|Open  ) = ',num2str(c(1))]);
disp(['P(Data|Closed) = ',num2str(c(2))]);
disp(['P(Data|Stuck ) = ',num2str(c(3))]);

function printWEIG(c)
% PRINT PMF
disp(['P(Data|Open  )*P(Open  ) = ',num2str(c(1))]);
disp(['P(Data|Closed)*P(Closed) = ',num2str(c(2))]);
disp(['P(Data|Stuck )*P(Stuck ) = ',num2str(c(3))]);
