% Summary: Beta distribution skews to low numbers
% (betadistribution.eps), so it spreads out
% the number of prey each predator has, so that even those very high
% in trophic level can also have just a few prey (see
% comparing_no_prey_orig_and_generalised_cascade.eps)
% The purpose of this script is to plot the distributions of ri for
% the Cascade Model and the Generalised Cascade Model (from Stouffer)
% to try and figure out what the significance of the exponential
% function is
%
% In the original Cascade Model, species j with nj < ni becomes a prey
% of i with fixed probability x0 = 2CS/(S - 1)
%
% In the Generalised Cascade Model, species i preys on
% species j with nj < ni, with a species-specific probability
% x drawn—from the beta distribution or an exponential
% distribution—from the interval [0, 1]
%
% So what we should see a difference in is the number of prey versus
% ni
%
% The Beta Distribution describes a distribution of x, 0 < x < 1,
% where the probability of a particular value for x is given by:
%
% p(x) = \beta(1-x)^{\beta-1}
%
% If y is a uniform random variable, you can generate x with a Beta
% Distribution using:
%
% x = 1-(1-y)^(1/beta)
%
% The Beta Distribution has the expected value E(x) = 1/(1+beta), so
% Beta can also be selected to satisfy the connectance criteria c, so
%
% beta = 1/(2*c) - 1;
%
% Let's use Little Rock Lake, which has 92 species and a connectance of 0.12
S = 92;
C = 0.12;
beta = 1/(2*C) - 1;
n = rand(S,1); % random niche value
% Niche Model
y = rand(S,1);
x = (1-(1-y).^(1/beta));
r = n.*x;
c = rand*(n-r/2)+r/2;
niche_no_prey = zeros(S,1);
for i = 1:S
ci = c(i);
ri = r(i);
low_bound = ci-ri/2;
upp_bound = ci+ri/2;
niche_no_prey(i) = length(find((n > low_bound) & (n < upp_bound)));
end
% Cascade Models:
potl_prey_V = zeros(S,1);
for i = 1:S
potl_prey_V(i) = length(find(n cascade_x
cascade_no_prey(i) += 1;
end
end
end
% In the Generalised Cascade Model, species i preys on
% species j with nj < ni, with a species-specific probability
% x drawn—from the beta distribution or an exponential
% distribution—from the interval [0, 1]
gen_cascade_no_prey = zeros(S,1);
for i = 1:S
potl_prey = find(n predator_specific_x
gen_cascade_no_prey(i) += 1;
end
end
end
hold off
plot(n,potl_prey_V,';maximum possible;ko')
hold on
plot(n,cascade_no_prey,';original cascade;bx')
plot(n,gen_cascade_no_prey,';generalised cascade;rx')
%plot(n,niche_no_prey,';niche;kx')
legend ({"maximum possible", "original cascade","generalised cascade"}, "location", "northwest")
hold off