package scientific.stats.continuous

import std.math.*
import std.unittest.*
import std.unittest.testmacro.*

import scientific.numbers.*
import scientific.stats.utils.*
import scientific.stats.random.*

/*
 * Probability density function for erlang distribution
 */
public func erlangPDF(x: Float64, m: Int64, beta!:Float64 = 1.0): Float64 {
    if (x < 0.0) {
		return 0.0
	} else {
	    return pow(beta, Float64(-m)) * pow(x, Float64(m - 1)) / Float64(fac(m - 1)) * pow(Float64.getE(), - x / beta)
	}
}

public func erlangSample(r: Random, m: Int64, beta!:Float64 = 1.0): Float64{
    var u: Float64 = 1.0
    var x: Float64 = 0.0
    for (i in 0..UInt32(m)) {
		u *= r.nextFloat64()
	}
    x = -beta * log(u)
    return x
}

@Test
public class TestErlang {
    @TestCase
    func testErlang(): Unit {
        @Assert(approxEqual(erlangPDF(0.0, 1), 1.0,      atol:1e-6))
        @Assert(approxEqual(erlangPDF(1.0, 1), 0.367879, atol:1e-6))
        @Assert(approxEqual(erlangPDF(2.0, 1), 0.135335, atol:1e-6))
        @Assert(approxEqual(erlangPDF(3.0, 1), 0.049787, atol:1e-6))
        @Assert(approxEqual(erlangPDF(4.0, 1), 0.018316, atol:1e-6))
    }
}