Brian J. McGill, Brian A. Maurer, and Michael D. Weisser. 2006. Empirical evaluation of neutral theory. Ecology 87:1411–1423.

Appendix D. A description of algorithms and usage of code for fitting analytical ZSM.

The analytical version is presented in "Neutral theory and relative species abundance in ecology" (Volkov et al. 2003. Nature 424:1035–1037.)

The authors refuse to provide source code for estimating the analytical model, just as for the Monte Carlo version.

The main trick is in evaluating their Eq. 7:

Calcluating <φ_n>

Where:

The first problem in evaluating this analytically is that the Γ() function blows up to infinity very quickly. As a result, most normal computers that naively start calculating Eq. 7, end up with ∞/∞. In short the equation blows up numerically. This is easily resolved using the gammaln function (denoted here lnΓ) defined as lnΓ(x)=ln(Γ(x)). This in turn implies Γ(x)=exp(lnΓ(x)) Now also recall that exp(a)exp(b)=exp(a+b). If we also take the constant in front of the integral in Eq. 7 to inside of the integral (so that very large numbers can cancel out very large numbers), then we get a rearranged from of Eq. 7

<φ_n>=θ∫exp(lnΓ(J+1)-lnΓ(n+1)-lnΓ(J-n+1)+lnΓ(γ)-lnΓ(J+γ)+lnΓ(n+x)-lnΓ(1+x)+lnΓ(J-n+γ-x)-lnΓ(γ-x)-xθ/γ)dx

The integrand (part inside the integral) is now calculable on most any math package that has the lnΓ function.

