-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathRv.jl
68 lines (54 loc) · 1.57 KB
/
Rv.jl
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
65
66
67
68
module Rv
export RV, sampleRV, sampleRVTrace, sampleRVTimes, rejectionSampler
using Distributions
#shortcut
flip(x) = Bernoulli(x)
type RV
distgenerator::Function
name #String or Char
parents
end
function sampleRV(r::RV)
# generate a single sample from scratch
return sampleRV(r, Dict())
end
function sampleRV(r::RV, d::Dict{Any,Any})
# generate a single sample but recycle already sampled values
key = r.name
if !haskey(d, key)
vals = map((s)->sampleRV(s,d), r.parents)
d[key] = rand(r.distgenerator(vals...))
end
return d[key]
end
function sampleRV(rs::Array{RV,1}, d=Dict())
# sample simultaneously from several random variables
return [sampleRV(r,d) for r in rs]
end
#samples a given RV n times. An array is returned, which containes all sampled values as Float64
function sampleRVTimes(r::RV, n :: Integer)
result = Array(Float64,n)
for i in 1:n
result[i] = sampleRV(r)
end
result
end
#a rejectionSampler that takes the maximum value of the samples as the envelope function.
# The Rv is sampled n-times. An array with accepted values will be returned
# M is a factor for tuning the envelope function c.
#It can be passed to the hist() function to see the PDF
function rejectionSampler(r::RV, M, c :: Function,n=1000::Int64)
f = sampleRVTimes(r,n)
#g funktion
result = Float64[]
#c * g(x) müssen immer >= f(x) sein
i = 1
while i < n
if rand() <= abs((val = sampleRV(r))/ (M*c(i)))
push!(result,val)
end
i += 1
end
result
end
end