Advertisement
Guest User

Untitled

a guest
Sep 12th, 2015
120
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 19.25 KB | None | 0 0
  1. - #=
  2. - DPM
  3. -
  4. - Dirichlet Process Mixture Models
  5. -
  6. - 25/08/2015
  7. - Adham Beyki, odinay@gmail.com
  8. -
  9. - =#
  10. -
  11. - type DPM{T}
  12. - bayesian_component::T
  13. - K::Int64
  14. - aa::Float64
  15. - a1::Float64
  16. - a2::Float64
  17. - K_hist::Vector{Int64}
  18. - K_zz_dict::Dict{Int64, Vector{Int64}}
  19. -
  20. - DPM{T}(c::T, K::Int64, aa::Float64, a1::Float64, a2::Float64) = new(c, K, aa, a1, a2,
  21. - Int64[], (Int64 => Vector{Int64})[])
  22. - end
  23. 1 DPM{T}(c::T, K::Int64, aa::Real, a1::Real, a2::Real) = DPM{typeof(c)}(c, K, convert(Float64, aa),
  24. - convert(Float64, a1), convert(Float64, a2))
  25. -
  26. - function Base.show(io::IO, dpm::DPM)
  27. - println(io, "Dirichlet Mixture Model with $(dpm.K) $(typeof(dpm.bayesian_component)) components")
  28. - end
  29. -
  30. - function initialize_gibbs_sampler!(dpm::DPM, zz::Vector{Int64})
  31. - # populates the cluster labels randomly
  32. 1 zz[:] = rand(1:dpm.K, length(zz))
  33. - end
  34. -
  35. - function DPM_sample_hyperparam(aa::Float64, a1::Float64, a2::Float64, K::Int64, NN::Int64, iters::Int64)
  36. -
  37. - # resampling concentration parameter based on Escobar and West 1995
  38. 352 for n = 1:iters
  39. 3504 eta = rand(Distributions.Beta(aa+1, NN))
  40. 3504 rr = (a1+K-1) / (NN*(a2-log(NN)))
  41. 3504 pi_eta = rr / (1+rr)
  42. -
  43. 3504 if rand() < pi_eta
  44. 0 aa = rand(Distributions.Gamma(a1+K, 1/(a2-log(eta))))
  45. - else
  46. 3504 aa = rand(Distributions.Gamma(a1+K-1, 1/(a2-log(eta))))
  47. - end
  48. - end
  49. 352 aa
  50. - end
  51. -
  52. - function DPM_sample_pp{T1, T2}(
  53. - bayesian_components::Vector{T1},
  54. - xx::T2,
  55. - nn::Vector{Float64},
  56. - pp::Vector{Float64},
  57. - aa::Float64)
  58. -
  59. 1760000 K = length(nn)
  60. 1760000 @inbounds for kk = 1:K
  61. 11384379 pp[kk] = log(nn[kk]) + logpredictive(bayesian_components[kk], xx)
  62. - end
  63. 1760000 pp[K+1] = log(aa) + logpredictive(bayesian_components[K+1], xx)
  64. 1760000 normalize_pp!(pp, K+1)
  65. 1760000 return sample(pp[1:K+1])
  66. - end
  67. -
  68. -
  69. - function collapsed_gibbs_sampler!{T1, T2}(
  70. - dpm::DPM{T1},
  71. - xx::Vector{T2},
  72. - zz::Vector{Int64},
  73. - n_burnins::Int64, n_lags::Int64, n_samples::Int64, n_internals::Int64; max_clusters::Int64=100)
  74. -
  75. -
  76. 2 NN = length(xx) # number of data points
  77. 2 nn = zeros(Float64, dpm.K) # count array
  78. 2 n_iterations = n_burnins + (n_samples)*(n_lags+1)
  79. 2 bayesian_components = [deepcopy(dpm.bayesian_component) for k = 1:dpm.K+1]
  80. 2 dpm.K_hist = zeros(Int64, n_iterations)
  81. 2 pp = zeros(Float64, max_clusters)
  82. -
  83. 2 tic()
  84. 2 for ii = 1:NN
  85. 10000 kk = zz[ii]
  86. 10000 additem(bayesian_components[kk], xx[ii])
  87. 10000 nn[kk] += 1
  88. - end
  89. 2 dpm.K_hist[1] = dpm.K
  90. 2 elapsed_time = toq()
  91. -
  92. 2 for iteration = 1:n_iterations
  93. -
  94. 352 println("iteration: $iteration, KK: $(dpm.K), KK mode: $(indmax(hist(dpm.K_hist,
  95. - 0.5:maximum(dpm.K_hist)+0.5)[2])), elapsed time: $elapsed_time")
  96. -
  97. 352 tic()
  98. 352 @inbounds for ii = 1:NN
  99. 1760000 kk = zz[ii]
  100. 1760000 nn[kk] -= 1
  101. 1760000 delitem(bayesian_components[kk], xx[ii])
  102. -
  103. - # remove the cluster if empty
  104. 1760000 if nn[kk] == 0
  105. 166 println("tcomponent $kk has become inactive")
  106. 166 splice!(nn, kk)
  107. 166 splice!(bayesian_components, kk)
  108. 166 dpm.K -= 1
  109. -
  110. - # shifting the labels one cluster back
  111. 830166 idx = find(x -> x>kk, zz)
  112. 166 zz[idx] -= 1
  113. - end
  114. -
  115. 1760000 kk = DPM_sample_pp(bayesian_components, xx[ii], nn, pp, dpm.aa)
  116. -
  117. 1760000 if kk == dpm.K+1
  118. 171 println("tcomponent $kk activated.")
  119. 171 push!(bayesian_components, deepcopy(dpm.bayesian_component))
  120. 171 push!(nn, 0)
  121. 171 dpm.K += 1
  122. - end
  123. -
  124. 1760000 zz[ii] = kk
  125. 1760000 nn[kk] += 1
  126. 1760000 additem(bayesian_components[kk], xx[ii])
  127. - end
  128. -
  129. 352 dpm.aa = DPM_sample_hyperparam(dpm.aa, dpm.a1, dpm.a2, dpm.K, NN, n_internals)
  130. 352 dpm.K_hist[iteration] = dpm.K
  131. 352 dpm.K_zz_dict[dpm.K] = deepcopy(zz)
  132. 352 elapsed_time = toq()
  133. - end
  134. - end
  135. -
  136. - function truncated_gibbs_sampler{T1, T2}(dpm::DPM{T1}, xx::Vector{T2}, zz::Vector{Int64},
  137. - n_burnins::Int64, n_lags::Int64, n_samples::Int64, n_internals::Int64, K_truncation::Int64)
  138. -
  139. - NN = length(xx) # number of data points
  140. - nn = zeros(Int64, K_truncation) # count array
  141. - bayesian_components = [deepcopy(dpm.bayesian_component) for k = 1:K_truncation]
  142. - n_iterations = n_burnins + (n_samples)*(n_lags+1)
  143. - dpm.K_hist = zeros(Int64, n_iterations)
  144. - states = (ASCIIString => Int64)[]
  145. - n_states = 0
  146. -
  147. - tic()
  148. - for ii = 1:NN
  149. - kk = zz[ii]
  150. - additem(bayesian_components[kk], xx[ii])
  151. - nn[kk] += 1
  152. - end
  153. - dpm.K_hist[1] = dpm.K
  154. -
  155. - # constructing the sticks
  156. - beta_VV = rand(Distributions.Beta(1.0, dpm.aa), K_truncation)
  157. - beta_VV[end] = 1.0
  158. - π = ones(Float64, K_truncation)
  159. - π[2:end] = 1 - beta_VV[1:K_truncation-1]
  160. - π = log(beta_VV) + log(cumprod(π))
  161. -
  162. - elapsed_time = toq()
  163. -
  164. - for iteration = 1:n_iterations
  165. -
  166. - println("iteration: $iteration, # active components: $(length(findn(nn)[1])), mode: $(indmax(hist(dpm.K_hist,
  167. - 0.5:maximum(dpm.K_hist)+0.5)[2])), elapsed time: $elapsed_time n", nn)
  168. -
  169. - tic()
  170. - for ii = 1:NN
  171. - kk = zz[ii]
  172. - nn[kk] -= 1
  173. - delitem(bayesian_components[kk], xx[ii])
  174. -
  175. - # resampling label
  176. - pp = zeros(Float64, K_truncation)
  177. - for kk = 1:K_truncation
  178. - pp[kk] = π[kk] + logpredictive(bayesian_components[kk], xx[ii])
  179. - end
  180. - pp = exp(pp - maximum(pp))
  181. - pp /= sum(pp)
  182. -
  183. - # sample from pp
  184. - kk = sampleindex(pp)
  185. - zz[ii] = kk
  186. - nn[kk] += 1
  187. - additem(bayesian_components[kk], xx[ii])
  188. -
  189. - for kk = 1:K_truncation-1
  190. - gamma1 = 1 + nn[kk]
  191. - gamma2 = dpm.aa + sum(nn[kk+1:end])
  192. - beta_VV[kk] = rand(Distributions.Beta(gamma1, gamma2))
  193. - end
  194. - beta_VV[end] = 1.0
  195. - π = ones(Float64, K_truncation)
  196. - π[2:end] = 1 - beta_VV[1:K_truncation-1]
  197. - π = log(beta_VV) + log(cumprod(π))
  198. -
  199. - # resampling concentration parameter based on Escobar and West 1995
  200. - for internal_iters = 1:n_internals
  201. - eta = rand(Distributions.Beta(dpm.aa+1, NN))
  202. - rr = (dpm.a1+dpm.K-1) / (NN*(dpm.a2-log(NN)))
  203. - pi_eta = rr / (1+rr)
  204. -
  205. - if rand() < pi_eta
  206. - dpm.aa = rand(Distributions.Gamma(dpm.a1+dpm.K, 1/(dpm.a2-log(eta))))
  207. - else
  208. - dpm.aa = rand(Distributions.Gamma(dpm.a1+dpm.K-1, 1/(dpm.a2-log(eta))))
  209. - end
  210. - end
  211. - end
  212. -
  213. - nn_string = nn2string(nn)
  214. - if !haskey(states, nn_string)
  215. - n_states += 1
  216. - states[nn_string] = n_states
  217. - end
  218. - dpm.K_hist[iteration] = states[nn_string]
  219. - dpm.K_zz_dict[states[nn_string]] = deepcopy(zz)
  220. - elapsed_time = toq()
  221. - end
  222. - return states
  223. - end
  224. -
  225. -
  226. - function posterior{T1, T2}(dpm::DPM{T1}, xx::Vector{T2}, K::Int64, K_truncation::Int64=0)
  227. 2 n_components = 0
  228. 1 if K_truncation == 0
  229. 1 n_components = K
  230. - else
  231. 0 n_components = K_truncation
  232. - end
  233. -
  234. 1 bayesian_components = [deepcopy(dpm.bayesian_component) for kk=1:n_components]
  235. 1 zz = dpm.K_zz_dict[K]
  236. -
  237. 1 NN = length(xx)
  238. 1 nn = zeros(Int64, n_components)
  239. -
  240. 1 for ii = 1:NN
  241. 5000 kk = zz[ii]
  242. 5000 additem(bayesian_components[kk], xx[ii])
  243. 5000 nn[kk] += 1
  244. - end
  245. -
  246. 1 return([posterior(bayesian_components[kk]) for kk=1:n_components], nn)
  247. - end
  248. -
  249.  
  250. - #=
  251. - DPM
  252. -
  253. - Dirichlet Process Mixture Models
  254. -
  255. - 25/08/2015
  256. - Adham Beyki, odinay@gmail.com
  257. -
  258. - =#
  259. -
  260. - type DPM{T}
  261. - bayesian_component::T
  262. - K::Int64
  263. - aa::Float64
  264. - a1::Float64
  265. - a2::Float64
  266. - K_hist::Vector{Int64}
  267. - K_zz_dict::Dict{Int64, Vector{Int64}}
  268. -
  269. - DPM{T}(c::T, K::Int64, aa::Float64, a1::Float64, a2::Float64) = new(c, K, aa, a1, a2,
  270. - Int64[], (Int64 => Vector{Int64})[])
  271. - end
  272. 0 DPM{T}(c::T, K::Int64, aa::Real, a1::Real, a2::Real) = DPM{typeof(c)}(c, K, convert(Float64, aa),
  273. - convert(Float64, a1), convert(Float64, a2))
  274. -
  275. - function Base.show(io::IO, dpm::DPM)
  276. - println(io, "Dirichlet Mixture Model with $(dpm.K) $(typeof(dpm.bayesian_component)) components")
  277. - end
  278. -
  279. - function initialize_gibbs_sampler!(dpm::DPM, zz::Vector{Int64})
  280. - # populates the cluster labels randomly
  281. 0 zz[:] = rand(1:dpm.K, length(zz))
  282. - end
  283. -
  284. - function DPM_sample_hyperparam(aa::Float64, a1::Float64, a2::Float64, K::Int64, NN::Int64, iters::Int64)
  285. -
  286. - # resampling concentration parameter based on Escobar and West 1995
  287. 0 for n = 1:iters
  288. 0 eta = rand(Distributions.Beta(aa+1, NN))
  289. 0 rr = (a1+K-1) / (NN*(a2-log(NN)))
  290. 0 pi_eta = rr / (1+rr)
  291. -
  292. 0 if rand() < pi_eta
  293. 0 aa = rand(Distributions.Gamma(a1+K, 1/(a2-log(eta))))
  294. - else
  295. 0 aa = rand(Distributions.Gamma(a1+K-1, 1/(a2-log(eta))))
  296. - end
  297. - end
  298. 0 aa
  299. - end
  300. -
  301. - function DPM_sample_pp{T1, T2}(
  302. - bayesian_components::Vector{T1},
  303. - xx::T2,
  304. - nn::Vector{Float64},
  305. - pp::Vector{Float64},
  306. - aa::Float64)
  307. -
  308. 0 K = length(nn)
  309. 0 @inbounds for kk = 1:K
  310. 0 pp[kk] = log(nn[kk]) + logpredictive(bayesian_components[kk], xx)
  311. - end
  312. 0 pp[K+1] = log(aa) + logpredictive(bayesian_components[K+1], xx)
  313. 0 normalize_pp!(pp, K+1)
  314. 0 return sample(pp[1:K+1])
  315. - end
  316. -
  317. -
  318. - function collapsed_gibbs_sampler!{T1, T2}(
  319. - dpm::DPM{T1},
  320. - xx::Vector{T2},
  321. - zz::Vector{Int64},
  322. - n_burnins::Int64, n_lags::Int64, n_samples::Int64, n_internals::Int64; max_clusters::Int64=100)
  323. -
  324. -
  325. 191688 NN = length(xx) # number of data points
  326. 96 nn = zeros(Float64, dpm.K) # count array
  327. 0 n_iterations = n_burnins + (n_samples)*(n_lags+1)
  328. 384 bayesian_components = [deepcopy(dpm.bayesian_component) for k = 1:dpm.K+1]
  329. 2864 dpm.K_hist = zeros(Int64, n_iterations)
  330. 176 pp = zeros(Float64, max_clusters)
  331. -
  332. 48 tic()
  333. 0 for ii = 1:NN
  334. 0 kk = zz[ii]
  335. 0 additem(bayesian_components[kk], xx[ii])
  336. 0 nn[kk] += 1
  337. - end
  338. 0 dpm.K_hist[1] = dpm.K
  339. 0 elapsed_time = toq()
  340. -
  341. 0 for iteration = 1:n_iterations
  342. -
  343. 5329296 println("iteration: $iteration, KK: $(dpm.K), KK mode: $(indmax(hist(dpm.K_hist,
  344. - 0.5:maximum(dpm.K_hist)+0.5)[2])), elapsed time: $elapsed_time")
  345. -
  346. 16800 tic()
  347. 28000000 @inbounds for ii = 1:NN
  348. 0 kk = zz[ii]
  349. 0 nn[kk] -= 1
  350. 0 delitem(bayesian_components[kk], xx[ii])
  351. -
  352. - # remove the cluster if empty
  353. 0 if nn[kk] == 0
  354. 161880 println("tcomponent $kk has become inactive")
  355. 0 splice!(nn, kk)
  356. 0 splice!(bayesian_components, kk)
  357. 0 dpm.K -= 1
  358. -
  359. - # shifting the labels one cluster back
  360. 69032 idx = find(x -> x>kk, zz)
  361. 42944 zz[idx] -= 1
  362. - end
  363. -
  364. 0 kk = DPM_sample_pp(bayesian_components, xx[ii], nn, pp, dpm.aa)
  365. -
  366. 0 if kk == dpm.K+1
  367. 158976 println("tcomponent $kk activated.")
  368. 14144 push!(bayesian_components, deepcopy(dpm.bayesian_component))
  369. 4872 push!(nn, 0)
  370. 0 dpm.K += 1
  371. - end
  372. -
  373. 0 zz[ii] = kk
  374. 0 nn[kk] += 1
  375. 0 additem(bayesian_components[kk], xx[ii])
  376. - end
  377. -
  378. 0 dpm.aa = DPM_sample_hyperparam(dpm.aa, dpm.a1, dpm.a2, dpm.K, NN, n_internals)
  379. 0 dpm.K_hist[iteration] = dpm.K
  380. 14140000 dpm.K_zz_dict[dpm.K] = deepcopy(zz)
  381. 0 elapsed_time = toq()
  382. - end
  383. - end
  384. -
  385. - function truncated_gibbs_sampler{T1, T2}(dpm::DPM{T1}, xx::Vector{T2}, zz::Vector{Int64},
  386. - n_burnins::Int64, n_lags::Int64, n_samples::Int64, n_internals::Int64, K_truncation::Int64)
  387. -
  388. - NN = length(xx) # number of data points
  389. - nn = zeros(Int64, K_truncation) # count array
  390. - bayesian_components = [deepcopy(dpm.bayesian_component) for k = 1:K_truncation]
  391. - n_iterations = n_burnins + (n_samples)*(n_lags+1)
  392. - dpm.K_hist = zeros(Int64, n_iterations)
  393. - states = (ASCIIString => Int64)[]
  394. - n_states = 0
  395. -
  396. - tic()
  397. - for ii = 1:NN
  398. - kk = zz[ii]
  399. - additem(bayesian_components[kk], xx[ii])
  400. - nn[kk] += 1
  401. - end
  402. - dpm.K_hist[1] = dpm.K
  403. -
  404. - # constructing the sticks
  405. - beta_VV = rand(Distributions.Beta(1.0, dpm.aa), K_truncation)
  406. - beta_VV[end] = 1.0
  407. - π = ones(Float64, K_truncation)
  408. - π[2:end] = 1 - beta_VV[1:K_truncation-1]
  409. - π = log(beta_VV) + log(cumprod(π))
  410. -
  411. - elapsed_time = toq()
  412. -
  413. - for iteration = 1:n_iterations
  414. -
  415. - println("iteration: $iteration, # active components: $(length(findn(nn)[1])), mode: $(indmax(hist(dpm.K_hist,
  416. - 0.5:maximum(dpm.K_hist)+0.5)[2])), elapsed time: $elapsed_time n", nn)
  417. -
  418. - tic()
  419. - for ii = 1:NN
  420. - kk = zz[ii]
  421. - nn[kk] -= 1
  422. - delitem(bayesian_components[kk], xx[ii])
  423. -
  424. - # resampling label
  425. - pp = zeros(Float64, K_truncation)
  426. - for kk = 1:K_truncation
  427. - pp[kk] = π[kk] + logpredictive(bayesian_components[kk], xx[ii])
  428. - end
  429. - pp = exp(pp - maximum(pp))
  430. - pp /= sum(pp)
  431. -
  432. - # sample from pp
  433. - kk = sampleindex(pp)
  434. - zz[ii] = kk
  435. - nn[kk] += 1
  436. - additem(bayesian_components[kk], xx[ii])
  437. -
  438. - for kk = 1:K_truncation-1
  439. - gamma1 = 1 + nn[kk]
  440. - gamma2 = dpm.aa + sum(nn[kk+1:end])
  441. - beta_VV[kk] = rand(Distributions.Beta(gamma1, gamma2))
  442. - end
  443. - beta_VV[end] = 1.0
  444. - π = ones(Float64, K_truncation)
  445. - π[2:end] = 1 - beta_VV[1:K_truncation-1]
  446. - π = log(beta_VV) + log(cumprod(π))
  447. -
  448. - # resampling concentration parameter based on Escobar and West 1995
  449. - for internal_iters = 1:n_internals
  450. - eta = rand(Distributions.Beta(dpm.aa+1, NN))
  451. - rr = (dpm.a1+dpm.K-1) / (NN*(dpm.a2-log(NN)))
  452. - pi_eta = rr / (1+rr)
  453. -
  454. - if rand() < pi_eta
  455. - dpm.aa = rand(Distributions.Gamma(dpm.a1+dpm.K, 1/(dpm.a2-log(eta))))
  456. - else
  457. - dpm.aa = rand(Distributions.Gamma(dpm.a1+dpm.K-1, 1/(dpm.a2-log(eta))))
  458. - end
  459. - end
  460. - end
  461. -
  462. - nn_string = nn2string(nn)
  463. - if !haskey(states, nn_string)
  464. - n_states += 1
  465. - states[nn_string] = n_states
  466. - end
  467. - dpm.K_hist[iteration] = states[nn_string]
  468. - dpm.K_zz_dict[states[nn_string]] = deepcopy(zz)
  469. - elapsed_time = toq()
  470. - end
  471. - return states
  472. - end
  473. -
  474. -
  475. - function posterior{T1, T2}(dpm::DPM{T1}, xx::Vector{T2}, K::Int64, K_truncation::Int64=0)
  476. 0 n_components = 0
  477. 0 if K_truncation == 0
  478. 0 n_components = K
  479. - else
  480. 0 n_components = K_truncation
  481. - end
  482. -
  483. 0 bayesian_components = [deepcopy(dpm.bayesian_component) for kk=1:n_components]
  484. 0 zz = dpm.K_zz_dict[K]
  485. -
  486. 0 NN = length(xx)
  487. 0 nn = zeros(Int64, n_components)
  488. -
  489. 0 for ii = 1:NN
  490. 0 kk = zz[ii]
  491. 0 additem(bayesian_components[kk], xx[ii])
  492. 0 nn[kk] += 1
  493. - end
  494. -
  495. 0 return([posterior(bayesian_components[kk]) for kk=1:n_components], nn)
  496. - end
  497. -
  498.  
  499. 2 NN = length(xx) # number of data points
  500.  
  501. 191688 NN = length(xx) # number of data points
  502.  
  503. 352 @inbounds for ii = 1:NN
  504.  
  505. 28000000 @inbounds for ii = 1:NN
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement