_sampling.rst 14.3 KB
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
.. _sampling:

Sampling from the posterior distribution
----------------------------------------

When analyzing empirical networks, one should be open to the possibility
that there will be more than one fit of the SBM with similar posterior
probabilities. In such situations, one should instead `sample`
partitions from the posterior distribution, instead of simply finding
its maximum. One can then compute quantities that are averaged over the
different model fits, weighted according to their posterior
probabilities.

Full support for model averaging is implemented in ``graph-tool`` via an
efficient `Markov chain Monte Carlo (MCMC)
<https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo>`_ algorithm
[peixoto-efficient-2014]_. It works by attempting to move nodes into
different groups with specific probabilities, and `accepting or
rejecting
<https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm>`_
such moves so that, after a sufficiently long time, the partitions will
be observed with the desired posterior probability. The algorithm is
designed so that its run-time (i.e. each sweep of the MCMC) is linear on
the number of edges in the network, and independent on the number of
groups being used in the model, and hence is suitable for use on very
large networks.

In order to perform such moves, one needs again to operate with
:class:`~graph_tool.inference.blockmodel.BlockState` or
Tiago Peixoto's avatar
Tiago Peixoto committed
30
31
32
33
34
35
36
37
38
39
40
41
:class:`~graph_tool.inference.nested_blockmodel.NestedBlockState`
instances, and calling either their
:meth:`~graph_tool.inference.blockmodel.BlockState.mcmc_sweep` or
:meth:`~graph_tool.inference.blockmodel.BlockState.multiflip_mcmc_sweep`
methods. The former implements a simpler MCMC where a single node is
moved at a time, where the latter is a more efficient version that
performs merges and splits [peixoto-merge-split-2020]_, which should be
in general preferred.

For example, the following will perform 1000 sweeps of the algorithm
with the network of characters in the novel Les Misérables, starting
from a random partition into 20 groups
42
43
44
45
46

.. testcode:: model-averaging

   g = gt.collection.data["lesmis"]

Tiago Peixoto's avatar
Tiago Peixoto committed
47
48
49
50
   state = gt.BlockState(g, B=1)   # This automatically initializes the state with a partition
                                   # into one group. The user could also pass a higher number
                                   # to start with a random partition of a given size, or pass a
                                   # specific initial partition using the 'b' parameter.
51
52
53

   # Now we run 1,000 sweeps of the MCMC. Note that the number of groups
   # is allowed to change, so it will eventually move from the initial
Tiago Peixoto's avatar
Tiago Peixoto committed
54
   # value of B=1 to whatever is most appropriate for the data.
55

Tiago Peixoto's avatar
Tiago Peixoto committed
56
   dS, nattempts, nmoves = state.multiflip_mcmc_sweep(niter=1000)
57
58
59
60
61
62

   print("Change in description length:", dS)
   print("Number of accepted vertex moves:", nmoves)

.. testoutput:: model-averaging

Tiago Peixoto's avatar
Tiago Peixoto committed
63
64
   Change in description length: -71.707346...
   Number of accepted vertex moves: 44055
65

Tiago Peixoto's avatar
Tiago Peixoto committed
66
67
Although the above is sufficient to implement sampling from the
posterior, there is a convenience function called
68
69
70
71
72
73
74
75
76
77
:func:`~graph_tool.inference.mcmc.mcmc_equilibrate` that is intend to
simplify the detection of equilibration, by keeping track of the maximum
and minimum values of description length encountered and how many sweeps
have been made without a "record breaking" event. For example,

.. testcode:: model-averaging

   # We will accept equilibration if 10 sweeps are completed without a
   # record breaking event, 2 consecutive times.

Tiago Peixoto's avatar
Tiago Peixoto committed
78
   gt.mcmc_equilibrate(state, wait=10, nbreaks=2, mcmc_args=dict(niter=10))
79
80
81
82
83

Note that the value of ``wait`` above was made purposefully low so that
the output would not be overly long. The most appropriate value requires
experimentation, but a typically good value is ``wait=1000``.

