Created
May 5, 2020 08:30
-
-
Save nightscape/07459a0f77466a8934aac735157db072 to your computer and use it in GitHub Desktop.
A not yet working Scala implementation of Sequential Random Sampling from Vitter
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// Translated from https://github.com/gliese1337/vitter-sample/blob/master/src/index.ts | |
// TODO: The returns were yields and the entire algorithm should create a Seq or Iterable of Ints | |
import scala.util.control.Breaks._ | |
object SequentialRandomSampling { | |
val negalphainv = -13 | |
def skip(_k: Int, _N: Int) { | |
var k = _k | |
var N = _N | |
var qu1 = N - k + 1 | |
var S = 0 | |
var threshold = -negalphainv * k | |
var kinv = 1 / k | |
var Vprime = Math.pow(Math.random(), kinv) | |
var X = 0.0 | |
while (k > 1 && threshold < N) { | |
println(s"calculating skip for k = ${k} and N = ${N} in Method D") | |
val kmin1inv = 1 / (k - 1) | |
breakable { | |
while (true) { | |
breakable { | |
while (true) { | |
X = N * (1 - Vprime) | |
S = Math.floor(X).toInt | |
if (S < qu1) { | |
break | |
} | |
Vprime = Math.pow(Math.random(), kinv) | |
} | |
} | |
val y1 = Math.pow(Math.random() * N / qu1, kmin1inv) | |
Vprime = y1 * (1 - X / N) * (qu1 / (qu1 - S)) | |
if (Vprime <= 1) { | |
break | |
} | |
var y2 = 1 | |
var top = N - 1 | |
var bottom: Int = 0 | |
var limit: Int = 0 | |
if (k - 1 > S) { | |
bottom = N - k | |
limit = N - S | |
} else { | |
bottom = N - S - 1 | |
limit = qu1 | |
} | |
for (t <- (N - 1 until limit step -1)) { | |
y2 = (y2 * top) / bottom | |
top -= 1 | |
bottom -= 1 | |
} | |
if (N / (N - X) >= y1 * Math.pow(y2, kmin1inv)) { | |
Vprime = Math.pow(Math.random(), kmin1inv) | |
break | |
} | |
Vprime = Math.pow(Math.random(), kinv) | |
} | |
return S | |
} | |
N -= S + 1 | |
k -= 1 | |
kinv = kmin1inv | |
qu1 -= S | |
threshold += negalphainv | |
} | |
if (k > 1) { // Method A | |
var top = N - k | |
while (k >= 2) { | |
val V = Math.random() | |
var quot = top / N | |
S = 0 | |
while (quot > V) { | |
S += 1 | |
top -= 1 | |
N -= 1 | |
quot *= top / N | |
} | |
return S | |
N -= 1 | |
k -= 1 | |
} | |
return Math.floor(N * Math.random()); | |
} else { | |
return Math.floor(N * Vprime) | |
} | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment