% SimHMC.m
%
% Neal Patwari, 4 Nov 2010
%
function [S,X] = SimHMC(pr_s_given_x, P, p_0, duration)

% Assume, for convenience, that X_n are in the set 1, 2, ..., possible_x
% and that S_n are in the set 1, 2, ..., possible_s.
[possible_x, possible_s] = size(pr_s_given_x);

% cumsum is an easy way to get a CDF.
S_CDF   = cumsum(pr_s_given_x,2);
X_0_CDF = cumsum(p_0);
P_CDF   = cumsum(P,2);

% The expression CDF < rand will be one until the random value falls below
% a CDF value.  The expression  "sum(CDF < rand) + 1" finds the index of
% that first CDF value greater than the random number.  This is an easy way
% to generate random variables with an arbitrary pmf.
X(1)    = sum(X_0_CDF < rand) + 1;
S(1)    = sum(S_CDF(X(1),:) < rand) + 1;
for n = 2:duration,
    X(n)    = sum(P_CDF(X(n-1),:) < rand) + 1;
    S(n)    = sum(S_CDF(X(n),:) < rand) + 1;
end
