% PROGRAM harvestopt.m
% This MATLAB code simulates a simple population growth situation
% in which a population increases by one individual each time period
% with probability p and decreases by one individual each time
% period with probability 1-p. The population is harvested when it
% reaches level b, and a fraction f of the population is harvested
% at that time, so that the population level then drops to the
% greatest integer less than b*(1-f). This code allows the user to
% investigate how large b and f can be made, for particular values of p,
% while maintaining the population above a certain minimum acceptable
% level with high probability (this is the minimum viable population
% size).
% Thus, this is the level-crossing
% calculation in which one wants to choose the level b at which to
% allow harvesting to occur so as to maximize the harvest while
% ensuring with high probability that the population level does not
% fall too low. The population is assumed to start with 100 individuals.
% The minimum viable population size is assumed to be .1 of this
% that is 10 individuals.
% The user is prompted to supply
% the probability p, the number n of simulations for each b value
% used, the farction f of the population harvested 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 number of harvests
% and a vector y giving the total amount of harvest
% for each simulation. The total simulation time is 300 time units.
% The mean and standard deviation of each of these vectors is saved for
% each b value in the vectors mtim,stim for x, and mtotharv, stotharv for y.
% Graphs of these means and std are plotted as a function of the
% b values.
%
p=input('p = probability of going up one unit: ')
f=input('f = fraction of population harvested when level b is reached: ')
n=input('n= number of simulations for each b:' )
blow=input('blow=lower value for range of b at which to harvest: ')
bhi=input('bhi=upper value for range of b at which to harvest: ')
x=[ ];
y=[ ];
c=[ ];
mtotharv=[ ];
stotharv=[ ];
mtim=[ ];
stim=[ ];
popinit=100;
minviable=.1*popinit;
% 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 number of harvests is
% ttot and total harvest is totharv
for j=1:n
ttot=0;
totharv=0;
itime=0;
xn=popinit;
while (itime < 300);
if (xn < minviable)
ttot=0;totharv=0;
break;
else
end;
if (xn < bb)
a=rand(1,1);
if a <= p
xn=xn+1;
else
xn=xn-1;
end;
itime=itime+1;
else
ttot=ttot+1;
totharv=totharv+bb*f;
xn=bb*(1-f);
end;
end;
% The vector x contains the total number of harvests
% the vector y contains the total harvest for each simulation
x=[x ttot];
y=[y totharv];
end;
mtim=[mtim mean(x)];
stim=[stim std(x)];
mtotharv=[mtotharv mean(y)];
stotharv=[stotharv std(y)];
x=[ ];
y=[ ];
end;
plot(bval,mtim,'o',bval,stim,'+')
title('mean number of harvests for various b values as o, std as +'),pause
plot(bval,mtotharv,'o',bval,stotharv,'+')
title('mean total harvest for various b values as o, std as +'),pause
disp('Values of harvest level b used')
bval,pause
disp('Mean number of harvests for each harvest level')
mtim,pause
disp('Mean total harvest for each harvest level')
mtotharv,pause