% Tutorial 8

% Next four weeks - work period for group assignment

% Question 1:

% A certain population of moose can modeled as follows:
% Depending on the amount of rain that occurs in the spring (and therefore 
% the abundance of plant-life), the moose can either have a "good" year or 
% a "bad" year.  If there is lots of rain, the moose do well, and the growth 
% of the population is increased by a factor of random variable r_g 
% which is sampled from a normal distribution with mean u_g and variance s_g.  
% If there is not enough rainfall, the population suffers, and 10% of the
% moose are lost. The probabilities of having lots or little rain, 
% respectively, are p_g, and p_b.

% Take u_g=2, s_g=0.1, M(0)=200, p_g=1/3, p_b=2/3

%%
% Run and plot a single trial up to t=10
u_g=2;
s_g=0.1;
p_g=1/3;
p_b=2/3;
t=0:10;
M=zeros(1,11); % remember needs to go up to 11, 0, 1, ..., 10
M(1)=200;
% now iterate the recursive formula
% remember that we need to reselect our random variables at each step
% since they are sampled from a normal distribution, recall that to sample
% from a normal distribution with mean u and st.dev. s, use r=randn*s+u
for(i=1:10)
    % first pick if it is a good or bad year by selecting a uniformly
    % random value between 0 and 1
    x=rand(1,1);
    if(x<=p_g) % then it is a good year, so select r_g
        r_g=randn*s_g+u_g; %take the abs since there is a small chance that
                                %this could be negative
        M(i+1)=(1+r_g)*M(i);
    else %it is a bad year
        M(i+1)=(0.9)*M(i); 
    end
end
% now plot M(t)
plot(t,M,'LineWidth',3)
    xlabel('year')
    ylabel('population size')
    title('Sample moose population over 10 years')
    %ylim([0 500000])
  
%%

% Now run 500 trials up to t=50.  Display a histogram at t=50
% this time we need to let M be a 500 by 51 matrix.  Each of the 500 rows
% stores the results from one of the trials
M=zeros(500,51);
M(:,1)=200; %set the same initial condition for all trials
% use a double for loop to calculate the results of all trials
% the outer loop goes over the trials, while the inner loop calculates the 
% results for each individual trial using the recursive formula
for i=1:500
    for j=1:50
        %pick a uniformly random value between 0 and 1
        x=rand(1,1);
        if(x<=p_g)
            %then it is a good year
            r_g=randn*s_g+u_g;
            M(i,j+1)=(1+r_g)*M(i,j);
        else
            %it is a bad year
            M(i,j+1)=0.9*M(i,j);
        end
    end
end
%make a histogram with 10 bins at t=50
figure();
hist(M(:,51),10)
    xlabel('Population size at t=50 years')
    ylabel('Amount of trials')
    title('Histogram at t=50')
    %ylim([0 100])
    
% How do you interpret this histogram? The number of trials is on the 
% y-axis, and the population size is on the x-axis. The size of the bar 
% represents the number of trials who had a final size within that x-range.


%%
% Estimate how many years it will take the population to reach 500,000
% moose.

% loop over the 500 trials we ran.  Record the first year where the
% population is at or above 500,000 for each of the trials, and then take
% the average
% store the years in an array of length 500
index=zeros(1,500);
% loop over the trials.  The first time M(i,j)>500,000 we can break the
% inner loop and move on to the next trial
for i=1:500
    for j=1:51 
        if M(i,j)>=500000
            index(i)=j-1;    %because j-1 is the current year
            break;           %we can stop comparing for this trial
        end
    end
end
averageYears=sum(index)/length(index);
disp('The average years taken for the moose to reach 1/2 million is:');
disp(averageYears);

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Further question:

% A man is walking along a four block stretch.  If he is at corner 1,2 or
% 3, he turns left or right with equal probability.  He continues this
% until he reaches his home which is located at corner 0, or the bar which
% is located at corner 4.  If he is at home or at the bar, he stays there.

% 1. What is the transition matrix for this problem?

% 2. What are the fixed points of the model?

% 3. What is the time-dependent solution?

%. 4. What is the steady state of this model?