Skip to content

Dirichlet Process simple example #125

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
gcgibson wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Dirichlet Process simple example #125

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
135 changes: 135 additions & 0 deletions doc/examples/dirichlet_process.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,135 @@
#!/Applications/Julia-0.6.app/Contents/Resources/julia/bin/julia



## Define a new univariate Distribution type for Mamba.
## The definition must be placed within an unevaluated quote block.

## Data
data = Dict{Symbol, Any}(
:y => randn(100) + 5
)

## Maximum number of active components
number_of_components = 40


function stick_breaking(beta)
cum_prod = 1 - beta[1]
weights = ones(number_of_components)
weights[1] = beta[1]
for i = 2:number_of_components
weights[i] = beta[i]*cum_prod
cum_prod = cum_prod*(1-beta[i])
end
weights
end




@everywhere extensions = quote

## Load needed packages and import methods to be extended
using Distributions
import Distributions: minimum, maximum, logpdf

## Type declaration
type NewUnivarDist <: ContinuousUnivariateDistribution
mu::Array{Float64,1}
sigma::Float64
w::Array{Float64,1}
end

## The following method functions must be implemented

## Minimum and maximum support values
minimum(d::NewUnivarDist) = -Inf
maximum(d::NewUnivarDist) = Inf

## Normalized or unnormalized log-density value
function logpdf(d::NewUnivarDist, x::Real)
logdens = 0
for i in 1:40
logdens = logdens + d.w[i]*(1.0/sqrt((2*pi*d.sigma)))*exp( - 0.5 * ((x - d.mu[i]) / d.sigma)^2)
end
log(logdens)
end


end


## Add the extensions
using Mamba
@everywhere eval(extensions)
## Implement a Mamba model using the new distribution
model = Model(

y = Stochastic(1,
(beta,mu,tau) ->
begin
weights = stick_breaking(beta)
UnivariateDistribution[
NewUnivarDist(mu, tau,weights) for i in 1:100
]
end,
false
),
ytilde = Logical(1,
(beta,mu,tau) ->
begin
weights = stick_breaking(beta)
idx = wsample(1:number_of_components,weights )
rands = rand(Normal(mu[idx], tau)
,1)
rands
end
),
beta = Stochastic(1,
(alpha) ->
Beta(1, alpha)
),

mu = Stochastic(1,
(tau) ->
begin
Normal(0,tau)
end
),

tau = Stochastic(
() -> InverseGamma(1,1)
),
alpha = Stochastic(
() -> Uniform(0,1)
)
)


## Sampling Scheme
## Initial Values
inits = [
Dict{Symbol, Any}(
:y => data[:y],
:beta =>rand(number_of_components)*1,
:mu => rand(number_of_components),
:tau => 1,
:alpha => .1,
) for i in 1:3
]


## Sampling Scheme

scheme = [Slice([:beta,:mu,:tau,:alpha],3.0)]

setsamplers!(model, scheme)


## MCMC Simulations
sim = mcmc(model, data, inits, 10000, burnin=2500, thin=2, chains=3)
describe(sim)



209 changes: 209 additions & 0 deletions doc/examples/dirichlet_process.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,209 @@
.. index:: Examples; Dirichlet Process: Simulaion Study

@bdeonovic bdeonovic Aug 28, 2017

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fix typos of "simulation" and check rest of the document as well.


.. _example-DP:

Dirichlet Process: Simulaion Study
------------------------------------------------------

An example of a dirichlet process using the beta stick-breaking construction. Right now the example is a simple simulation study where we try and recover the parameters of a single normal distribution.

Model
^^^^^


.. math::
\y \sim \sum_{i=1}^\infty \w_i * N(\mu_i,tau)
\w_i = \beta_i \prod_{j=1}^i(1-\beta_j)
\beta_i &\sim \text{Beta}(1,alpha)
\mu_i &\sim \text{Normal}(0,\sigma^2)
\sigma^2 &\sim \text{InverseGamma}(1,1)\\
\alpha &\sim \text{Uniform}(0,1),



