-
Notifications
You must be signed in to change notification settings - Fork 0
/
synmodps.jl
50 lines (40 loc) · 1.44 KB
/
synmodps.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
using Distributions
using StatsBase
using ArgCheck
mutable struct SyntheticModelPS
tt::Vector{Float64}
mixtures::Matrix{NegativeBinomialMixture}
end
function SyntheticModelPS(tt::AbstractVector{<:Number}, ncomps::AbstractVector{Int})
mixtures = [ NegativeBinomialMixture(fill(NegativeBinomial(0.001, 0.999), ncomps[i]),
Categorical(ncomps[i]))
for i in 1:length(ncomps), j in 1:length(tt) ]
SyntheticModelPS(tt, mixtures)
end
function StatsBase.fit!(syn::SyntheticModelPS, data::AbstractArray{Int,3}; threads=true, kwargs...)
d, T = size(syn.mixtures)
@argcheck d == size(data,2)
@argcheck T == size(data,1)
if !threads
for i in 1:d
for t in 1:T
fit!(syn.mixtures[i,t], @view data[t,i,:]; kwargs...)
isvalidmixture(syn.mixtures[i,t]) || return false
end
end
else
Threads.@threads for i in 1:d
for t in 1:T
fit!(syn.mixtures[i,t], @view data[t,i,:]; kwargs...)
isvalidmixture(syn.mixtures[i,t]) || return false
end
end
end
return true
end
function StatsBase.loglikelihood(syn::SyntheticModelPS, data::AbstractArray{Int,3})
d, T = size(syn.mixtures)
@argcheck d == size(data,2)
@argcheck T == size(data,1)
sum(sum.(logpdf(syn.mixtures[i,t], (@view data[t,i,:])) for i in 1:d, t in 1:T))
end