From 16385b5f4e8e15fcb27774ff41b5443be0cd9642 Mon Sep 17 00:00:00 2001 From: mrFendel Date: Fri, 5 May 2023 18:45:54 +0300 Subject: [PATCH] -- refactoring --- .../{series => distributions}/zSNormalCdf.kt | 12 +-- .../kscience/kmath/series/SeriesAlgebra.kt | 2 +- .../kmath/series/VarianceRatioTest.kt | 89 +++++++++---------- 3 files changed, 51 insertions(+), 52 deletions(-) rename kmath-stat/src/commonMain/kotlin/space/kscience/kmath/{series => distributions}/zSNormalCdf.kt (79%) diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/zSNormalCdf.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/zSNormalCdf.kt similarity index 79% rename from kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/zSNormalCdf.kt rename to kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/zSNormalCdf.kt index 1eb733ea9..adfbbc9cd 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/zSNormalCdf.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/zSNormalCdf.kt @@ -3,19 +3,19 @@ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. */ -package space.kscience.kmath.series +package space.kscience.kmath.distributions import space.kscience.kmath.operations.DoubleField.pow import kotlin.math.PI import kotlin.math.absoluteValue import kotlin.math.exp -public fun zSNormalCDF(x: Double): Double { - /** - * Zelen & Severo approximation for the standard normal CDF. - * The error is bounded by 7.5 * 10e-8. - * */ +/** + * Zelen & Severo approximation for the standard normal CDF. + * The error is bounded by 7.5 * 10e-8. + * */ +public fun zSNormalCDF(x: Double): Double { val t = 1 / (1 + 0.2316419 * x.absoluteValue) val summ = 0.319381530*t - 0.356563782*t.pow(2) + 1.781477937*t.pow(3) - 1.821255978*t.pow(4) + 1.330274429*t.pow(5) diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt index 72b9ca7a2..cabff25e6 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt @@ -216,7 +216,7 @@ public open class SeriesAlgebra, out BA : BufferAlgebra override fun multiply(left: Buffer, right: Buffer): Buffer = left.zip(right) { l, r -> l * r } - public fun Buffer.diff(shift: Int=1): Buffer = this.zipWithShift(shift) {l, r -> r - l} + public fun Buffer.difference(shift: Int=1): Buffer = this.zipWithShift(shift) {l, r -> r - l} public companion object } diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatioTest.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatioTest.kt index 280a9db64..3c6747ed2 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatioTest.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatioTest.kt @@ -5,6 +5,7 @@ package space.kscience.kmath.series +import space.kscience.kmath.distributions.zSNormalCDF import space.kscience.kmath.operations.DoubleField.pow import space.kscience.kmath.operations.algebra import space.kscience.kmath.operations.bufferAlgebra @@ -12,62 +13,60 @@ import space.kscience.kmath.operations.fold import kotlin.math.absoluteValue -// TODO: add p-value with formula: 2*(1 - cdf(|zScore|)) +/** + * Container class for Variance Ratio Test result: + * ratio itself, corresponding Z-score, also it's p-value + * **/ public data class VarianceRatioTestResult(val varianceRatio: Double=1.0, val zScore: Double=0.0, val pValue: Double=0.5) - /** - * Container class for Variance Ratio Test result: - * ratio itself, corresponding Z-score, also it's p-value - * **/ -public fun varianceRatioTest(series: Series, shift: Int, homoscedastic: Boolean=true): VarianceRatioTestResult { - /** - * Calculates the Z-statistic and the p-value for the Lo and MacKinlay's Variance Ratio test (1987) - * under Homoscedastic or Heteroscedstic assumptions - * with two-sided p-value test - * https://ssrn.com/abstract=346975 - * **/ +/** + * Calculates the Z-statistic and the p-value for the Lo and MacKinlay's Variance Ratio test (1987) + * under Homoscedastic or Heteroscedstic assumptions + * with two-sided p-value test + * https://ssrn.com/abstract=346975 + * **/ +public fun SeriesAlgebra.varianceRatioTest(series: Series, shift: Int, homoscedastic: Boolean=true): VarianceRatioTestResult { require(shift > 1) {"Shift must be greater than one"} require(shift < series.size) {"Shift must be smaller than sample size"} val sum = { x: Double, y: Double -> x + y } - with(Double.algebra.bufferAlgebra.seriesAlgebra()) { - val mean = series.fold(0.0, sum) / series.size - val demeanedSquares = series.map { power(it - mean, 2) } - val variance = demeanedSquares.fold(0.0, sum) - if (variance == 0.0) return VarianceRatioTestResult() + + val mean = series.fold(0.0, sum) / series.size + val demeanedSquares = series.map { (it - mean).pow(2) } + val variance = demeanedSquares.fold(0.0, sum) + if (variance == 0.0) return VarianceRatioTestResult() - var seriesAgg = series - for (i in 1.. v1 + v2 } - } - - val demeanedSquaresAgg = seriesAgg.map { power(it - shift * mean, 2) } - val varianceAgg = demeanedSquaresAgg.fold(0.0, sum) - - val varianceRatio = - varianceAgg * (series.size.toDouble() - 1) / variance / (series.size.toDouble() - shift.toDouble() + 1) / (1 - shift.toDouble()/series.size.toDouble()) / shift.toDouble() - - - // calculating asymptotic variance - val phi = if (homoscedastic) { // under homoscedastic null hypothesis - 2 * (2 * shift - 1.0) * (shift - 1.0) / (3 * shift * series.size) - } else { // under heteroscedastic null hypothesis - var accumulator = 0.0 - for (j in 1.. v1 * v2 }.fold(0.0, sum) / variance.pow(2) - accumulator += delta * 4 * (shift - j).toDouble().pow(2) / shift.toDouble().pow(2) - } - accumulator - } - - val zScore = (varianceRatio - 1) / phi.pow(0.5) - val pValue = 2*(1 - zSNormalCDF(zScore.absoluteValue)) - return VarianceRatioTestResult(varianceRatio, zScore, pValue) + var seriesAgg = series + for (i in 1.. v1 + v2 } } + + val demeanedSquaresAgg = seriesAgg.map { (it - shift * mean).pow(2) } + val varianceAgg = demeanedSquaresAgg.fold(0.0, sum) + + val varianceRatio = + varianceAgg * (series.size.toDouble() - 1) / variance / (series.size.toDouble() - shift.toDouble() + 1) / (1 - shift.toDouble()/series.size.toDouble()) / shift.toDouble() + + + // calculating asymptotic variance + val phi = if (homoscedastic) { // under homoscedastic null hypothesis + 2 * (2 * shift - 1.0) * (shift - 1.0) / (3 * shift * series.size) + } else { // under heteroscedastic null hypothesis + var accumulator = 0.0 + for (j in 1.. v1 * v2 }.fold(0.0, sum) / variance.pow(2) + accumulator += delta * 4 * (shift - j).toDouble().pow(2) / shift.toDouble().pow(2) + } + accumulator + } + + val zScore = (varianceRatio - 1) / phi.pow(0.5) + val pValue = 2*(1 - zSNormalCDF(zScore.absoluteValue)) + return VarianceRatioTestResult(varianceRatio, zScore, pValue) }