84
85
The function :func:`~graph_tool.inference.mcmc.mcmc_equilibrate` accepts
a ``callback`` argument that takes an optional function to be invoked
86
after each call to
87
88
89
90
91
92
93
:meth:`~graph_tool.inference.blockmodel.BlockState.multiflip_mcmc_sweep`. This
function should accept a single parameter which will contain the actual
:class:`~graph_tool.inference.blockmodel.BlockState` instance. We will
use this in the example below to collect the posterior vertex marginals
(via :class:`~graph_tool.inference.partition_modes.PartitionModeState`,
which disambiguates group labels [peixoto-revealing-2020]_), i.e. the
posterior probability that a node belongs to a given group:
94
95
96
97
98
99

.. testcode:: model-averaging

   # We will first equilibrate the Markov chain
   gt.mcmc_equilibrate(state, wait=1000, mcmc_args=dict(niter=10))

100
   bs = [] # collect some partitions
101

102
103
104
   def collect_partitions(s):
      global bs
      bs.append(s.b.a.copy())
105

106
107
   # Now we collect partitions for exactly 100,000 sweeps, at intervals
   # of 10 sweeps:
108
   gt.mcmc_equilibrate(state, force_niter=10000, mcmc_args=dict(niter=10),
109
                       callback=collect_partitions)
110

111
112
113
114
   # Disambiguate partitions and obtain marginals
   pmode = gt.PartitionModeState(bs, converge=True)
   pv = pmode.get_marginal(g)
                       
115
116
117
   # Now the node marginals are stored in property map pv. We can
   # visualize them as pie charts on the nodes:
   state.draw(pos=g.vp.pos, vertex_shape="pie", vertex_pie_fractions=pv,
118
              output="lesmis-sbm-marginals.svg")
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

.. figure:: lesmis-sbm-marginals.*
   :align: center
   :width: 450px

   Marginal probabilities of group memberships of the network of
   characters in the novel Les Misérables, according to the
   degree-corrected SBM. The `pie fractions
   <https://en.wikipedia.org/wiki/Pie_chart>`_ on the nodes correspond
   to the probability of being in group associated with the respective
   color.

We can also obtain a marginal probability on the number of groups
itself, as follows.

.. testcode:: model-averaging

   h = np.zeros(g.num_vertices() + 1)

   def collect_num_groups(s):
       B = s.get_nonempty_B()
       h[B] += 1

142
143
   # Now we collect partitions for exactly 100,000 sweeps, at intervals
   # of 10 sweeps:
144
145
146
147
148
149
150
151
152
153
154
155
   gt.mcmc_equilibrate(state, force_niter=10000, mcmc_args=dict(niter=10),
                       callback=collect_num_groups)

.. testcode:: model-averaging
   :hide:

   figure()
   Bs = np.arange(len(h))
   idx = h > 0
   bar(Bs[idx], h[idx] / h.sum(), width=1, color="#ccb974")
   gca().set_xticks([6,7,8,9])
   xlabel("$B$")
156
   ylabel(r"$P(B|\boldsymbol A)$")
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
   savefig("lesmis-B-posterior.svg")

.. figure:: lesmis-B-posterior.*
   :align: center

   Marginal posterior probability of the number of nonempty groups for
   the network of characters in the novel Les Misérables, according to
   the degree-corrected SBM.


Hierarchical partitions
+++++++++++++++++++++++

We can also perform model averaging using the nested SBM, which will
give us a distribution over hierarchies. The whole procedure is fairly
analogous, but now we make use of
:class:`~graph_tool.inference.nested_blockmodel.NestedBlockState` instances.

Here we perform the sampling of hierarchical partitions using the same
network as above.

.. testcode:: nested-model-averaging

   g = gt.collection.data["lesmis"]

Tiago Peixoto's avatar
Tiago Peixoto committed
182
183
   state = gt.NestedBlockState(g)   # By default this creates a state with an initial single-group
                                    # hierarchy of depth ceil(log2(g.num_vertices()).
184
185
186

   # Now we run 1000 sweeps of the MCMC

