%DISCRETE DISTRIBUTIONS'
%in ecology, very often, we make discrete meaurements
%say, food level is 'low', 'medium' or 'high'
%call these values f(1), f(2), f(3)-- a vector of ENVIRONMENTAL STATES


%Model 1. In each period there is the same independent probability
%of getting low, medium or high food level
%these are elements of a vector of probabilities
%p = (p(1), p(2), p(3)) -- the elements of p must add to 1 -- why?


%Example

f = [0.5 1 1.5]; % low, med, hi
p = [0.4 0.2 0.4]; %note that here lo & hi are more probable than med
%we'll contrast another case, this is NOT uniform or normal

p1 = [0.2 0.6 0.2];

%Given an f and a p we want to
%SIMULATE T values of F(t) by choosing for eact t=1, 2, ..., T
% a random number distributed according to p (or p1, as the case may be)

%How to do this?
% We would like to choose State 1 with probability p(1), State 2
% with prob. p(2), etc

%DRAWS FROM DISCRETE DISTRIBUTION

%Step 1 -- Make the cumulative probability vector

psum = cumsum(p) %look up help cumsum

%Note that Prob{that a state is 1 or 2} will be psum(2), whilst
%Prob(that a state is 1 or 2 or 3} will be psum(3).


%Step 2 -- Choose a random number that is Uniformly distributed
% between 0 and 1

X = rand(1,1);


%Step 3 -- Find the State. 
% The Logic: we know that Prob{X <= w} = w for any number w between 0 & 1.  
% Thus 
% Prob{that X LIES between psum(2)and psum(3)} = (psum(3) - psum(2) = p(3).

%The Algorithm


jj = min(find(psum > X))


%This selects the environmental state j
% two steps here: first find, then min
% type 'help find'
%  i=find(a < 100), returns the indices of A where there is a number < 100
%  so... 
%  X is your random number and psum is a vector like [0.3 0.5 1]
%   if the random number is , for example, 0.7
%               then "find (psum > 0.7)" would give "3"
%  NOW the second step is
%  min[3], which clearly is 3

%now try the psum of the above example, where 
psum= [0.4 0.6 1]
find(psum>.2)
%then take the min of that to find out which interval it is
min(ans)

%try another
find(psum>.5)
min(ans)

%and then
find(psum>.7)
min(ans)


%TRY THIS -- for T = 10 periods, find the environmental state according
%to distribution p -- as follows


T = 10;

EnvState = zeros(1,10); 
%sets up the vector to receive the results of the loop
%we will put the states into this vector when we find them

for i = 1:T
   X = rand(1,1)
	j = min(find(psum > X))
   EnvState(i) = j
end

clear EnvState

%Now repeat it for T=100, but add semicolons to the code to mask the action of the
% loop

T = 100;

EnvState = zeros(1,100); 
%we will put the states into this vector when w find them

for i = 1:T
   X = rand(1,1);	
   j = min(find(psum > X));
   EnvState(i) = j;
end

%Explore the result
EnvState

figure(1)
plot(1:T, EnvState, '*-') %See how the State jumps around over time
ylabel('EnvState')
xlabel('time')

figure(2)
hist(EnvState)  % Compare this to p



%Try it again with p1 in place of p
%just go into the script and alter it and run it again...


%What about the actual food level?

for i = 1:T
   FoodLevel(i) = f(EnvState(i));%this uses the state to index the food
end


figure(3)
plot(1:T, FoodLevel,'o-')
ylabel('FoodLevel')
xlabel('Time')

figure(4)
hist(FoodLevel)
ylabel('Frequency')
xlabel('FoodLevel')

%Moments -- We are usually interested in Food level
%Theoretically,

fbar = sum(p.*f) %this is the average sum(p(i) f(i))
f2bar = sum(p.*f.*f) %this is average of F^2
Varf = f2bar - fbar; %the variance
CVf = sqrt(Varf/fbar/fbar) %this is the Coefficient of Variation

%Compute these for p & p1; compare

%How would you compute these from the simulated values?
%Example --

meanf = sum(FoodLevel)/T;

%END