-
-
Notifications
You must be signed in to change notification settings - Fork 300
/
discrete_rv.jl
80 lines (57 loc) · 1.83 KB
/
discrete_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
69
70
71
72
73
74
75
76
77
78
79
80
#=
Generates an array of draws from a discrete random variable with a
specified vector of probabilities.
@author : Spencer Lyon <spencer.lyon@nyu.edu>
@date: 2014-07-10
References
----------
https://lectures.quantecon.org/jl/finite_markov.html
TODO: as of 07/10/2014 it is not possible to define the property
interface we see in the python version. Once issue 1974 from
the main Julia repository is resolved, we can revisit this and
implement that feature.
=#
"""
Generates an array of draws from a discrete random variable with
vector of probabilities given by `q`.
##### Fields
- `q::AbstractVector`: A vector of non-negative probabilities that sum to 1
- `Q::AbstractVector`: The cumulative sum of `q`
"""
mutable struct DiscreteRV{TV1<:AbstractVector, TV2<:AbstractVector}
q::TV1
Q::TV2
function DiscreteRV{TV1,TV2}(q, Q) where {TV1,TV2}
abs(Q[end] - 1) > 1e-10 && error("q should sum to 1")
new{TV1,TV2}(q, Q)
end
end
function DiscreteRV(q::TV) where TV<:AbstractVector
Q = cumsum(q)
DiscreteRV{TV,typeof(Q)}(q, Q)
end
"""
Make a single draw from the discrete distribution.
##### Arguments
- `d::DiscreteRV`: The `DiscreteRV` type represetning the distribution
##### Returns
- `out::Int`: One draw from the discrete distribution
"""
Random.rand(d::DiscreteRV) = searchsortedfirst(d.Q, rand())
"""
Make multiple draws from the discrete distribution represented by a
`DiscreteRV` instance
##### Arguments
- `d::DiscreteRV`: The `DiscreteRV` type representing the distribution
- `k::Int`
##### Returns
- `out::Vector{Int}`: `k` draws from `d`
"""
Random.rand(d::DiscreteRV, k::Int) = Int[rand(d) for i=1:k]
function Random.rand!(out::AbstractArray{T}, d::DiscreteRV) where T<:Integer
@inbounds for I in eachindex(out)
out[I] = rand(d)
end
out
end
@deprecate draw Random.rand