% PROGRAM levelopt.m
% This MATLAB code simulates a random walk with probability p
% of going up one unit, r of remaining at the same state, and q
% of going down one unit. It assigns a cost of being in each
% state according to a function c state^2, and assumes there is a
% reinitialization cost K for moving the system back to state 0
% once the state hits a distance b from the origin.
% Thus, this is the level-crossing
% calculation in which one wants to choose the level b at which to
% reinitialize the system so as to minimize the long-term cost of
% (K+ A)/B
% where A is the sum of costs in the various states until level b
% is reached, and B is the time to hit level b.
% The user is prompted to supply
% the probabilities p, and q (if p+q <1 then the walk stays in the
% current state with a positive probability), the cost K, the number n of
% simulations for each b value used, and the range of b-values
% to use (the code partitions this range into mm equal pieces).
% At the end, output will be a vector, x, giving the times til
% the level b is reached, a vector y giving the costs along each
% sample path associated with that b value, and rat which contains the
% cost/time ratio along each path. Each of these is of length n.
% The mean and standard deviation of each of these vectors is saved for
% each b value in the vectors mtim,stim for x, mcost, scost for y, and
% mrat,srat for rat. Finally, the long-term cost calculated as the
% ratio of the means as indicated above is calculated.
% Graphs of these means and std are plotted as a function of the
% b values and a graph showing the relationship between the two methods
% of calculating long-term cost/time is shown.
%
!c:
p=input('p = probability of going up one unit: ')
q=input('q = probability of going down one unit: ')
K=input('K = reinitialization cost: ')
c=input('c = constant in cost function: ')
n=input('n= number of simulations for each b:' )
blow=input('blow=lower value for range of b to reset system to 0: ')
bhi=input('bhi=upper value for range of b to reset system to 0: ')
% nh will count the number of moves up, nt the number of moves down
x=[ ];
y=[ ];
rat=[ ];
mrat=[ ];
srat=[ ];
mcost=[ ];
scost=[ ];
mtim=[ ];
stim=[ ];
% mm is the number of cases for b to run
mm=10;
bval(1)=blow;
for i=2:(mm+1)
bval(i)=blow+(i-1)*(bhi-blow)/mm;
end;
% iterate over the mm+1 b values
bval
for ii=1:(mm+1)
bb=bval(ii);
% for each b value, make n simulations
% here the current state is xn, total time is
% ttot and total cost is totcost
for j=1:n
ttot=0;
totcost=0;
xn=0;
z=[ ];
i=0;
while (xn < bb & xn > -bb)
i=i+1;
a=rand(1,1);
if a <= p
xn=xn+1;
elseif a >= 1-q
xn=xn-1;
else
xn=xn;
end;
ttot=ttot+1;
totcost=totcost+c*xn*xn;
z=[z xn];
end;
% The vector x contains the times til reaching level b
% the vector y contains the total costs for each simulation
% the vector z contains the states of the system
% the vector rat contains the sample path cost/unit time
rr=(K+totcost)/ttot;
x=[x ttot];
y=[y totcost];
rat=[rat rr];
%plot(z,'-')
%title('plot of one sample path'),pause
end;
mtim=[mtim mean(x)];
stim=[stim std(x)];
mcost=[mcost mean(y)];
scost=[scost std(y)];
mrat=[mrat mean(rat)];
srat=[srat std(rat)];
x=[ ];
y=[ ];
rat=[ ];
end;
rate=(K+mcost)./mtim;
plot(bval,mtim,'o',bval,stim,'+')
title('mean time for various b values as o, std as +'),pause
plot(bval,mcost,'o',bval,scost,'+')
title('mean cost for various b values as o, std as +'),pause
plot(bval,mrat,'o',bval,srat,'+')
title('sample path cost/time ratio as o, std as +'),pause
rate=(K+mcost)./mtim;
plot(bval,rate,'o',bval,mrat,'+')
title('ratio of means of cost/time as o, sample path costs/time as +'),pause
plot(rate,mrat,'o')
title('ratio of means of cost/time on x-axis, sample path costs/time on y-axis'),pause