The next trick is performing the integral. The integrand looks something like the delta function - an extremely narrow spike with a value close to zero everywhere else. The height of the spike is also quite small (often less than 0.0001). Fig. D1 is an example for J=20000, θ=50, m=0.1, and n=500 (note the vertical axis goes up to 3.5E-5.

MATLAB's quad8 function works fine on this. Be sure not to use the inline function form but define an actual function to represent the integrand. However, quad8 is very slow and about to decommissioned. It took 45 minutes to calculate the SAD for the BCI data, given all the parameter.

Thus another numerical approximation needs to be invoked. This approximation is two part. For small n (where the integrand is veiled and only the right part of the spike appears for x≥0, MATLAB's quadl function works well and fast. For n "large", the method of steepest descent as suggested by Volkov et al works well and is fast. The code as written transitions from quadl to steepest descent when the value of the integrand (Eq. 7) at y=0 is "small". It is possible that for some parameters, the quadl function will "miss" the spike and report values for too low for some φ_n. This is usually obvious as a big "chunk" out of the SAD distribution.

The final numerical trick is to not evaluate <φ_n> for each n. As currently coded, the routine estimates for each n up to 4096. After that we estimate every 100 or so n (or more often) and interpolate inbetween. This greatly speeds up code for very large J/maxabundances with little loss in accuracy.

Alternatives

Vallade and Houchmandzadeh. 2003. Analytical solution of a nuetral model on biodiversity. Physical Review E:061902.

McKane, Alonso, and Soule. 2004. Analytical solution of Hubbell's model of local community dynamics. Theoretical Population Biology 65:67–73 (Note: this paper actually just presents a special case result of a more general paper by the same authors in 2000 in Physical Review E 62:8466)

The routines provided implement four different analytical solutions:

1. Volkov et al without any numerical tricks (slow integration) is by far the slowest but also safest.
2. Volkov et al with numerical tricks as described above is the fastest but most likely to run into numerical instability (fancy word for produce wrong results).
3. Vallade and Houcmandzadeh fairly fast without any numerical tricks.
4. McKane et al. also fairly fast without any numerical tricks.

1. Run estimate_anal_ZSM with bVolkov=1; and the others zero.
2. Run getphi with the parameter estimates produced in step #1 with bReallyAccurateAndSlow=0; and the others set to 0.
3. Compare the resulting outputs: plot(1:length(phivolkov),phivolkov,1:length(phivolkov),phiaccandslow);set(gca,'YScale','log'); is a quick and dirty way to check in Matlab assuming you saved the phi arrays for the two runs.
4. You may also want to run either the V&H (bValladeHouch=1) and/or the McKane et al (bMcKane=1) methods.

1. Take J & S directly from the abundance data.
2. Vary m.
3. Calculate θ from the summed version of Eq. 7 (S=∫…∑…dx) through iterative numerical solution.
4. Calculate <φ_n> for each n.
5. Calculate the Χ2 (chi-squared) goodness of fit using Preston log2 binning of the current SAD vs. the empirical SAD.
6. Repeat over various parameters of m to minimize the Χ2 value.

This approach has the net effect of replacing the parameter θ with the parameter S. In many ways this is a net gain, as S is much more easily estimated than θ. However, statistically, it does not matter - the analytical ZSM remains a three parameter model (despite claims to the contrary by Volkov et al).

This method breaks down (becomes extremely slow) if J is very large (say J>100,000) or m is very small (say m<0.01).In those conditions, I have provided a routine (estimate_anal_ZSM_theta) that estimates m and θ directly by optimizing fit for those two parameters. Although it does not exactly match the Volkov et al estimation method, it is much more feasible for the aforementioned cases and an entirely reasonable and accurate approach.

The code provided

[theta,m,phin,prestonpdf]=estimate_anal_ZSM(data)

estimates θ and m for the data. The data must be all positive integers. This function also returns the value of φ_n and the PDF on preston log2 bins

theta=solvetheta(data,m)

solves for θ given J, S and a guess at m as per the method described in Parameter estimation

plotexpphi

provides plots of the integrand for various values of parameters (not required but educational)

[phin,prestonpdf]=getphi(theta,m,J,maxabundance,S)

figures how high to vary n (based on rounding up log2 of maxabundance observed). Iterates over n=1..max and calls expphi. Assembles into an array dimensioned by n and a log2 preston pdf

phi=expphi(n,theta,J,m[,S])

evaluates Eq. 7 for one set of parameters. S is only used if m=1.0 which is the same as the logseries distribution

sadhist,myhistc

utility routines to calculate a preston log2 binning

[theta,m,phin,prestonpdf]=estimate_anal_ZSM_theta(data)

like estimate_anal_ZSM but estimates theta and m directly rather than the method published in Volkov etal which derives theta from S at considerable computational cost. This method works better if J is very large (>100,000?) or m is very small (<0.01)

Practicalities

1. Calling expphi is very fast.
2. Calling getphi takes about 15 seconds on a Pentium 500 machine.
3. Calling estimate_anal_zsm takes about 20 minutes on a Pentium 500 machine for the BCI data.
4. Times will go up or down mainly as the maximum abundance goes up or down.
5. If you have MATLAB 6.1 or earlier, you need to do two changes. First delete the calls to 'hidewarn' and 'showwarn' in expphi.m - you will now get a few warnings but you can ignore them. Second you will need an additional package of five files (mypsi.m,digamma.m,digammas.m, trigamma.m,trigammas.m) to compensate for the fact that the psi function was only added to matlab in version 6.5. The last four files come from StatBox. Just add these files to your working directory and rename mypsi.m to psi.m Email me if you need these files.
6. If you have J>>20,000 or m<0.01 you are probably best to use the estimate_anal_ZSM_theta routine (see previous section).

Steepest descent

1. Let y=f(x) have a peak at x=x0 and be very close to zero everywhere else (as does our integrand for n large).
2. Take a Taylor series expansion of f about x=x0 and recall that the first derivative f'(x0)=0 so:
   ∫f(z)dz≈f(x0)∫(1+(z-x0)^2/2! f''(x0)/f(x0) + … dz
   Now perform a change of variables with s=[(z-x0) *sqrt(f(x0)/f''(x0)/2)]^2.
   Recalling the Taylor series for exp(x) we now have
   ≈ sqrt(2*f(x0)^3/-f''(x0)) ∫ exp(-s^2) ds
   We can either evaluate the &int over the appropriate range using numerical methods for the error function (erf) or assume that the integral from -∞ to +∞ which evaluates to sqrt(π).

FIG. D1. An integrand example for J=20000, θ=50, m=0.1, and n=500 (note the vertical axis goes up to 3.5E-5).