package scientific.stats.continuous
import std.math.*
import std.unittest.*
import std.unittest.testmacro.*
import scientific.numbers.*
import scientific.stats.random.*
foreign func malloc(size: UIntNative): CPointer<Unit>
foreign func beta(a: Float64, b: Float64): Float64
foreign func cdfbet(which: CPointer<Int32>, p: CPointer<Float64>, q: CPointer<Float64>,
x: CPointer<Float64>, y: CPointer<Float64>, a: CPointer<Float64>,
b: CPointer<Float64>, status: CPointer<Int32>, bound: CPointer<Float64>): Float64
/*
* Log beta value
*/
public func betaLog(alpha: Float64, beta: Float64): Float64 {
var ans: Float64 = gammaLog(alpha) + gammaLog(beta) - gammaLog(alpha + beta)
return ans
}
/*
* Probability density function for beta distribution
*/
public func betaPDF(x: Float64, a: Float64, b: Float64): Float64 {
let ans: Float64 = pow(1.0 - x, b - 1.0) * pow(x, a - 1.0) / unsafe {beta(a, b)}
return ans
}
/*
* Cumulative probability density function for beta distribution
*/
// Input, int *WHICH, indicates which of the next four argument
// values is to be calculated from the others.
// 1: Calculate P and Q from X, Y, A and B;
// 2: Calculate X and Y from P, Q, A and B;
// 3: Calculate A from P, Q, X, Y and B;
// 4: Calculate B from P, Q, X, Y and A.
//
// Input/output, double *P, the integral from 0 to X of the
// chi-square distribution. Input range: [0, 1].
//
// Input/output, double *Q, equals 1-P. Input range: [0, 1].
//
// Input/output, double *X, the upper limit of integration
// of the beta density. If it is an input value, it should lie in
// the range [0,1]. If it is an output value, it will be searched for
// in the range [0,1].
//
// Input/output, double *Y, equal to 1-X. If it is an input
// value, it should lie in the range [0,1]. If it is an output value,
// it will be searched for in the range [0,1].
//
// Input/output, double *A, the first parameter of the beta
// density. If it is an input value, it should lie in the range
// (0, +infinity). If it is an output value, it will be searched
// for in the range [1D-300,1D300].
//
// Input/output, double *B, the second parameter of the beta
// density. If it is an input value, it should lie in the range
// (0, +infinity). If it is an output value, it will be searched
// for in the range [1D-300,1D300].
//
// Output, int *STATUS
//
// Output, double *BOUND
public func betaCDF(xx: Float64, aa: Float64, bb: Float64): Float64 {
var t1 = unsafe {malloc(UIntNative(8))}
var x = unsafe { CPointer<Float64>(t1) }
unsafe {x.write(xx)}
var t2 = unsafe {malloc(UIntNative(8))}
var a = unsafe { CPointer<Float64>(t2) }
unsafe {a.write(aa)}
var t3 = unsafe {malloc(UIntNative(8))}
var b = unsafe { CPointer<Float64>(t3) }
unsafe {b.write(bb)}
var t4 = unsafe {malloc(UIntNative(8))}
var which = unsafe { CPointer<Int32>(t4) }
unsafe {which.write(1)}
var t5 = unsafe {malloc(UIntNative(8))}
var p = unsafe { CPointer<Float64>(t5) }
unsafe {p.write(0.0)}
var t6 = unsafe {malloc(UIntNative(8))}
var q = unsafe { CPointer<Float64>(t6) }
unsafe {q.write(0.0)}
var t7 = unsafe {malloc(UIntNative(8))}
var y = unsafe { CPointer<Float64>(t7) }
unsafe {y.write(1.0 - xx)}
var t8 = unsafe {malloc(UIntNative(8))}
var status = unsafe { CPointer<Int32>(t8) }
unsafe {status.write(0)}
var t9 = unsafe {malloc(UIntNative(8))}
var bound = unsafe { CPointer<Float64>(t9) }
unsafe {bound.write(0.0)}
unsafe {cdfbet (which, p, q, x, y, a, b, status, bound)}
return unsafe {p.read()}
}
public func betaPPF(pp: Float64, aa: Float64, bb: Float64): Float64 {
var t1 = unsafe {malloc(UIntNative(8))}
var x = unsafe { CPointer<Float64>(t1) }
unsafe {x.write(0.0)}
var t2 = unsafe {malloc(UIntNative(8))}
var a = unsafe { CPointer<Float64>(t2) }
unsafe {a.write(aa)}
var t3 = unsafe {malloc(UIntNative(8))}
var b = unsafe { CPointer<Float64>(t3) }
unsafe {b.write(bb)}
var t4 = unsafe {malloc(UIntNative(8))}
var which = unsafe { CPointer<Int32>(t4) }
unsafe {which.write(2)}
var t5 = unsafe {malloc(UIntNative(8))}
var p = unsafe { CPointer<Float64>(t5) }
unsafe {p.write(pp)}
var t6 = unsafe {malloc(UIntNative(8))}
var q = unsafe { CPointer<Float64>(t6) }
unsafe {q.write(1.0 - pp)}
var t7 = unsafe {malloc(UIntNative(8))}
var y = unsafe { CPointer<Float64>(t7) }
unsafe {y.write(0.0)}
var t8 = unsafe {malloc(UIntNative(8))}
var status = unsafe { CPointer<Int32>(t8) }
unsafe {status.write(0)}
var t9 = unsafe {malloc(UIntNative(8))}
var bound = unsafe { CPointer<Float64>(t9) }
unsafe {bound.write(0.0)}
unsafe {cdfbet (which, p, q, x, y, a, b, status, bound)}
return unsafe {x.read()}
}
/*
* Sampling from the beta distribution
*/
public func betaSample(r: Random, a: Float64, b: Float64): Float64 {
let x = r.nextFloat64()
return betaPPF(x, a, b)
}
@Test
public class TestBeta {
@TestCase
func testPDF(): Unit {
@Assert(approxEqual(betaPDF(0.1, 1.0, 2.0), 1.7999999999999998, atol:1e-13))
@Assert(approxEqual(betaPDF(0.3, 1.0, 2.0), 1.4000000000000001, atol:1e-13))
@Assert(approxEqual(betaPDF(0.5, 1.0, 2.0), 1.0 , atol:1e-13))
@Assert(approxEqual(betaPDF(0.7, 1.0, 2.0), 0.6000000000000001, atol:1e-13))
@Assert(approxEqual(betaPDF(0.9, 1.0, 2.0), 0.19999999999999998, atol:1e-13))
}
@TestCase
func testCDF(): Unit {
@Assert(approxEqual(betaCDF(0.1, 1.0, 2.0), 0.19, atol:1e-13))
@Assert(approxEqual(betaCDF(0.3, 1.0, 2.0), 0.51, atol:1e-13))
@Assert(approxEqual(betaCDF(0.5, 1.0, 2.0), 0.75 , atol:1e-13))
@Assert(approxEqual(betaCDF(0.7, 1.0, 2.0), 0.9099999999999999, atol:1e-13))
@Assert(approxEqual(betaCDF(0.9, 1.0, 2.0), 0.99, atol:1e-13))
}
@TestCase
func testPPF(): Unit {
@Assert(approxEqual(betaPPF(0.1, 1.0, 2.0), 0.051316701949486204, atol:1e-6))
@Assert(approxEqual(betaPPF(0.3, 1.0, 2.0), 0.16333997346592444, atol:1e-6))
@Assert(approxEqual(betaPPF(0.5, 1.0, 2.0), 0.2928932188134525, atol:1e-6))
@Assert(approxEqual(betaPPF(0.7, 1.0, 2.0), 0.45227744249483387, atol:1e-6))
@Assert(approxEqual(betaPPF(0.9, 1.0, 2.0), 0.683772233983162, atol:1e-6))
}
}