% Solve the Gambler's Problem (Example 4.3 page 84) by value
% iteration.

% Set up the p(s',r|s,a) function
ph = 0.4; % probability of heads
p = zeros(101,2,101,100);
for s=1:99
for a=0:min([s,100-s])
sp = s + a; % if heads
if(sp >= 100)
r = 1; % win
else
r = 0; % continue or loose
end
p(sp+1,r+1,s+1,a+1) = ph; % prob heads
sp = s - a; % if tails
p(sp+1,r+1,s+1,a+1) = 1-ph; % prob tails
end
end

V = zeros(101,1);
V([0,100]+1) = 0; % zero value in terminal states
PI = zeros(100,1); % actions
count = 0;
while true
delta = 0;
Vp = zeros(101,1);
for s=1:99
% Maximize over actions
va = zeros(min([s,100-s])+1,1);
for a=0:min([s,100-s])
for sp=0:100
for r=0:1
va(a+1) = va(a+1) + p(sp+1,r+1,s+1,a+1)*[r + V(sp+1)];
end
end
end
[vmax,ap] = max(va);
ind = find(abs(va-vmax)<1e-6); % A bit of numerical stuff
ap = min(ind); % A bit of numerical stuff
PI(s+1) = ap-1;
delta = max([delta,abs(vmax-V(s+1))]);
Vp(s+1) = vmax;
end
V = Vp;
count = count + 1;
fprintf('%4d, delta = %20.16f\n',count,delta);
subplot(211); plot([0:100],V);
subplot(212); bar([0:99],PI); drawnow; pause;
if(delta < 1e-8)
break;
end
end