incomplete gamma, chi square, uniform distribution

chi2.go

package chi2

import "math"

func Chi2p(dof int, distance float64) float64 {
    return IncGammaQ(.5*float64(dof), .5*distance)
}

func IncGammaQ(a, x float64) float64 {
    if x < a+1 {
        return 1 - gser(a, x)
    }
    return gcf(a, x)
}

func gser(a, x float64) float64 {
    gln, _ := math.Lgamma(a)
    ap := a
    sum := 1 / a
    del := sum
    for n := int(20 * (1 + math.Sqrt(a))); n > 0; n-- {
        ap++
        del *= x / ap
        sum += del
        if math.Abs(del) < math.Abs(sum)*1e-14 {
            return sum * math.Exp(-x+a*math.Log(x)-gln)
        }
    }
    // I think it's a bug if it doesn't converge
    // I don't think this should be turned into an error
    panic("gser: no convergence")
}

func gcf(a, x float64) float64 {
    var gold, fac, b1 float64 = 0, 1, 1
    var b0, a0 float64 = 0, 1
    a1 := x
    gln, _ := math.Lgamma(a)
    for n := int(20 * (1 + math.Sqrt(a))); n > 0; n-- {
        an := float64(n)
        ana := an - a
        a0 = (a1 + a0*ana) * fac
        b0 = (b1 + a0*ana) * fac
        anf := an * fac
        a1 = x*a0 + anf*a1
        b1 = x*b0 + anf*b1
        if a1 != 0 {
            fac = 1 / a1
            g := b1 * fac
            if math.Abs((g-gold)/g) < 1e-14 {
                return math.Exp(-x+a*math.Log(x)-gln) * g
            }
            gold = g
        }
    }
    // I think it's a bug if it doesn't converge
    // I don't think this should be turned into an error
    panic("gcf: no convergence")
}

uniform_test.go

package chi2

import "testing"

func TestUniform(t *testing.T) {
    for _, dset := range [][]int{
        {199809, 200665, 199607, 200270, 199649},
        {522573, 244456, 139979, 71531, 21461},
    } {
        utest(t, dset)
    }
}

const sigLevel = .05

func utest(t *testing.T, dset []int) {
    var sum int
    for _, c := range dset {
        sum += c
    }
    t.Log(" dataset:", dset)
    t.Log(" samples:                      ", sum)
    t.Log(" categories:                   ", len(dset))

    dof := len(dset) - 1
    t.Log(" degrees of freedom:           ", dof)

    dist := chi2ud(dset)
    t.Log(" chi square test statistic:    ", dist)

    p := Chi2p(dof, dist)
    t.Log(" p-value of test statistic:    ", p)

    sig := p < sigLevel
    t.Logf(" significant at %2.0f%% level?      %t\n", sigLevel*100, sig)
    t.Log(" uniform?                      ", !sig, "\n")
}

func chi2ud(ds []int) float64 {
    var sum, expected float64
    for _, d := range ds {
        expected += float64(d)
    }
    expected /= float64(len(ds))
    for _, d := range ds {
        x := float64(d) - expected
        sum += x * x
    }
    return sum / expected
}

So, this isn't much. The test data is from Rosetta Code. I don't seem to have the original CF code. This is pretty much straight from NR. Needs much work before committing.