package scientific.stats.continuous

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

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

/*
 * Probability density function for cauchy distribution
 */
public func cauchyPDF(x: Float64, alpha: Float64, beta: Float64): Float64 {
    return beta / (Float64.getPI() * (pow(beta, Float64(2)) + pow(x - alpha, Float64(2))))
}

/*
 * Cumulative density function for cauchy distribution
 */
public func cauchyCDF(x: Float64, alpha: Float64, beta: Float64): Float64 {
    return 0.5 + 1.0 / Float64.getPI() * atan(x - alpha / beta)
}

/*
 * Sampling from the cauchy distribution
 */
public func cauchySample(m: Random, alpha: Float64, beta: Float64): Float64 {
    var u: Float64 = m.nextFloat64()
    var x = (alpha - beta) / tan(Float64.getPI() * u)
    return x
}

@Test
public class TestCauchy {
    @TestCase
    func testCauchyPDF(): Unit {
        @Assert(approxEqual(cauchyPDF(-2.0, 0.0, 1.0), 0.06366197723675814, atol:1e-13))
        @Assert(approxEqual(cauchyPDF(-1.0, 0.0, 1.0), 0.15915494309189535, atol:1e-13))
        @Assert(approxEqual(cauchyPDF(0.0, 0.0, 1.0),  0.3183098861837907,  atol:1e-13))
        @Assert(approxEqual(cauchyPDF(1.0, 0.0, 1.0),  0.15915494309189535, atol:1e-13))
        @Assert(approxEqual(cauchyPDF(2.0, 0.0, 1.0),  0.06366197723675814, atol:1e-13))
    }

    @TestCase
    func testCauchyCDF(): Unit {
        @Assert(approxEqual(cauchyCDF(-2.0, 0.0, 1.0), 0.14758361765043326, atol:1e-13))
        @Assert(approxEqual(cauchyCDF(-1.0, 0.0, 1.0), 0.25,                atol:1e-13))
        @Assert(approxEqual(cauchyCDF(0.0, 0.0, 1.0),  0.5,                 atol:1e-13))
        @Assert(approxEqual(cauchyCDF(1.0, 0.0, 1.0),  0.75,                atol:1e-13))
        @Assert(approxEqual(cauchyCDF(2.0, 0.0, 1.0),  0.8524163823495667,  atol:1e-13))
    }
}