Newer
Older
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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
# Auto covariance
#=
function autocov(x::AbstractVector, lag::Int)
m = mean(x)
cov = (x[lag+1:end] .- m) .* (x[1:end-lag] .- m)
cor = (x .- m) .^ 2
return sum(cov)/sum(cor)
end
=#
function autocov(c::RejectionChains)
X = states(c)
C = [ autocov(X[i]) for i=1:nchains(c) ]
return C
end
function autocov(c::RejectionChains, lags::AbstractVector{<:AbstractVector})
X = states(c)
C = [ autocov(X[i], lags[i]) for i=1:nchains(c) ]
return C
end
function autocov(c::RejectionChains, lags::AbstractVector{<:Integer})
X = states(c)
C = [ autocov(X[i], lags) for i=1:nchains(c) ]
return C
end
# Log partition function, constant factor for normalize log density
function logpart(logp::AbstractVector)
return -log( sum( exp.(-logp) ) )
end
function logpart(c::RejectionChains)
# single level
if ! is_multilevel(c)
reps = repetitions(c)
return [
(
logprobs = getfield.(c.samples[i], :logprob);
logprobs = vcat( fill.(logprobs, reps[i])... );
logpart(logprobs)
)
for i=1:nchains(c)
]
# multi level
else
reps = repetitions(c)
lvls = levels(c); maxlvl = maximum(lvls)
return [
[ # return one constant per level
(
idx = lvls[i] <= l;
logprobs = getindex.(getfield.(c.samples[i][idx], :logprob), l);
logprobs = vcat( fill.(logprobs, reps[i][idx])... );
logpart(logprobs)
)
for l = 1:maxlvl
]
for i=1:nchains(c)
]
end
end
function empirical_cdf(c::RejectionChains, x)
x = convert(eltype(c), x)
if ! is_multilevel(c)
reps = repetitions(c)
return [
(
s = getfield.(c.samples[i], :state);
sum( reps[i] .* [ all(X .<= x) for X in s] ) / sum(reps[i])
)
for i=1:nchains(c)
]
else
reps = repetitions(c)
lvls = levels(c); maxlvl = maximum(lvls)
cdfs = [
(
idx = lvls[i] .<= 1;
s = getfield.(c.samples[i], :state)[idx];
[ sum(reps[i][idx] .* [ all(X .<= x) for X in s]) / sum(reps[i][idx]) ]
)
for i=1:nchains(c)
]
sizehint!.(cdfs, maxlvl)
logparts = logpart(c)
#Z = mean(logparts) # average over chain, otherwise for each chain seperately
for l=2:maxlvl
for i=1:nchains(c)
idx = lvls[i] <= l
s = getfield.(c.samples[i], :state)[idx]; # ?????
f = getfield.(c.samples[i], :logprob)[idx];
f0 = getindex.(f, l-1)
f1 = getindex.(f, l)
Z0,Z1 = logparts[i][l-1:l]
_cdf = sum(
reps[i][idx] .* [ all(X .<= x) for x in s ]
.* (1 .- exp.(f1 - f0 - Z1 + Z0))
) / sum(reps[i][idx])
push!(cdfs[i], _cdf)
end
end
return cdfs
end
end