function [value, confInt] = PW311_Simulate(probs,payoff,discount,...
    startState,numRepl,numSteps,policy)
numStates = length(probs);
sumDiscRewards = zeros(numRepl,1);
discFactors = discount.^(0:(numSteps-1));
for k = 1:numRepl
    state = startState;
    Rewards = zeros(numSteps,1);
    for t = 1:numSteps
        if (state == 1) % state 1, no choice
            if rand > probs(state)
                state = 2;
            end
        elseif state == numStates
            % last state, get reward and back to square 1
            state = 1;
            Rewards(t) = payoff(numStates);
        else % at interior states, apply policy
            if policy(state) == 1
                % collect reward and go back to square 1
                Rewards(t) = payoff(state);
                state = 1;
            else % random transition
                if rand > probs(state)
                    state = state + 1;
                end
            end
        end % end if on states
    end % each replication
    sumDiscRewards(k) = dot(Rewards,discFactors);
end % all replications
[value,~,confInt] = normfit(sumDiscRewards);
