kscience/kmath
9
0
Fork 1
Code Issues Pull Requests 6 Packages Projects Releases 16 Wiki Activity

Dev #280

Merged
altavir merged 99 commits from dev into master 2021-04-16 19:45:55 +03:00
Conversation 0 Commits 99 Files Changed 283 +334 -31
4 changed files with 334 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files
Showing only changes of commit d4aa4587a9 - Show all commits

11
kmath-prob/src/commonMain/kotlin/kscience/kmath/prob/internal/InternalErf.kt Normal file
View File

@ -0,0 +1,11 @@
package kscience.kmath.prob.internal
import kotlin.math.abs
internal object InternalErf {
fun erfc(x: Double): Double {
if (abs(x) > 40) return if (x > 0) 0.0 else 2.0
val ret = InternalGamma.regularizedGammaQ(0.5, x * x, 10000)
return if (x < 0) 2 - ret else ret
}
}

245
kmath-prob/src/commonMain/kotlin/kscience/kmath/prob/internal/InternalGamma.kt Normal file
View File

@ -0,0 +1,245 @@
package kscience.kmath.prob.internal
import kotlin.math.*
private abstract class ContinuedFraction protected constructor() {
protected abstract fun getA(n: Int, x: Double): Double
protected abstract fun getB(n: Int, x: Double): Double
fun evaluate(x: Double, maxIterations: Int): Double {
val small = 1e-50
var hPrev = getA(0, x)
if (hPrev == 0.0 || abs(0.0 - hPrev) <= small) hPrev = small
var n = 1
var dPrev = 0.0
var cPrev = hPrev
var hN = hPrev
while (n < maxIterations) {
val a = getA(n, x)
val b = getB(n, x)
var dN = a + b * dPrev
if (dN == 0.0 || abs(0.0 - dN) <= small) dN = small
var cN = a + b / cPrev
if (cN == 0.0 || abs(0.0 - cN) <= small) cN = small
dN = 1 / dN
val deltaN = cN * dN
hN = hPrev * deltaN
check(!hN.isInfinite()) { "hN is infinite" }
check(!hN.isNaN()) { "hN is NaN" }
if (abs(deltaN - 1.0) < 10e-9) break
dPrev = dN
cPrev = cN
hPrev = hN
n++
}
check(n < maxIterations) { "n is more than maxIterations" }
return hN
}
}
internal object InternalGamma {
const val LANCZOS_G = 607.0 / 128.0
private val LANCZOS = doubleArrayOf(
0.99999999999999709182,
57.156235665862923517,
-59.597960355475491248,
14.136097974741747174,
-0.49191381609762019978,
.33994649984811888699e-4,
.46523628927048575665e-4,
-.98374475304879564677e-4,
.15808870322491248884e-3,
-.21026444172410488319e-3,
.21743961811521264320e-3,
-.16431810653676389022e-3,
.84418223983852743293e-4,
-.26190838401581408670e-4,
.36899182659531622704e-5
)
private val HALF_LOG_2_PI = 0.5 * ln(2.0 * PI)
private const val INV_GAMMA1P_M1_A0 = .611609510448141581788E-08
private const val INV_GAMMA1P_M1_A1 = .624730830116465516210E-08
private const val INV_GAMMA1P_M1_B1 = .203610414066806987300E+00
private const val INV_GAMMA1P_M1_B2 = .266205348428949217746E-01
private const val INV_GAMMA1P_M1_B3 = .493944979382446875238E-03
private const val INV_GAMMA1P_M1_B4 = -.851419432440314906588E-05
private const val INV_GAMMA1P_M1_B5 = -.643045481779353022248E-05
private const val INV_GAMMA1P_M1_B6 = .992641840672773722196E-06
private const val INV_GAMMA1P_M1_B7 = -.607761895722825260739E-07
private const val INV_GAMMA1P_M1_B8 = .195755836614639731882E-09
private const val INV_GAMMA1P_M1_P0 = .6116095104481415817861E-08
private const val INV_GAMMA1P_M1_P1 = .6871674113067198736152E-08
private const val INV_GAMMA1P_M1_P2 = .6820161668496170657918E-09
private const val INV_GAMMA1P_M1_P3 = .4686843322948848031080E-10
private const val INV_GAMMA1P_M1_P4 = .1572833027710446286995E-11
private const val INV_GAMMA1P_M1_P5 = -.1249441572276366213222E-12
private const val INV_GAMMA1P_M1_P6 = .4343529937408594255178E-14
private const val INV_GAMMA1P_M1_Q1 = .3056961078365221025009E+00
private const val INV_GAMMA1P_M1_Q2 = .5464213086042296536016E-01
private const val INV_GAMMA1P_M1_Q3 = .4956830093825887312020E-02
private const val INV_GAMMA1P_M1_Q4 = .2692369466186361192876E-03
private const val INV_GAMMA1P_M1_C = -.422784335098467139393487909917598E+00
private const val INV_GAMMA1P_M1_C0 = .577215664901532860606512090082402E+00
private const val INV_GAMMA1P_M1_C1 = -.655878071520253881077019515145390E+00
private const val INV_GAMMA1P_M1_C2 = -.420026350340952355290039348754298E-01
private const val INV_GAMMA1P_M1_C3 = .166538611382291489501700795102105E+00
private const val INV_GAMMA1P_M1_C4 = -.421977345555443367482083012891874E-01
private const val INV_GAMMA1P_M1_C5 = -.962197152787697356211492167234820E-02
private const val INV_GAMMA1P_M1_C6 = .721894324666309954239501034044657E-02
private const val INV_GAMMA1P_M1_C7 = -.116516759185906511211397108401839E-02
private const val INV_GAMMA1P_M1_C8 = -.215241674114950972815729963053648E-03
private const val INV_GAMMA1P_M1_C9 = .128050282388116186153198626328164E-03
private const val INV_GAMMA1P_M1_C10 = -.201348547807882386556893914210218E-04
private const val INV_GAMMA1P_M1_C11 = -.125049348214267065734535947383309E-05
private const val INV_GAMMA1P_M1_C12 = .113302723198169588237412962033074E-05
private const val INV_GAMMA1P_M1_C13 = -.205633841697760710345015413002057E-06
fun logGamma(x: Double): Double {
val ret: Double
when {
x.isNaN() || x <= 0.0 -> ret = Double.NaN
x < 0.5 -> return logGamma1p(x) - ln(x)
x <= 2.5 -> return logGamma1p(x - 0.5 - 0.5)
x <= 8.0 -> {
val n = floor(x - 1.5).toInt()
var prod = 1.0
(1..n).forEach { i -> prod *= x - i }
return logGamma1p(x - (n + 1)) + ln(prod)
}
else -> {
val tmp = x + LANCZOS_G + .5
ret = (x + .5) * ln(tmp) - tmp + HALF_LOG_2_PI + ln(lanczos(x) / x)
}
}
return ret
}
private fun regularizedGammaP(
a: Double,
x: Double,
maxIterations: Int = Int.MAX_VALUE
): Double = when {
a.isNaN() || x.isNaN() || a <= 0.0 || x < 0.0 -> Double.NaN
x == 0.0 -> 0.0
x >= a + 1 -> 1.0 - regularizedGammaQ(a, x, maxIterations)
else -> {
// calculate series
var n = 0.0 // current element index
var an = 1.0 / a // n-th element in the series
var sum = an // partial sum
while (abs(an / sum) > 10e-15 && n < maxIterations && sum < Double.POSITIVE_INFINITY) {
// compute next element in the series
n += 1.0
an *= x / (a + n)
// update partial sum
sum += an
}
when {
n >= maxIterations -> throw error("Maximal iterations is exceeded $maxIterations")
sum.isInfinite() -> 1.0
else -> exp(-x + a * ln(x) - logGamma(a)) * sum
}
}
}
fun regularizedGammaQ(
a: Double,
x: Double,
maxIterations: Int = Int.MAX_VALUE
): Double = when {
a.isNaN() || x.isNaN() || a <= 0.0 || x < 0.0 -> Double.NaN
x == 0.0 -> 1.0
x < a + 1.0 -> 1.0 - regularizedGammaP(a, x, maxIterations)
else -> 1.0 / object : ContinuedFraction() {
override fun getA(n: Int, x: Double): Double = 2.0 * n + 1.0 - a + x
override fun getB(n: Int, x: Double): Double = n * (a - n)
}.evaluate(x, maxIterations) * exp(-x + a * ln(x) - logGamma(a))
}
private fun lanczos(x: Double): Double =
(LANCZOS.size - 1 downTo 1).sumByDouble { LANCZOS[it] / (x + it) } + LANCZOS[0]
private fun invGamma1pm1(x: Double): Double {
require(x >= -0.5)
require(x <= 1.5)
val ret: Double
val t = if (x <= 0.5) x else x - 0.5 - 0.5
if (t < 0.0) {
val a = INV_GAMMA1P_M1_A0 + t * INV_GAMMA1P_M1_A1
var b = INV_GAMMA1P_M1_B8
b = INV_GAMMA1P_M1_B7 + t * b
b = INV_GAMMA1P_M1_B6 + t * b
b = INV_GAMMA1P_M1_B5 + t * b
b = INV_GAMMA1P_M1_B4 + t * b
b = INV_GAMMA1P_M1_B3 + t * b
b = INV_GAMMA1P_M1_B2 + t * b
b = INV_GAMMA1P_M1_B1 + t * b
b = 1.0 + t * b
var c = INV_GAMMA1P_M1_C13 + t * (a / b)
c = INV_GAMMA1P_M1_C12 + t * c
c = INV_GAMMA1P_M1_C11 + t * c
c = INV_GAMMA1P_M1_C10 + t * c
c = INV_GAMMA1P_M1_C9 + t * c
c = INV_GAMMA1P_M1_C8 + t * c
c = INV_GAMMA1P_M1_C7 + t * c
c = INV_GAMMA1P_M1_C6 + t * c
c = INV_GAMMA1P_M1_C5 + t * c
c = INV_GAMMA1P_M1_C4 + t * c
c = INV_GAMMA1P_M1_C3 + t * c
c = INV_GAMMA1P_M1_C2 + t * c
c = INV_GAMMA1P_M1_C1 + t * c
c = INV_GAMMA1P_M1_C + t * c
ret = (if (x > 0.5) t * c / x else x * (c + 0.5 + 0.5))
} else {
var p = INV_GAMMA1P_M1_P6
p = INV_GAMMA1P_M1_P5 + t * p
p = INV_GAMMA1P_M1_P4 + t * p
p = INV_GAMMA1P_M1_P3 + t * p
p = INV_GAMMA1P_M1_P2 + t * p
p = INV_GAMMA1P_M1_P1 + t * p
p = INV_GAMMA1P_M1_P0 + t * p
var q = INV_GAMMA1P_M1_Q4
q = INV_GAMMA1P_M1_Q3 + t * q
q = INV_GAMMA1P_M1_Q2 + t * q
q = INV_GAMMA1P_M1_Q1 + t * q
q = 1.0 + t * q
var c = INV_GAMMA1P_M1_C13 + p / q * t
c = INV_GAMMA1P_M1_C12 + t * c
c = INV_GAMMA1P_M1_C11 + t * c
c = INV_GAMMA1P_M1_C10 + t * c
c = INV_GAMMA1P_M1_C9 + t * c
c = INV_GAMMA1P_M1_C8 + t * c
c = INV_GAMMA1P_M1_C7 + t * c
c = INV_GAMMA1P_M1_C6 + t * c
c = INV_GAMMA1P_M1_C5 + t * c
c = INV_GAMMA1P_M1_C4 + t * c
c = INV_GAMMA1P_M1_C3 + t * c
c = INV_GAMMA1P_M1_C2 + t * c
c = INV_GAMMA1P_M1_C1 + t * c
c = INV_GAMMA1P_M1_C0 + t * c
ret = (if (x > 0.5) t / x * (c - 0.5 - 0.5) else x * c)
}
return ret
}
private fun logGamma1p(x: Double): Double {
require(x >= -0.5)
require(x <= 1.5)
return -ln1p(invGamma1pm1(x))
}
}

78
kmath-prob/src/commonMain/kotlin/kscience/kmath/prob/internal/InternalUtils.kt Normal file
View File

@ -0,0 +1,78 @@
package kscience.kmath.prob.internal
import kotlin.math.ln
import kotlin.math.min
internal object InternalUtils {
private val FACTORIALS = longArrayOf(
1L, 1L, 2L,
6L, 24L, 120L,
720L, 5040L, 40320L,
362880L, 3628800L, 39916800L,
479001600L, 6227020800L, 87178291200L,
1307674368000L, 20922789888000L, 355687428096000L,
6402373705728000L, 121645100408832000L, 2432902008176640000L
)
private const val BEGIN_LOG_FACTORIALS = 2
fun factorial(n: Int): Long = FACTORIALS[n]
fun validateProbabilities(probabilities: DoubleArray?): Double {
require(!(probabilities == null || probabilities.isEmpty())) { "Probabilities must not be empty." }
var sumProb = 0.0
probabilities.forEach { prob ->
validateProbability(prob)
sumProb += prob
}
require(!(sumProb.isInfinite() || sumProb <= 0)) { "Invalid sum of probabilities: $sumProb" }
return sumProb
}
private fun validateProbability(probability: Double): Unit =
require(!(probability < 0 || probability.isInfinite() || probability.isNaN())) { "Invalid probability: $probability" }
class FactorialLog private constructor(
numValues: Int,
cache: DoubleArray?
) {
private val logFactorials: DoubleArray = DoubleArray(numValues)
init {
val endCopy: Int
if (cache != null && cache.size > BEGIN_LOG_FACTORIALS) {
// Copy available values.
endCopy = min(cache.size, numValues)
cache.copyInto(
logFactorials,
BEGIN_LOG_FACTORIALS,
BEGIN_LOG_FACTORIALS, endCopy
)
}
// All values to be computed
else endCopy = BEGIN_LOG_FACTORIALS
// Compute remaining values.
(endCopy until numValues).forEach { i ->
if (i < FACTORIALS.size)
logFactorials[i] = ln(FACTORIALS[i].toDouble())
else
logFactorials[i] = logFactorials[i - 1] + ln(i.toDouble())
}
}
fun value(n: Int): Double {
if (n < logFactorials.size)
return logFactorials[n]
return if (n < FACTORIALS.size) ln(FACTORIALS[n].toDouble()) else InternalGamma.logGamma(n + 1.0)
}
companion object {
fun create(): FactorialLog = FactorialLog(0, null)
}
}
}

31
kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/RandomSourceGenerator.kt
View File

@ -1,31 +0,0 @@
package scientifik.kmath.prob
import kscience.kmath.prob.RandomGenerator
import org.apache.commons.rng.simple.RandomSource
public class RandomSourceGenerator(private val source: RandomSource, seed: Long?) :
RandomGenerator {
private val random = seed?.let {
RandomSource.create(source, seed)
} ?: RandomSource.create(source)
public override fun nextBoolean(): Boolean = random.nextBoolean()
public override fun nextDouble(): Double = random.nextDouble()
public override fun nextInt(): Int = random.nextInt()
public override fun nextInt(until: Int): Int = random.nextInt(until)
public override fun nextLong(): Long = random.nextLong()
public override fun nextLong(until: Long): Long = random.nextLong(until)
public override fun fillBytes(array: ByteArray, fromIndex: Int, toIndex: Int) {
require(toIndex > fromIndex)
random.nextBytes(array, fromIndex, toIndex - fromIndex)
}
override fun fork(): RandomGenerator = RandomSourceGenerator(source, nextLong())
}
public fun RandomGenerator.Companion.fromSource(source: RandomSource, seed: Long? = null): RandomSourceGenerator =
RandomSourceGenerator(source, seed)
public fun RandomGenerator.Companion.mersenneTwister(seed: Long? = null): RandomSourceGenerator =
fromSource(RandomSource.MT, seed)