Tiago Peixoto's avatar
Tiago Peixoto committed
187
188
189
190
191
   dS, nmoves = 0, 0
   for i in range(100):
       ret = state.multiflip_mcmc_sweep(niter=10)
       dS += ret[0]
       nmoves += ret[1]
192
193
194
195
196
197

   print("Change in description length:", dS)
   print("Number of accepted vertex moves:", nmoves)

.. testoutput:: nested-model-averaging

Tiago Peixoto's avatar
Tiago Peixoto committed
198
199
200
201
   Change in description length: -73.716766...
   Number of accepted vertex moves: 366160

.. warning::
202

Tiago Peixoto's avatar
Tiago Peixoto committed
203
204
205
206
207
208
209
210
211
212
213
   When using
   :class:`~graph_tool.inference.nested_blockmodel.NestedBlockState`, a
   single call to
   :meth:`~graph_tool.inference.nested_blockmodel.NestedBlockState.multiflip_mcmc_sweep`
   or
   :meth:`~graph_tool.inference.nested_blockmodel.NestedBlockState.mcmc_sweep`
   performs ``niter`` sweeps at each hierarchical level once. This means
   that in order for the chain to equilibrate, we need to call these
   functions several times, i.e. it is not enough to call it once with a
   large value of ``niter``.
   
214
215
216
217
218
219
220
221
222
Similarly to the the non-nested case, we can use
:func:`~graph_tool.inference.mcmc.mcmc_equilibrate` to do most of the boring
work, and we can now obtain vertex marginals on all hierarchical levels:

.. testcode:: nested-model-averaging

   # We will first equilibrate the Markov chain
   gt.mcmc_equilibrate(state, wait=1000, mcmc_args=dict(niter=10))

223
224
   # collect nested partitions
   bs = []
225

226
227
228
   def collect_partitions(s):
      global bs
      bs.append(s.get_bs())
229
230
231

   # Now we collect the marginals for exactly 100,000 sweeps
   gt.mcmc_equilibrate(state, force_niter=10000, mcmc_args=dict(niter=10),
232
233
234
235
236
237
238
239
240
241
                       callback=collect_partitions)

   # Disambiguate partitions and obtain marginals
   pmode = gt.PartitionModeState(bs, nested=True, converge=True)
   pv = pmode.get_marginal(g)

   # Get consensus estimate
   bs = pmode.get_max_nested()

   state = state.copy(bs=bs)
242

243
244
245
   # We can visualize the marginals as pie charts on the nodes:
   state.draw(vertex_shape="pie", vertex_pie_fractions=pv,
              output="lesmis-nested-sbm-marginals.svg")
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283

.. figure:: lesmis-nested-sbm-marginals.*
   :align: center
   :width: 450px

   Marginal probabilities of group memberships of the network of
   characters in the novel Les Misérables, according to the nested
   degree-corrected SBM. The pie fractions on the nodes correspond to
   the probability of being in group associated with the respective
   color.

We can also obtain a marginal probability of the number of groups
itself, as follows.

.. testcode:: nested-model-averaging

   h = [np.zeros(g.num_vertices() + 1) for s in state.get_levels()]

   def collect_num_groups(s):
       for l, sl in enumerate(s.get_levels()):
          B = sl.get_nonempty_B()
          h[l][B] += 1

   # Now we collect the marginal distribution for exactly 100,000 sweeps
   gt.mcmc_equilibrate(state, force_niter=10000, mcmc_args=dict(niter=10),
                       callback=collect_num_groups)

.. testcode:: nested-model-averaging
   :hide:

   figure()
   f, ax = plt.subplots(1, 5, figsize=(10, 3))
   for i, h_ in enumerate(h[:5]):
       Bs = np.arange(len(h_))
       idx = h_ > 0
       ax[i].bar(Bs[idx], h_[idx] / h_.sum(), width=1, color="#ccb974")
       ax[i].set_xticks(Bs[idx])
       ax[i].set_xlabel("$B_{%d}$" % i)
284
       ax[i].set_ylabel(r"$P(B_{%d}|\boldsymbol A)$" % i)
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
       locator = MaxNLocator(prune='both', nbins=5)
       ax[i].yaxis.set_major_locator(locator)
   tight_layout()
   savefig("lesmis-nested-B-posterior.svg")

