%Programs to learn about stochastic (random) environments
%This one is stochastic environments 1, uniform probability (sev1uni.m)
%

%Item 1. Simplest case. Environment is a one-dimensional thing, say
%the amount of total food available to an individual in a year. Time unit
% is years, so time t = 0, 1, 2 etc means years 0, 1, 2, etc
% The environment is a variable F(t) indicating food supply in each year.
% Every year is unpredictable in the same way, food ranges from a miminum
% level f1 to a maximum level f2, and there is a probability distribution
% of the food level, call it g(f).

% In year t, the food level F(t) is randomly chosen according to the
%distribution g(f). 
%Empirically, you would measure F(t) over several years
%and make a histogram -- this would estimate g(f).
%Theoretically, you would use a model distribution, and we consider two
%examples below.

%Example 1.1 The distribution is uniform, meaning that g(f) is constant
%(the same)for every f, food level, between f1, minimum food level, 
% and f2, the maximum food level. 
% Let us suppose that average food level (in some arbitrary relative units) is



fbar = 1;



%and that the range is given by a number D, so that
%f1 = fbar - D, and f2 = fbar + D. The reason for using D here
%is simply that we can change its values -- small D means small variation
% large D means large variation.

%So start with small D = 0.2

D = 0.2;

f1 = fbar - D;  f2 = fbar + D;

%so food level ranges from 0.8 to 1.2, with a mean of 1.

%Because g(f) is constant and Integral {g(f) df} (total probability density)
% must be 1, we have

%g(f) = (1/(2*D)), in the range f1 <= f <= f2, 
% At each food level is, g = 1/2*.2 = 1/.4 = 2.5
% 2.5*.4 = 1, gives the area under the g-curve from 0.8 to 1.2.
% and g(f) = 0 outside this range.

%To represent this case in Matlab we will make a vector of 100 equally spaced
%values between 0.5*f1 & 1.5*f2 and find g at each point.
%the code below illustrates a useful Matlab device, linspace
% for more information on it, type 'help linspace'

fvector = linspace(0.8*f1,1.2*f2,100); %generates 100 equally spaced points on the interval
inrange = ((fvector >= f1) & (fvector <= f2));
outofrange = ((fvector < f1) | (fvector > f2));

figure(1)  %see how in range is inverse of outofrange
subplot(3,1,1)
plot(inrange)

subplot(3,1,2)
plot(outofrange)

gvector(inrange) = (1/2/D); %assign probabilities accordingly
gvector(outofrange) = 0;
subplot(3,1,3)
plot(gvector)

figure(2)
plot(fvector, gvector, '*')
xlabel('f = Food supply');
ylabel('g(f) = Probability Density');

%So much for the shape of the probability distribution

%

% To simulate stochastically means that we must generate a sample of values
%  of food level that MIGHT occur in some period of years,
%  historically called Monte Carlo simulation (casinos-gambling-etc).

%

%Set a period of T years (t = 1, 2, ..., T)
%Matlab can generate random values from many distributions

%The basic Matlab function
% X = rand(1,1) produces a vector with 1 row & 1 column
% where X is drawn according to a uniform distribution
% with values between 0 and 1

% Type 'help rand' and read up on this.


% We want T values chosen from a uniform distribution
% on [f1, f2]. So we should set
% F = f1 + (f2 - f1)*X 
%    [ why?]

%

%Let us take

T = 100;

X = rand (1,T); %meaning that we get a vector of T values

Fvector = f1 + (f2 - f1)*X; %which are the randomly chosen food levels



%Explore these values
figure(3)
plot(1:T, Fvector) %shows the values you got
xlabel('time')
ylabel('food supply')

figure(4)
hist(Fvector) %is a histogram of these 100 values -- does it look like g(f)?
title('T=100')
%

%Find moments -- first moment is mean
meanf = sum(Fvector)/T; %this is the 'longhand' way
meanf = mean(Fvector) %uses a Matlab function

%Is meanf = fbar?

%second moment is expectation(F^2) -- but we are more interested in
%variance(F) = E(F^2) - {E(F)}^2

% We can get the Standard Deviation = (Variance)^0.5 by
% another Matlab function

SDevF = std(Fvector)


%

%The differences between g & hist, meanf & fbar, etc are caused by
%the fact that we are looking a SAMPLE of size T. As T gets LARGE we
%get closer to the true values (with an error like 1/sqrt(T))

%

%TRY T = 1000 and SEE WHAT YOU GET

T = 1000;
X = rand (1,T); %meaning that we get a vector of T values
Fvector = f1 + (f2 - f1)*X; %which are the food levels

figure(5)
hist(Fvector)
title('T=1000')

meanf1000 = mean(Fvector) %uses a Matlab function

%Is meanf = fbar?

%second moment is expectation(F^2) -- but we are more interested in
%variance(F) = E(F^2) - {E(F)}^2

% We can get the Standard Deviation = (Variance)^0.5 by
% another Matlab function

SDevF1000 = std(Fvector)


%

%In this model the value of F(t) has NO effect on the value of F(t+1)
%How do we test this in a sample?
%One way is a Scatterplot

figure(5)
plot(Fvector(1:(T-1)), Fvector(2:T), '.')

%Another way is to use Matlab function 'corrcoef' or 'cov' -- TRY IT

corrcoef(Fvector(1:(T-1)), Fvector(2:T))
cov(Fvector(1:(T-1)), Fvector(2:T))


%2. Using a uniform distribution to choose students at random each morning.
%  We will generate generate the code for sampling students at random
%  there are k students... we can divide the space between 0 and 1 into equal intervals
k=12
studentspaces = linspace(0,1,k+1)

X = rand(1,10)

%who should go?