package scientific.stats.continuous

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

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

/*
 * Log of Probability density function
 */
public func moyalLogPDF(x: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale   

    let res = - 0.5 * (y + exp(-y)) - 0.5 * log(2.0 * Float64.getPI())
    return res - log(scale)
}

/*
 * Probability density function
 */
public func moyalPDF(x: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale

    let temp = moyalLogPDF(x, loc: loc, scale: scale)
    return exp(temp)
}


/*
 * Cumulative probability density function
 */
public func moyalCDF(x: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale

    let res = unsafe{ erfc(exp(- 0.5 * y) / sqrt(2.0)) }
    return res
}


/*
 * Cumulative probability density function
 */
public func moyalLogCDF(x: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale

    let temp = moyalCDF(x, loc: loc, scale: scale)
    if (temp < 0.000001) {
        throw IllegalArgumentException("moyalLogCDF: return-value too small.")
    }
    return log(temp)
}


/*
 * PPF
 */
public func moyalPPF(q: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    if (q <= 0.0 || q >= 1.0) {
        throw IllegalArgumentException("moyalPPF: quantile out of bound.")
    }

    let temp = - 1.0 / sqrt(2.0) * normalPPF(0.5 * q)
    let res = -log(2.0 * temp * temp)
    return res * scale + loc
}

@Test
public class TestMoyal {
    @TestCase
    func testMoyal(): Unit {
        @Assert(approxEqual(moyalLogPDF(3.0, loc: 2.0, scale: 1.0), -1.6028782537903938, atol:1e-13))
        @Assert(approxEqual(moyalLogCDF(3.0, loc: 2.0, scale: 1.0), -0.6085074933067561, atol:1e-13))
        @Assert(approxEqual(moyalPPF(0.7, loc: 2.0, scale: 1.0),     3.9073598213026015, atol:1e-8))
    }
}