package scientific.stats.continuous

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

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


/*
 * Log of Probability density function
 */
public func alphaLogPDF(x: Float64, a: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

    let y = (x - loc) / scale 

    if(y <= 0.0 || a <= 0.0) {
        throw IllegalArgumentException("alphaLogPDF: parameters out of bound.")
    }

    if(scale == 0.0) {
        throw IllegalArgumentException("alphaPDF: scale cannot be 0.")
    }

    let t1 = normalCDF(a)

    return - 2.0 * log(y) - log(t1) - 0.5 * log(2.0 * Float64.getPI()) - 0.5 * (a - 1.0 / y) * (a - 1.0 / y) - log(scale)
}


/*
 * Probability density function
 */
public func alphaPDF(x: Float64, a: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

    let y = (x - loc) / scale 

    if(x <= 0.0 || a <= 0.0) {
        throw IllegalArgumentException("alphaPDF: parameters out of bound.")
    }

    if(scale == 0.0) {
        throw IllegalArgumentException("alphaPDF: scale cannot be 0.")
    }

    return exp(alphaLogPDF(x, a, loc: loc, scale: scale))
}



/*
 * Cumulative probability density function
 */
public func alphaCDF(x: Float64, a: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

    let y = (x - loc) / scale 

    if(y <= 0.0 || a <= 0.0) {
        throw IllegalArgumentException("alphaPDF: parameters out of bound.")
    }

    return normalCDF(a - 1.0 / y) / normalCDF(a) / scale
}


/*
 * Log of cumulative probability density function
 */
public func alphaLogCDF(x: Float64, a: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

    let y = (x - loc) / scale 

    if(y <= 0.0 || a <= 0.0) {
        throw IllegalArgumentException("alphaPDF: parameters out of bound.")
    }

    let temp = alphaCDF(x, a, loc: loc, scale: scale)

    if(temp < 0.000001) {
        throw IllegalArgumentException("alphaPDF: return-value too small.")
    }

    return log(temp)
}

@Test
public class TestAlpha {
    @TestCase
    func testAlpha(): Unit {
        @Assert(approxEqual(alphaLogPDF(1.0, 2.0), -1.3959256238757092,  atol:1e-13))
        @Assert(approxEqual(alphaPDF(1.0, 2.0),     0.24760374232796756, atol:1e-13))
        @Assert(approxEqual(alphaCDF(1.0, 2.0),     0.8609310408460744,  atol:1e-13))
    }
}