-
Notifications
You must be signed in to change notification settings - Fork 4
/
arsampling.py
64 lines (53 loc) · 2.13 KB
/
arsampling.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
import numpy as np
__all__ = ['rejectionSampling']
def rejectionSampling(envPDF, envRV, targetPDF, n):
"""Generates n random variates from the target distribution
using acceptance/rejection sampling with the envelope distribution, envPDF.
Parameters
----------
envPDF : func
Takes an array of x and returns their probabilities
times a constant, such that envPDF(x) > targetPDF(x), always.
envRV : func
Takes single arg n and returns n random variates from envPDF.
targetPDF : func
Takes array of x and returns their probabilities
n : int
Number of total samples to draw from targetPDF
Returns
-------
arr : array
Array of samples distributed proportionately to target PDF
acceptanceProbability : float
Estimate for diagnosing the envelope distribution."""
i = 0.
arr = np.zeros(n)
acceptanceProbability = 1.
allDraws = None
while i < n:
"""Draw as many samples as we expect to need,
based on the acceptance probability (initially 1)"""
neededN = n-i
drawingN = int(np.ceil(neededN * (1./acceptanceProbability)))
samples = envRV(drawingN)
envPr = envPDF(samples)
targetPr = targetPDF(samples)
unifRVs = np.random.rand(drawingN)
"""Keep the samples that meet the acceptance criteria"""
keepInds = unifRVs < (targetPr/envPr)
actualN = keepInds.sum()
if actualN < neededN:
"""Fill arr with all the samples that were accepted"""
arr[i:(i+actualN)] = samples[keepInds]
else:
"""Fill the rest of the arr with as many samples as needed"""
arr[i:] = samples[keepInds][:neededN]
"""Keep track of all draws to compute acceptanceProbability"""
if allDraws is None:
allDraws = targetPr/envPr
else:
allDraws = np.concatenate((allDraws, targetPr/envPr))
"""Estimate the acceptance probability"""
acceptanceProbability = allDraws.mean()
i += actualN
return arr, acceptanceProbability