.. figure:: lesmis-nested-B-posterior.*
   :align: center

   Marginal posterior probability of the number of nonempty groups
   :math:`B_l` at each hierarchy level :math:`l` for the network of
   characters in the novel Les Misérables, according to the nested
   degree-corrected SBM.

Below we obtain some hierarchical partitions sampled from the posterior
distribution.

.. testcode:: nested-model-averaging

   for i in range(10):
Tiago Peixoto's avatar
Tiago Peixoto committed
304
305
       for j in range(100):
           state.multiflip_mcmc_sweep(niter=10)
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
       state.draw(output="lesmis-partition-sample-%i.svg" % i, empty_branches=False)

.. image:: lesmis-partition-sample-0.svg
   :width: 200px
.. image:: lesmis-partition-sample-1.svg
   :width: 200px
.. image:: lesmis-partition-sample-2.svg
   :width: 200px
.. image:: lesmis-partition-sample-3.svg
   :width: 200px
.. image:: lesmis-partition-sample-4.svg
   :width: 200px
.. image:: lesmis-partition-sample-5.svg
   :width: 200px
.. image:: lesmis-partition-sample-6.svg
   :width: 200px
.. image:: lesmis-partition-sample-7.svg
   :width: 200px
.. image:: lesmis-partition-sample-8.svg
   :width: 200px
.. image:: lesmis-partition-sample-9.svg
   :width: 200px
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403

Characterizing the posterior distribution
+++++++++++++++++++++++++++++++++++++++++

The posterior distribution of partitions can have an elaborate
structure, containing multiple possible explanations for the data. In
order to summarize it, we can infer the modes of the distribution using
:class:`~graph_tool.inference.partition_modes.ModeClusterState`, as
described in [peixoto-revealing-2020]_. This amounts to identifying
clusters of partitions that are very similar to each other, but
sufficiently different from those that belong to other
clusters. Collective, such "modes" represent the different stories that
the data is telling us through the model. Here is an example using again
the Les Misérables network:

.. testcode:: partition-modes

   g = gt.collection.data["lesmis"]

   state = gt.NestedBlockState(g)

   # Equilibration
   gt.mcmc_equilibrate(state, force_niter=1000, mcmc_args=dict(niter=10))

   bs = []
   
   def collect_partitions(s):
      global bs
      bs.append(s.get_bs())

   # We will collect only partitions 1000 partitions. For more accurate
   # results, this number should be increased.
   gt.mcmc_equilibrate(state, force_niter=1000, mcmc_args=dict(niter=10),
                       callback=collect_partitions)

   # Infer partition modes
   pmode = gt.ModeClusterState(bs, nested=True)

   # Minimize the mode state itself
   gt.mcmc_equilibrate(pmode, wait=1, mcmc_args=dict(niter=1, beta=np.inf))

   # Get inferred modes
   modes = pmode.get_modes()

   for i, mode in enumerate(modes):
       b = mode.get_max_nested()    # mode's maximum
       pv = mode.get_marginal(g)    # mode's marginal distribution

       print(f"Mode {i} with size {mode.get_M()/len(bs)}")
       state = state.copy(bs=b)
       state.draw(vertex_shape="pie", vertex_pie_fractions=pv,
                  output="lesmis-partition-mode-%i.svg" % i)

Running the above code gives us the relative size of each mode,
corresponding to their collective posterior probability.

.. testoutput:: partition-modes

    Mode 0 with size 0.389389...
    Mode 1 with size 0.352352...
    Mode 2 with size 0.129129...
    Mode 3 with size 0.117117...
    Mode 4 with size 0.012012...
                  
Below are the marginal node distributions representing the partitions that belong to each inferred mode:
       
.. image:: lesmis-partition-mode-0.svg
   :width: 200px
.. image:: lesmis-partition-mode-1.svg
   :width: 200px
.. image:: lesmis-partition-mode-2.svg
   :width: 200px
.. image:: lesmis-partition-mode-3.svg
   :width: 200px
.. image:: lesmis-partition-mode-4.svg
   :width: 200px