JuliaStats/Distributions.jl — pull request #254 by bicycle1885 (2014-06-28 16:13:15+00:00)
Title: WIP: more accurate KSDist function

This is a work-in-progress implementation of Simard and L'Ecuyer (2011) meta-algorithm .
The current commit contains debug code and duplicated functions in order to check the differences between the current and the proposed implementation.

I've used KolmogorovSmirnovDist.c as a reference, which is an implementation of the meta-algorithm.
At some inputs, there seems to be serious dispersion between the current one and the reference.
For example when n = 5 and x = 2.256281e-01 the CDFs are:

reference | current | proposed
-------------|-----------|--------------
8.852118e-02 | 9.711784e-02 | 9.711784e-02

In this case, the `cdf_durbin` function is in charge, which is reused in the proposed implementation.

For other example when n = 1000000 and x = -3.366834e-02 the CDFs are completely different:

reference | current | proposed
-------------|----------|-------------
1.000000e+00 | 0.000000e+00 | 1.000000e+00

What I'm wondering are:

* in the first example, is the current implementation wrong?
* in the second example, returning zero looks to be sensible, am I right?

Here is the complete log of comparison: [out.txt](https://dl.dropboxusercontent.com/u/25281592/out.txt)
(**CAUTION - it's 696KB**)

And here is the way for quick checking and reproducing it:

```
# check out my branch
curl -O -O http://simul.iro.umontreal.ca/ksdir/KolmogorovSmirnovDist.{c,h}
clang -fPIC -shared -dynamiclib KolmogorovSmirnovDist.c -o libksd.dylib
julia check.jl > out.txt
```

check.jl:
```julia
using Distributions

function _cdf(d, x)
ccall((:KScdf, "libksd"), Float64, (Int, Float64), d.n, x)
end

begin
ns = vcat([1:50], [100, 200, 500, 1000, 2000, 10_000, 50_000, 100_000, 500_000, 1_000_000])
xs = linspace(-0.1, 1.1, 200)
n_all = 0
n_eq1 = 0
n_eq2 = 0
println("n\tx\treference\tcurrent\tproposed\tcurrent_ok\tproposed_ok")
for n in ns
ksd = Distributions.KSDist(n)
for x in xs
v = _cdf(ksd, x)
v1 = Distributions.cdf(ksd, x)
v2 = Distributions.cdf2(ksd, x)
eq1 = isapprox(v, v1)
eq2 = isapprox(v, v2)
@printf "%d\t%12.6e\t%12.6e\t%12.6e\t%12.6e\t%s\t%s\n" n x v v1 v2 (eq1 ? "o" : "x") (eq2 ? "o" : "x")

n_all += 1
if eq1
n_eq1 += 1
end
if eq2
n_eq2 += 1
end
end
end

println(STDERR, "1: $n_eq1/$n_all")
println(STDERR, "2: $n_eq2/$n_all")
end
```

Thank you.