Analysis Program
^^^^^^^^^^^^^^^^

.. literalinclude:: dirichlet_process.jl
:language: julia

Results
^^^^^^^

.. code-block:: julia
Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2,3
Samples per chain = 3750

Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
tau 1.3409189107 0.09878433 0.0009313476 0.0077409625 162.849404
mu[1] 2.9507581422 2.74822898 0.0259105513 0.2597197646 111.968491
mu[2] 1.5383252240 1.59387997 0.0150272445 0.1507064742 111.852941
mu[3] 1.4356038154 0.93831604 0.0088465285 0.0886178529 112.112964
mu[4] 1.3054031597 0.39711734 0.0037440582 0.0368472602 116.152164
mu[5] 1.8843424785 2.23844766 0.0211042870 0.2111706883 112.363852
mu[6] 0.0134000782 0.65001929 0.0061284407 0.0613912271 112.108841
mu[7] 0.5272592879 0.80933985 0.0076305292 0.0764694066 112.017636
mu[8] 0.3309302588 0.71639841 0.0067542690 0.0675074715 112.617413
mu[9] 0.3406224923 0.82915095 0.0078173101 0.0783005727 112.134003
mu[10] -0.2856097066 0.48947873 0.0046148497 0.0455508607 115.471411
mu[11] -0.3974026638 1.31933175 0.0124387790 0.1249654487 111.462331
mu[12] 0.3524125280 0.98674188 0.0093030917 0.0932934619 111.867699
mu[13] 0.6800066800 0.43417612 0.0040934517 0.0405786288 114.481978
mu[14] 1.3087578174 0.40174375 0.0037876764 0.0369926573 117.941651
mu[15] 0.6190732220 1.06769592 0.0100663336 0.1003609395 113.178967
mu[16] 1.2617279963 0.93387787 0.0088046850 0.0882955983 111.867017
mu[17] 0.5585084538 1.04807087 0.0098813069 0.0986727064 112.820291
mu[18] 1.0793786870 0.82559848 0.0077838172 0.0779408862 112.203678
mu[19] -1.0364919319 1.20745372 0.0113839828 0.1136159762 112.943744
mu[20] 0.9353933268 0.26817106 0.0025283410 0.0240000199 124.853470
mu[21] 0.2352632813 1.21570703 0.0114617958 0.1152908155 111.190693
mu[22] 0.9728287689 0.69657195 0.0065673433 0.0649557797 114.999609
mu[23] -0.0389055206 0.79911028 0.0075340840 0.0751938841 112.940160
mu[24] -0.3929568143 1.09929455 0.0103642484 0.1037916252 112.176923
mu[25] 0.7916696179 0.75688650 0.0071359944 0.0714485757 112.221061
mu[26] 0.6316871871 0.61808325 0.0058273448 0.0573451898 116.171598
mu[27] 0.9210721135 0.60139803 0.0056700350 0.0567541181 112.286940
mu[28] -0.1707555961 1.11959255 0.0105556197 0.1056255114 112.352398
mu[29] 0.5939725280 0.45884565 0.0043260383 0.0420944705 118.818247
mu[30] 0.3610823839 0.58278815 0.0054945794 0.0545514756 114.132260
mu[31] -0.0034088155 0.57585081 0.0054291735 0.0539313212 114.008663
mu[32] 0.5597342110 0.93409245 0.0088067081 0.0877286420 113.369678
mu[33] 0.3532020355 1.02386796 0.0096531197 0.0964785946 112.622709
mu[34] 0.7120790620 0.95197558 0.0089753118 0.0902459246 111.274697
mu[35] 0.2170323486 0.97318929 0.0091753166 0.0918117041 112.356589
mu[36] 0.0361374512 0.71235889 0.0067161841 0.0670234976 112.964901
mu[37] 0.5665202983 0.60091594 0.0056654898 0.0553486928 117.872560
mu[38] 0.1005468489 0.64079609 0.0060414834 0.0587461394 118.981938
mu[39] 0.8293418044 0.69574775 0.0065595727 0.0654664486 112.944752
mu[40] 0.0516182381 1.02251515 0.0096403653 0.0961265247 113.149600
alpha 0.5996175545 0.17032649 0.0016058536 0.0153680322 122.836626
beta[1] 0.6485027561 0.45696991 0.0043083536 0.0431494988 112.156305
beta[2] 0.3819325195 0.38110623 0.0035931040 0.0358857281 112.784280
beta[3] 0.5448986051 0.40975705 0.0038632265 0.0385678290 112.876226
beta[4] 0.4625597345 0.36942924 0.0034830123 0.0345141907 114.569005
beta[5] 0.8273723394 0.19664935 0.0018540278 0.0180990510 118.052018
beta[6] 0.6819997468 0.27308705 0.0025746894 0.0252242596 117.210186
beta[7] 0.6815609901 0.24248609 0.0022861808 0.0222511296 118.759821
beta[8] 0.6058963239 0.34258825 0.0032299530 0.0318404985 115.767117
beta[9] 0.8014695333 0.20561826 0.0019385876 0.0183924372 124.981235
beta[10] 0.8068019819 0.19501205 0.0018385912 0.0175398671 123.614745
beta[11] 0.8020523214 0.19029161 0.0017940865 0.0169717762 125.714297
beta[12] 0.7475659216 0.22581892 0.0021290412 0.0203467778 123.176928
beta[13] 0.5759680227 0.30120561 0.0028397938 0.0280358438 115.424729
beta[14] 0.5960398327 0.33577643 0.0031657305 0.0313266495 114.887432
beta[15] 0.5642678629 0.30322934 0.0028588737 0.0279210265 117.945043
beta[16] 0.5213643445 0.29587529 0.0027895389 0.0275516842 115.324366
beta[17] 0.4686632235 0.29314178 0.0027637672 0.0267749938 119.866185
beta[18] 0.8026484576 0.20468072 0.0019297483 0.0180230163 128.973033
beta[19] 0.5759805704 0.33128440 0.0031233793 0.0306338936 116.949285
beta[20] 0.7186536019 0.24857832 0.0023436189 0.0228825955 118.009219
beta[21] 0.7209482274 0.30569537 0.0028821236 0.0285080772 114.985178
beta[22] 0.7830750801 0.20366640 0.0019201853 0.0180991106 126.626418
beta[23] 0.3827724460 0.28938600 0.0027283573 0.0264295556 119.887886
beta[24] 0.4536858291 0.27588549 0.0026010734 0.0252356361 119.516871
beta[25] 0.6038923328 0.33841513 0.0031906084 0.0314900240 115.492439
beta[26] 0.6883597856 0.27555475 0.0025979551 0.0256074193 115.793496
beta[27] 0.5379037321 0.32986998 0.0031100440 0.0303822776 117.881304
beta[28] 0.5981627541 0.32447866 0.0030592142 0.0298287237 118.332200
beta[29] 0.5669150500 0.31079445 0.0029301982 0.0288280351 116.229516
beta[30] 0.5492389198 0.31908958 0.0030084054 0.0295162653 116.869838
beta[31] 0.6875947580 0.31701557 0.0029888514 0.0294560877 115.827331
beta[32] 0.4226177579 0.31049135 0.0029273405 0.0286612739 117.356736
beta[33] 0.7043622405 0.31687618 0.0029875373 0.0292196208 117.606149
beta[34] 0.6447207268 0.32895731 0.0031014393 0.0307792659 114.225370
beta[35] 0.8017960783 0.18609618 0.0017545316 0.0163827900 129.032516
beta[36] 0.6710025627 0.31450626 0.0029651934 0.0293433723 114.878441
beta[37] 0.6004724338 0.31537729 0.0029734056 0.0292087724 116.582734
beta[38] 0.5419600723 0.24576053 0.0023170525 0.0226226363 118.015169
beta[39] 0.6201564297 0.32215754 0.0030373304 0.0298175456 116.732777
beta[40] 0.5148273475 0.31884911 0.0030061383 0.0294777659 116.998770
ytilde[1] 4.8056180090 1.49118492 0.0140590263 0.0231083291 3750.000000

Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
tau 1.16042118422 1.2701693552 1.335733733 1.414383072 1.54058874
mu[1] -1.33182698945 -0.7494224304 4.796129701 4.949127828 5.10299937
mu[2] -0.04871993014 0.2677637896 0.744801643 3.198133975 4.91095265
mu[3] 0.01187530042 0.7032553311 1.206200537 2.385028842 3.02910454
mu[4] 0.47060834370 1.0692001995 1.311619456 1.601211363 1.88862573
mu[5] -1.02205054351 -0.3776356227 1.576130475 4.732470982 4.97238447
mu[6] -1.02740582685 -0.4801561189 -0.019554284 0.500408706 1.18503208
mu[7] -0.87331669565 -0.1692517064 0.755224195 1.097845398 1.90942121
mu[8] -0.79714435039 -0.2539386158 0.335614253 0.890023656 1.63412466
mu[9] -1.12531648350 -0.2792846942 0.391547304 0.994752520 1.72236390
mu[10] -1.02618791009 -0.6621847626 -0.324801710 -0.034875655 0.81445265
mu[11] -2.94644972076 -1.7651287177 -0.058093989 0.619417734 1.53584660
mu[12] -1.65596065247 -0.5625243569 0.706033389 1.031338981 1.60104164
mu[13] -0.20299795689 0.4497656880 0.595315742 0.868967606 1.68293747
mu[14] 0.51217380086 0.9291731765 1.352409376 1.666990551 1.95930339
mu[15] -1.30532669405 -0.7866210752 1.127850246 1.499380750 1.77592532
mu[16] -0.97555236545 0.9063499195 1.465763050 1.905559607 2.68703197
mu[17] -0.80398265988 -0.2908288262 0.148363999 1.801527812 2.32018859
mu[18] -0.12495669368 0.4034646606 1.108662113 1.685774690 2.56801335
mu[19] -2.94845886794 -2.4110893705 -1.084634123 0.162467104 0.86814407
mu[20] 0.47513603459 0.7967962863 0.928243140 1.029432640 1.61916216
mu[21] -2.11585492353 -0.4958998135 0.174147414 0.850931068 2.53080626
mu[22] -0.27090087307 0.5253248392 0.916552693 1.588960187 2.24102759
mu[23] -1.38569900920 -0.6272128525 -0.263769891 0.660453705 1.45641011
mu[24] -2.47500898757 -1.2035093958 -0.633408205 0.817512933 1.43654857
mu[25] -0.32347392367 0.1879206477 0.680439070 1.281831460 2.32490258
mu[26] -0.43384467583 0.1588895598 0.494592611 1.041763374 1.93110394
mu[27] -0.52686318525 0.6983878530 0.832808676 1.437281535 1.85840559
mu[28] -2.36801259849 -0.6868855613 -0.274197820 0.814621912 1.51794163
mu[29] -0.30669744936 0.3416763876 0.554549372 0.896192830 1.44214579
mu[30] -0.79322316019 -0.0780938296 0.336263577 0.767247618 1.51016299
mu[31] -0.95915359104 -0.4497734217 0.011036136 0.423099821 1.00074174
mu[32] -1.27731901119 -0.3166590983 0.891897837 1.262508489 1.82831551
mu[33] -1.49185054300 -0.4501366364 0.707317545 1.166123581 1.65456290
mu[34] -0.33996121446 -0.0068196364 0.339092603 1.206825174 2.90582045
mu[35] -1.46730749443 -0.9920965157 0.640784686 0.966854101 1.32624251
mu[36] -1.18059022411 -0.6048175874 0.126344774 0.676624725 1.18465365
mu[37] -0.76688599902 0.1199904973 0.643165038 1.021848505 1.54376759
mu[38] -0.78194378762 -0.4665750629 0.141273617 0.491593478 1.22264839
mu[39] -0.26384636705 0.3578969449 0.608569322 1.473687928 2.28388565
mu[40] -1.87801752828 -0.8492529286 0.239535376 0.995290414 1.32744807
alpha 0.35773476132 0.4695608010 0.556217267 0.728709785 0.96227410
beta[1] 0.00079437523 0.0140449676 0.984260739 0.997832550 0.99997849
beta[2] 0.00066983577 0.0121149202 0.236813959 0.785859863 0.99908110
beta[3] 0.00066506903 0.0137153032 0.750645209 0.929638750 0.99871662
beta[4] 0.00084541361 0.0133814772 0.519673884 0.790865363 0.99369005
beta[5] 0.34422329999 0.6607202185 0.918521258 0.992551760 0.99975077
beta[6] 0.05723860386 0.4953541936 0.738717802 0.922496836 0.99884857
beta[7] 0.13626820238 0.4880257409 0.742279225 0.879321173 0.99686108
beta[8] 0.03413244317 0.2119677384 0.748879264 0.900505040 0.99956953
beta[9] 0.24763314330 0.6705026171 0.884527041 0.967268193 0.99931901
beta[10] 0.36196640709 0.6877008105 0.880889568 0.961814735 0.99956630
beta[11] 0.39108647489 0.7007951119 0.850000420 0.968725204 0.99987275
beta[12] 0.18835131302 0.6092682467 0.796914009 0.948810494 0.99933958
beta[13] 0.04401429693 0.3394642571 0.625190041 0.849704388 0.99619359
beta[14] 0.04530348842 0.2770910774 0.624937290 0.932065764 0.99864989
beta[15] 0.02546858787 0.2734590170 0.606785081 0.823325853 0.99658484
beta[16] 0.01143276922 0.3232026336 0.562250102 0.765277378 0.99371241
beta[17] 0.02778886147 0.2130092092 0.422787631 0.742362781 0.98479894
beta[18] 0.24808132995 0.7076860617 0.872123087 0.960297538 0.99930554
beta[19] 0.02481856633 0.2828764615 0.622204558 0.915038126 0.99712793
beta[20] 0.15532773437 0.5247423971 0.797133505 0.938406639 0.99907142
beta[21] 0.04159568220 0.5008176279 0.857625918 0.970940408 0.99979031
beta[22] 0.26235886597 0.6808012189 0.832507391 0.952236595 0.99936737
beta[23] 0.01314179563 0.1494234163 0.319678783 0.523364813 0.99853921
beta[24] 0.04193052816 0.2061918878 0.442045969 0.675879340 0.96535570
beta[25] 0.03716832994 0.2795827559 0.731415088 0.908349165 0.99977847
beta[26] 0.08405175969 0.5331120740 0.749977201 0.925260257 0.99862839
beta[27] 0.02267531531 0.2129242163 0.509712881 0.855671926 0.99944064
beta[28] 0.00838615212 0.3049402746 0.680379164 0.888892769 0.99761993
beta[29] 0.06568169227 0.2274457374 0.663187284 0.834694160 0.99089437
beta[30] 0.05884160295 0.2646533128 0.522027238 0.890599837 0.99871650
beta[31] 0.01357712731 0.4868519798 0.825675994 0.947383005 0.99958704
beta[32] 0.02309517030 0.1497337317 0.316474852 0.718867895 0.99648922
beta[33] 0.03331411219 0.5150319331 0.828142891 0.977327108 0.99995811
beta[34] 0.05033238559 0.3875906231 0.717680826 0.950242485 0.99983431
beta[35] 0.31479129593 0.7294004834 0.845587879 0.951520172 0.99959759
beta[36] 0.03791703780 0.3654268365 0.802488909 0.954852016 0.99920767
beta[37] 0.02887787368 0.3382716381 0.611321031 0.908120190 0.99796258
beta[38] 0.06341878500 0.3831618173 0.535211955 0.729341173 0.98872744
beta[39] 0.02108327846 0.3070900494 0.690057656 0.915977866 0.99980189
beta[40] 0.03189814221 0.2416422111 0.436136126 0.872523610 0.99845551
ytilde[1] 1.70354908344 3.8979349274 4.850851814 5.771041129 7.55571318