% 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