% compute -log(P(K|Z)) for K=0..K_max
%
% using x23.m to generate the data
%
% define L and K_max before running this script.
%
% Due to the different pseudo-random number generators used in Matlab and
% Octave, the initial pseudo-random values (noise, missed detections, and
% false alarms) will differ between the two simulation environments.  In
% either case, the startup.m and .octaverc startup scripts here set things so
% that the pseudo-random number generators are NOT initialized randomly based
% on the time of day.  So one should be able to duplicate the "random" data
% used in the examples here in multiple separate simulation sessions.

% Wmin will be used as a scale factor to prevent underflow in P
%
% Note that for unreachable final states, exp(-HUGE_VAL) == 0
% so these states do not affect the probability calculations here.

% P will include all final cost results vs. K
% P1 will include just the single lowest costs vs. K
%
P  = zeros( K_max, 1);
P1 = zeros( K_max, 1);

for K = 1:K_max
  x23			% setup data
  vit			% run Viterbi algorithm
  W = sort(w(:));
  if( K == 1)		% set scale factor
    Wmin = W(1);
  end
  P(K) = sum( exp( -(W - Wmin)));
  P1(K) = exp( -(W(1) - Wmin));
end


% compute P(K=0|Z) manually, Viterbi not needed
%
if( N > 1)
  W0 = missed_cost * sum(MM(2:N));
else
  W0 = 0
end

for i=1:N
  W0 = W0 - log( p_poisson( MM(i), false_alarm_rate) );
end

P0=exp(-(W0-Wmin));


% Compute sum over tau' and normalize.
%
% The normalization is done only for convenience and does not affect the
% analysis of which P(K|Z) is maximum.
%
PKZ = [ P0; P];

PKZ = -log(PKZ/sum(PKZ));

disp( 'K    PKZ');

[ [0:K_max]' PKZ]


% if L>1, compute PKZ1 using just the single lowest costs vs. K
%
if( L > 1)

  PKZ1 = [P0; P1];

  PKZ1 = -log(PKZ1/sum(PKZ1));

  disp( 'K    PKZ1');

  [ [0:K_max]' PKZ1]
end
