ewens_sampler
=============

A simple perl script to generate the expected site frequency spectrum for a
sample of size n under a population mutation rate theta, following the Ewens
sampling formula for a neutral population in demographic equilibrium under the
infinite alleles mutation model.


USAGE

    ewens.pl [options] theta n

    theta   population mutation rate (4*N*mu)
    n       sample size

    options:
      -u, --unconditional   use the unconditional expected SFS (see MODES)
      --log                 log-scale the chart bars (the SFS spans several
                            orders of magnitude)
      --rows N              show at most N allele classes in the chart
                            (default 50; N < 1 shows all n)
      --no-plot             skip the chart (values still go to file + stdout)
      -h, --help            show help and exit

Run "ewens.pl" with no arguments (or "ewens.pl --help") to see the full usage.

The expected number of allele classes of size j (for j = 1 .. n) is written,
one value per line, BOTH to a file called "ewens.out" and to stdout. An ASCII
bar chart of the spectrum is drawn to stderr, so it appears on screen even when
stdout is redirected, e.g.

    ewens.pl 5 100 > ewens.out      # values to the file, chart still on screen

For large n the chart shows the first 50 allele classes; all values still go to
the file and stdout.


MODES

    (default)             The expected frequency spectrum conditional on
                          k = the rounded expected number of alleles. This is
                          computed from Stirling numbers of the first kind and
                          reproduces the output of the original script.

    -u, --unconditional   The unconditional expected site frequency spectrum,

                              E[a_j] = (theta / j)
                                       * n! / (n-j)!
                                       * Gamma(theta + n - j) / Gamma(theta + n)

                          This is the standard closed-form Ewens result. It
                          satisfies sum_j E[a_j] = E[K] and sum_j j*E[a_j] = n
                          exactly, and is what you usually want for a point
                          estimate of the expected SFS.


NOTES

  - Both modes report point expectations, not full distributions or confidence
    intervals.

  - All arithmetic is done in log space with native double-precision floats, so
    the script no longer depends on Math::BigFloat and does not overflow, even
    though the Stirling numbers and factorials involved are astronomically
    large. Results match an exact arbitrary-precision computation to about 15
    significant digits.

  - This is a rewrite of an earlier Math::BigFloat version. It is far faster
    (hundreds of times faster for n in the hundreds, and it scales to n in the
    tens of thousands in a fraction of a second, where the original was
    effectively intractable), while producing the same numbers.


REQUIREMENTS

  - perl (core modules only; Math::BigFloat is no longer required).