Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/docs/usage.html @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 <HEAD>
wolffd@0 2 <TITLE>How to use the Bayes Net Toolbox</TITLE>
wolffd@0 3 </HEAD>
wolffd@0 4
wolffd@0 5 <BODY BGCOLOR="#FFFFFF">
wolffd@0 6 <!-- white background is better for the pictures and equations -->
wolffd@0 7
wolffd@0 8 <h1>How to use the Bayes Net Toolbox</h1>
wolffd@0 9
wolffd@0 10 This documentation was last updated on 29 October 2007.
wolffd@0 11 <br>
wolffd@0 12 Click
wolffd@0 13 <a href="http://bnt.insa-rouen.fr/">here</a>
wolffd@0 14 for a French version of this documentation (last updated in 2005).
wolffd@0 15 <p>
wolffd@0 16
wolffd@0 17 <ul>
wolffd@0 18 <li> <a href="install.html">Installation</a>
wolffd@0 19
wolffd@0 20 <li> <a href="#basics">Creating your first Bayes net</a>
wolffd@0 21 <ul>
wolffd@0 22 <li> <a href="#basics">Creating a model by hand</a>
wolffd@0 23 <li> <a href="#file">Loading a model from a file</a>
wolffd@0 24 <li> <a href="#GUI">Creating a model using a GUI</a>
wolffd@0 25 <li> <a href="graphviz.html">Graph visualization</a>
wolffd@0 26 </ul>
wolffd@0 27
wolffd@0 28 <li> <a href="#inference">Inference</a>
wolffd@0 29 <ul>
wolffd@0 30 <li> <a href="#marginal">Computing marginal distributions</a>
wolffd@0 31 <li> <a href="#joint">Computing joint distributions</a>
wolffd@0 32 <li> <a href="#soft">Soft/virtual evidence</a>
wolffd@0 33 <li> <a href="#mpe">Most probable explanation</a>
wolffd@0 34 </ul>
wolffd@0 35
wolffd@0 36 <li> <a href="#cpd">Conditional Probability Distributions</a>
wolffd@0 37 <ul>
wolffd@0 38 <li> <a href="#tabular">Tabular (multinomial) nodes</a>
wolffd@0 39 <li> <a href="#noisyor">Noisy-or nodes</a>
wolffd@0 40 <li> <a href="#deterministic">Other (noisy) deterministic nodes</a>
wolffd@0 41 <li> <a href="#softmax">Softmax (multinomial logit) nodes</a>
wolffd@0 42 <li> <a href="#mlp">Neural network nodes</a>
wolffd@0 43 <li> <a href="#root">Root nodes</a>
wolffd@0 44 <li> <a href="#gaussian">Gaussian nodes</a>
wolffd@0 45 <li> <a href="#glm">Generalized linear model nodes</a>
wolffd@0 46 <li> <a href="#dtree">Classification/regression tree nodes</a>
wolffd@0 47 <li> <a href="#nongauss">Other continuous distributions</a>
wolffd@0 48 <li> <a href="#cpd_summary">Summary of CPD types</a>
wolffd@0 49 </ul>
wolffd@0 50
wolffd@0 51 <li> <a href="#examples">Example models</a>
wolffd@0 52 <ul>
wolffd@0 53 <li> <a
wolffd@0 54 href="http://www.media.mit.edu/wearables/mithril/BNT/mixtureBNT.txt">
wolffd@0 55 Gaussian mixture models</a>
wolffd@0 56 <li> <a href="#pca">PCA, ICA, and all that</a>
wolffd@0 57 <li> <a href="#mixep">Mixtures of experts</a>
wolffd@0 58 <li> <a href="#hme">Hierarchical mixtures of experts</a>
wolffd@0 59 <li> <a href="#qmr">QMR</a>
wolffd@0 60 <li> <a href="#cg_model">Conditional Gaussian models</a>
wolffd@0 61 <li> <a href="#hybrid">Other hybrid models</a>
wolffd@0 62 </ul>
wolffd@0 63
wolffd@0 64 <li> <a href="#param_learning">Parameter learning</a>
wolffd@0 65 <ul>
wolffd@0 66 <li> <a href="#load_data">Loading data from a file</a>
wolffd@0 67 <li> <a href="#mle_complete">Maximum likelihood parameter estimation from complete data</a>
wolffd@0 68 <li> <a href="#prior">Parameter priors</a>
wolffd@0 69 <li> <a href="#bayes_learn">(Sequential) Bayesian parameter updating from complete data</a>
wolffd@0 70 <li> <a href="#em">Maximum likelihood parameter estimation with missing values (EM)</a>
wolffd@0 71 <li> <a href="#tying">Parameter tying</a>
wolffd@0 72 </ul>
wolffd@0 73
wolffd@0 74 <li> <a href="#structure_learning">Structure learning</a>
wolffd@0 75 <ul>
wolffd@0 76 <li> <a href="#enumerate">Exhaustive search</a>
wolffd@0 77 <li> <a href="#K2">K2</a>
wolffd@0 78 <li> <a href="#hill_climb">Hill-climbing</a>
wolffd@0 79 <li> <a href="#mcmc">MCMC</a>
wolffd@0 80 <li> <a href="#active">Active learning</a>
wolffd@0 81 <li> <a href="#struct_em">Structural EM</a>
wolffd@0 82 <li> <a href="#graphdraw">Visualizing the learned graph structure</a>
wolffd@0 83 <li> <a href="#constraint">Constraint-based methods</a>
wolffd@0 84 </ul>
wolffd@0 85
wolffd@0 86
wolffd@0 87 <li> <a href="#engines">Inference engines</a>
wolffd@0 88 <ul>
wolffd@0 89 <li> <a href="#jtree">Junction tree</a>
wolffd@0 90 <li> <a href="#varelim">Variable elimination</a>
wolffd@0 91 <li> <a href="#global">Global inference methods</a>
wolffd@0 92 <li> <a href="#quickscore">Quickscore</a>
wolffd@0 93 <li> <a href="#belprop">Belief propagation</a>
wolffd@0 94 <li> <a href="#sampling">Sampling (Monte Carlo)</a>
wolffd@0 95 <li> <a href="#engine_summary">Summary of inference engines</a>
wolffd@0 96 </ul>
wolffd@0 97
wolffd@0 98
wolffd@0 99 <li> <a href="#influence">Influence diagrams/ decision making</a>
wolffd@0 100
wolffd@0 101
wolffd@0 102 <li> <a href="usage_dbn.html">DBNs, HMMs, Kalman filters and all that</a>
wolffd@0 103 </ul>
wolffd@0 104
wolffd@0 105 </ul>
wolffd@0 106
wolffd@0 107
wolffd@0 108
wolffd@0 109 <h1><a name="basics">Creating your first Bayes net</h1>
wolffd@0 110
wolffd@0 111 To define a Bayes net, you must specify the graph structure and then
wolffd@0 112 the parameters. We look at each in turn, using a simple example
wolffd@0 113 (adapted from Russell and
wolffd@0 114 Norvig, "Artificial Intelligence: a Modern Approach", Prentice Hall,
wolffd@0 115 1995, p454).
wolffd@0 116
wolffd@0 117
wolffd@0 118 <h2>Graph structure</h2>
wolffd@0 119
wolffd@0 120
wolffd@0 121 Consider the following network.
wolffd@0 122
wolffd@0 123 <p>
wolffd@0 124 <center>
wolffd@0 125 <IMG SRC="Figures/sprinkler.gif">
wolffd@0 126 </center>
wolffd@0 127 <p>
wolffd@0 128
wolffd@0 129 <P>
wolffd@0 130 To specify this directed acyclic graph (dag), we create an adjacency matrix:
wolffd@0 131 <PRE>
wolffd@0 132 N = 4;
wolffd@0 133 dag = zeros(N,N);
wolffd@0 134 C = 1; S = 2; R = 3; W = 4;
wolffd@0 135 dag(C,[R S]) = 1;
wolffd@0 136 dag(R,W) = 1;
wolffd@0 137 dag(S,W)=1;
wolffd@0 138 </PRE>
wolffd@0 139 <P>
wolffd@0 140 We have numbered the nodes as follows:
wolffd@0 141 Cloudy = 1, Sprinkler = 2, Rain = 3, WetGrass = 4.
wolffd@0 142 <b>The nodes must always be numbered in topological order, i.e.,
wolffd@0 143 ancestors before descendants.</b>
wolffd@0 144 For a more complicated graph, this is a little inconvenient: we will
wolffd@0 145 see how to get around this <a href="usage_dbn.html#bat">below</a>.
wolffd@0 146 <p>
wolffd@0 147 In Matlab 6, you can use logical arrays instead of double arrays,
wolffd@0 148 which are 4 times smaller:
wolffd@0 149 <pre>
wolffd@0 150 dag = false(N,N);
wolffd@0 151 dag(C,[R S]) = true;
wolffd@0 152 ...
wolffd@0 153 </pre>
wolffd@0 154 However, <b>some graph functions (eg acyclic) do not work on
wolffd@0 155 logical arrays</b>!
wolffd@0 156 <p>
wolffd@0 157 You can visualize the resulting graph structure using
wolffd@0 158 the methods discussed <a href="#graphdraw">below</a>.
wolffd@0 159 For details on GUIs,
wolffd@0 160 click <a href="#GUI">here</a>.
wolffd@0 161
wolffd@0 162 <h2>Creating the Bayes net shell</h2>
wolffd@0 163
wolffd@0 164 In addition to specifying the graph structure,
wolffd@0 165 we must specify the size and type of each node.
wolffd@0 166 If a node is discrete, its size is the
wolffd@0 167 number of possible values
wolffd@0 168 each node can take on; if a node is continuous,
wolffd@0 169 it can be a vector, and its size is the length of this vector.
wolffd@0 170 In this case, we will assume all nodes are discrete and binary.
wolffd@0 171 <PRE>
wolffd@0 172 discrete_nodes = 1:N;
wolffd@0 173 node_sizes = 2*ones(1,N);
wolffd@0 174 </pre>
wolffd@0 175 If the nodes were not binary, you could type e.g.,
wolffd@0 176 <pre>
wolffd@0 177 node_sizes = [4 2 3 5];
wolffd@0 178 </pre>
wolffd@0 179 meaning that Cloudy has 4 possible values,
wolffd@0 180 Sprinkler has 2 possible values, etc.
wolffd@0 181 Note that these are cardinal values, not ordinal, i.e.,
wolffd@0 182 they are not ordered in any way, like 'low', 'medium', 'high'.
wolffd@0 183 <p>
wolffd@0 184 We are now ready to make the Bayes net:
wolffd@0 185 <pre>
wolffd@0 186 bnet = mk_bnet(dag, node_sizes, 'discrete', discrete_nodes);
wolffd@0 187 </PRE>
wolffd@0 188 By default, all nodes are assumed to be discrete, so we can also just
wolffd@0 189 write
wolffd@0 190 <pre>
wolffd@0 191 bnet = mk_bnet(dag, node_sizes);
wolffd@0 192 </PRE>
wolffd@0 193 You may also specify which nodes will be observed.
wolffd@0 194 If you don't know, or if this not fixed in advance,
wolffd@0 195 just use the empty list (the default).
wolffd@0 196 <pre>
wolffd@0 197 onodes = [];
wolffd@0 198 bnet = mk_bnet(dag, node_sizes, 'discrete', discrete_nodes, 'observed', onodes);
wolffd@0 199 </PRE>
wolffd@0 200 Note that optional arguments are specified using a name/value syntax.
wolffd@0 201 This is common for many BNT functions.
wolffd@0 202 In general, to find out more about a function (e.g., which optional
wolffd@0 203 arguments it takes), please see its
wolffd@0 204 documentation string by typing
wolffd@0 205 <pre>
wolffd@0 206 help mk_bnet
wolffd@0 207 </pre>
wolffd@0 208 See also other <a href="matlab_tips.html">useful Matlab tips</a>.
wolffd@0 209 <p>
wolffd@0 210 It is possible to associate names with nodes, as follows:
wolffd@0 211 <pre>
wolffd@0 212 bnet = mk_bnet(dag, node_sizes, 'names', {'cloudy','S','R','W'}, 'discrete', 1:4);
wolffd@0 213 </pre>
wolffd@0 214 You can then refer to a node by its name:
wolffd@0 215 <pre>
wolffd@0 216 C = bnet.names('cloudy'); % bnet.names is an associative array
wolffd@0 217 bnet.CPD{C} = tabular_CPD(bnet, C, [0.5 0.5]);
wolffd@0 218 </pre>
wolffd@0 219 This feature uses my own associative array class.
wolffd@0 220
wolffd@0 221
wolffd@0 222 <h2><a name="cpt">Parameters</h2>
wolffd@0 223
wolffd@0 224 A model consists of the graph structure and the parameters.
wolffd@0 225 The parameters are represented by CPD objects (CPD = Conditional
wolffd@0 226 Probability Distribution), which define the probability distribution
wolffd@0 227 of a node given its parents.
wolffd@0 228 (We will use the terms "node" and "random variable" interchangeably.)
wolffd@0 229 The simplest kind of CPD is a table (multi-dimensional array), which
wolffd@0 230 is suitable when all the nodes are discrete-valued. Note that the discrete
wolffd@0 231 values are not assumed to be ordered in any way; that is, they
wolffd@0 232 represent categorical quantities, like male and female, rather than
wolffd@0 233 ordinal quantities, like low, medium and high.
wolffd@0 234 (We will discuss CPDs in more detail <a href="#cpd">below</a>.)
wolffd@0 235 <p>
wolffd@0 236 Tabular CPDs, also called CPTs (conditional probability tables),
wolffd@0 237 are stored as multidimensional arrays, where the dimensions
wolffd@0 238 are arranged in the same order as the nodes, e.g., the CPT for node 4
wolffd@0 239 (WetGrass) is indexed by Sprinkler (2), Rain (3) and then WetGrass (4) itself.
wolffd@0 240 Hence the child is always the last dimension.
wolffd@0 241 If a node has no parents, its CPT is a column vector representing its
wolffd@0 242 prior.
wolffd@0 243 Note that in Matlab (unlike C), arrays are indexed
wolffd@0 244 from 1, and are layed out in memory such that the first index toggles
wolffd@0 245 fastest, e.g., the CPT for node 4 (WetGrass) is as follows
wolffd@0 246 <P>
wolffd@0 247 <P><IMG ALIGN=BOTTOM SRC="Figures/CPTgrass.gif"><P>
wolffd@0 248 <P>
wolffd@0 249 where we have used the convention that false==1, true==2.
wolffd@0 250 We can create this CPT in Matlab as follows
wolffd@0 251 <PRE>
wolffd@0 252 CPT = zeros(2,2,2);
wolffd@0 253 CPT(1,1,1) = 1.0;
wolffd@0 254 CPT(2,1,1) = 0.1;
wolffd@0 255 ...
wolffd@0 256 </PRE>
wolffd@0 257 Here is an easier way:
wolffd@0 258 <PRE>
wolffd@0 259 CPT = reshape([1 0.1 0.1 0.01 0 0.9 0.9 0.99], [2 2 2]);
wolffd@0 260 </PRE>
wolffd@0 261 In fact, we don't need to reshape the array, since the CPD constructor
wolffd@0 262 will do that for us. So we can just write
wolffd@0 263 <pre>
wolffd@0 264 bnet.CPD{W} = tabular_CPD(bnet, W, 'CPT', [1 0.1 0.1 0.01 0 0.9 0.9 0.99]);
wolffd@0 265 </pre>
wolffd@0 266 The other nodes are created similarly (using the old syntax for
wolffd@0 267 optional parameters)
wolffd@0 268 <PRE>
wolffd@0 269 bnet.CPD{C} = tabular_CPD(bnet, C, [0.5 0.5]);
wolffd@0 270 bnet.CPD{R} = tabular_CPD(bnet, R, [0.8 0.2 0.2 0.8]);
wolffd@0 271 bnet.CPD{S} = tabular_CPD(bnet, S, [0.5 0.9 0.5 0.1]);
wolffd@0 272 bnet.CPD{W} = tabular_CPD(bnet, W, [1 0.1 0.1 0.01 0 0.9 0.9 0.99]);
wolffd@0 273 </PRE>
wolffd@0 274
wolffd@0 275
wolffd@0 276 <h2><a name="rnd_cpt">Random Parameters</h2>
wolffd@0 277
wolffd@0 278 If we do not specify the CPT, random parameters will be
wolffd@0 279 created, i.e., each "row" of the CPT will be drawn from the uniform distribution.
wolffd@0 280 To ensure repeatable results, use
wolffd@0 281 <pre>
wolffd@0 282 rand('state', seed);
wolffd@0 283 randn('state', seed);
wolffd@0 284 </pre>
wolffd@0 285 To control the degree of randomness (entropy),
wolffd@0 286 you can sample each row of the CPT from a Dirichlet(p,p,...) distribution.
wolffd@0 287 If p << 1, this encourages "deterministic" CPTs (one entry near 1, the rest near 0).
wolffd@0 288 If p = 1, each entry is drawn from U[0,1].
wolffd@0 289 If p >> 1, the entries will all be near 1/k, where k is the arity of
wolffd@0 290 this node, i.e., each row will be nearly uniform.
wolffd@0 291 You can do this as follows, assuming this node
wolffd@0 292 is number i, and ns is the node_sizes.
wolffd@0 293 <pre>
wolffd@0 294 k = ns(i);
wolffd@0 295 ps = parents(dag, i);
wolffd@0 296 psz = prod(ns(ps));
wolffd@0 297 CPT = sample_dirichlet(p*ones(1,k), psz);
wolffd@0 298 bnet.CPD{i} = tabular_CPD(bnet, i, 'CPT', CPT);
wolffd@0 299 </pre>
wolffd@0 300
wolffd@0 301
wolffd@0 302 <h2><a name="file">Loading a network from a file</h2>
wolffd@0 303
wolffd@0 304 If you already have a Bayes net represented in the XML-based
wolffd@0 305 <a href="http://www.cs.cmu.edu/afs/cs/user/fgcozman/www/Research/InterchangeFormat/">
wolffd@0 306 Bayes Net Interchange Format (BNIF)</a> (e.g., downloaded from the
wolffd@0 307 <a
wolffd@0 308 href="http://www.cs.huji.ac.il/labs/compbio/Repository">
wolffd@0 309 Bayes Net repository</a>),
wolffd@0 310 you can convert it to BNT format using
wolffd@0 311 the
wolffd@0 312 <a href="http://www.digitas.harvard.edu/~ken/bif2bnt/">BIF-BNT Java
wolffd@0 313 program</a> written by Ken Shan.
wolffd@0 314 (This is not necessarily up-to-date.)
wolffd@0 315 <p>
wolffd@0 316 <b>It is currently not possible to save/load a BNT matlab object to
wolffd@0 317 file</b>, but this is easily fixed if you modify all the constructors
wolffd@0 318 for all the classes (see matlab documentation).
wolffd@0 319
wolffd@0 320 <h2><a name="GUI">Creating a model using a GUI</h2>
wolffd@0 321
wolffd@0 322 <ul>
wolffd@0 323 <li>Senthil Nachimuthu
wolffd@0 324 has started (Oct 07) an open source
wolffd@0 325 GUI for BNT called
wolffd@0 326 <a href="http://projeny.sourceforge.net">projeny</a>
wolffd@0 327 using Java. This is a successor to BNJ.
wolffd@0 328
wolffd@0 329 <li>
wolffd@0 330 Philippe LeRay has written (Sep 05)
wolffd@0 331 a
wolffd@0 332 <a href="http://banquiseasi.insa-rouen.fr/projects/bnt-editor/">
wolffd@0 333 BNT GUI</a> in matlab.
wolffd@0 334
wolffd@0 335 <li>
wolffd@0 336 <a
wolffd@0 337 href="http://www.dataonstage.com/BNT/PACKAGES/LinkStrength/index.html">
wolffd@0 338 LinkStrength</a>,
wolffd@0 339 a package by Imme Ebert-Uphoff for visualizing the strength of
wolffd@0 340 dependencies between nodes.
wolffd@0 341
wolffd@0 342 </ul>
wolffd@0 343
wolffd@0 344 <h2>Graph visualization</h2>
wolffd@0 345
wolffd@0 346 Click <a href="graphviz.html">here</a> for more information
wolffd@0 347 on graph visualization.
wolffd@0 348
wolffd@0 349
wolffd@0 350 <h1><a name="inference">Inference</h1>
wolffd@0 351
wolffd@0 352 Having created the BN, we can now use it for inference.
wolffd@0 353 There are many different algorithms for doing inference in Bayes nets,
wolffd@0 354 that make different tradeoffs between speed,
wolffd@0 355 complexity, generality, and accuracy.
wolffd@0 356 BNT therefore offers a variety of different inference
wolffd@0 357 "engines". We will discuss these
wolffd@0 358 in more detail <a href="#engines">below</a>.
wolffd@0 359 For now, we will use the junction tree
wolffd@0 360 engine, which is the mother of all exact inference algorithms.
wolffd@0 361 This can be created as follows.
wolffd@0 362 <pre>
wolffd@0 363 engine = jtree_inf_engine(bnet);
wolffd@0 364 </pre>
wolffd@0 365 The other engines have similar constructors, but might take
wolffd@0 366 additional, algorithm-specific parameters.
wolffd@0 367 All engines are used in the same way, once they have been created.
wolffd@0 368 We illustrate this in the following sections.
wolffd@0 369
wolffd@0 370
wolffd@0 371 <h2><a name="marginal">Computing marginal distributions</h2>
wolffd@0 372
wolffd@0 373 Suppose we want to compute the probability that the sprinker was on
wolffd@0 374 given that the grass is wet.
wolffd@0 375 The evidence consists of the fact that W=2. All the other nodes
wolffd@0 376 are hidden (unobserved). We can specify this as follows.
wolffd@0 377 <pre>
wolffd@0 378 evidence = cell(1,N);
wolffd@0 379 evidence{W} = 2;
wolffd@0 380 </pre>
wolffd@0 381 We use a 1D cell array instead of a vector to
wolffd@0 382 cope with the fact that nodes can be vectors of different lengths.
wolffd@0 383 In addition, the value [] can be used
wolffd@0 384 to denote 'no evidence', instead of having to specify the observation
wolffd@0 385 pattern as a separate argument.
wolffd@0 386 (Click <a href="cellarray.html">here</a> for a quick tutorial on cell
wolffd@0 387 arrays in matlab.)
wolffd@0 388 <p>
wolffd@0 389 We are now ready to add the evidence to the engine.
wolffd@0 390 <pre>
wolffd@0 391 [engine, loglik] = enter_evidence(engine, evidence);
wolffd@0 392 </pre>
wolffd@0 393 The behavior of this function is algorithm-specific, and is discussed
wolffd@0 394 in more detail <a href="#engines">below</a>.
wolffd@0 395 In the case of the jtree engine,
wolffd@0 396 enter_evidence implements a two-pass message-passing scheme.
wolffd@0 397 The first return argument contains the modified engine, which
wolffd@0 398 incorporates the evidence. The second return argument contains the
wolffd@0 399 log-likelihood of the evidence. (Not all engines are capable of
wolffd@0 400 computing the log-likelihood.)
wolffd@0 401 <p>
wolffd@0 402 Finally, we can compute p=P(S=2|W=2) as follows.
wolffd@0 403 <PRE>
wolffd@0 404 marg = marginal_nodes(engine, S);
wolffd@0 405 marg.T
wolffd@0 406 ans =
wolffd@0 407 0.57024
wolffd@0 408 0.42976
wolffd@0 409 p = marg.T(2);
wolffd@0 410 </PRE>
wolffd@0 411 We see that p = 0.4298.
wolffd@0 412 <p>
wolffd@0 413 Now let us add the evidence that it was raining, and see what
wolffd@0 414 difference it makes.
wolffd@0 415 <PRE>
wolffd@0 416 evidence{R} = 2;
wolffd@0 417 [engine, loglik] = enter_evidence(engine, evidence);
wolffd@0 418 marg = marginal_nodes(engine, S);
wolffd@0 419 p = marg.T(2);
wolffd@0 420 </PRE>
wolffd@0 421 We find that p = P(S=2|W=2,R=2) = 0.1945,
wolffd@0 422 which is lower than
wolffd@0 423 before, because the rain can ``explain away'' the
wolffd@0 424 fact that the grass is wet.
wolffd@0 425 <p>
wolffd@0 426 You can plot a marginal distribution over a discrete variable
wolffd@0 427 as a barchart using the built 'bar' function:
wolffd@0 428 <pre>
wolffd@0 429 bar(marg.T)
wolffd@0 430 </pre>
wolffd@0 431 This is what it looks like
wolffd@0 432
wolffd@0 433 <p>
wolffd@0 434 <center>
wolffd@0 435 <IMG SRC="Figures/sprinkler_bar.gif">
wolffd@0 436 </center>
wolffd@0 437 <p>
wolffd@0 438
wolffd@0 439 <h2><a name="observed">Observed nodes</h2>
wolffd@0 440
wolffd@0 441 What happens if we ask for the marginal on an observed node, e.g. P(W|W=2)?
wolffd@0 442 An observed discrete node effectively only has 1 value (the observed
wolffd@0 443 one) --- all other values would result in 0 probability.
wolffd@0 444 For efficiency, BNT treats observed (discrete) nodes as if they were
wolffd@0 445 set to 1, as we see below:
wolffd@0 446 <pre>
wolffd@0 447 evidence = cell(1,N);
wolffd@0 448 evidence{W} = 2;
wolffd@0 449 engine = enter_evidence(engine, evidence);
wolffd@0 450 m = marginal_nodes(engine, W);
wolffd@0 451 m.T
wolffd@0 452 ans =
wolffd@0 453 1
wolffd@0 454 </pre>
wolffd@0 455 This can get a little confusing, since we assigned W=2.
wolffd@0 456 So we can ask BNT to add the evidence back in by passing in an optional argument:
wolffd@0 457 <pre>
wolffd@0 458 m = marginal_nodes(engine, W, 1);
wolffd@0 459 m.T
wolffd@0 460 ans =
wolffd@0 461 0
wolffd@0 462 1
wolffd@0 463 </pre>
wolffd@0 464 This shows that P(W=1|W=2) = 0 and P(W=2|W=2) = 1.
wolffd@0 465
wolffd@0 466
wolffd@0 467
wolffd@0 468 <h2><a name="joint">Computing joint distributions</h2>
wolffd@0 469
wolffd@0 470 We can compute the joint probability on a set of nodes as in the
wolffd@0 471 following example.
wolffd@0 472 <pre>
wolffd@0 473 evidence = cell(1,N);
wolffd@0 474 [engine, ll] = enter_evidence(engine, evidence);
wolffd@0 475 m = marginal_nodes(engine, [S R W]);
wolffd@0 476 </pre>
wolffd@0 477 m is a structure. The 'T' field is a multi-dimensional array (in
wolffd@0 478 this case, 3-dimensional) that contains the joint probability
wolffd@0 479 distribution on the specified nodes.
wolffd@0 480 <pre>
wolffd@0 481 >> m.T
wolffd@0 482 ans(:,:,1) =
wolffd@0 483 0.2900 0.0410
wolffd@0 484 0.0210 0.0009
wolffd@0 485 ans(:,:,2) =
wolffd@0 486 0 0.3690
wolffd@0 487 0.1890 0.0891
wolffd@0 488 </pre>
wolffd@0 489 We see that P(S=1,R=1,W=2) = 0, since it is impossible for the grass
wolffd@0 490 to be wet if both the rain and sprinkler are off.
wolffd@0 491 <p>
wolffd@0 492 Let us now add some evidence to R.
wolffd@0 493 <pre>
wolffd@0 494 evidence{R} = 2;
wolffd@0 495 [engine, ll] = enter_evidence(engine, evidence);
wolffd@0 496 m = marginal_nodes(engine, [S R W])
wolffd@0 497 m =
wolffd@0 498 domain: [2 3 4]
wolffd@0 499 T: [2x1x2 double]
wolffd@0 500 >> m.T
wolffd@0 501 m.T
wolffd@0 502 ans(:,:,1) =
wolffd@0 503 0.0820
wolffd@0 504 0.0018
wolffd@0 505 ans(:,:,2) =
wolffd@0 506 0.7380
wolffd@0 507 0.1782
wolffd@0 508 </pre>
wolffd@0 509 The joint T(i,j,k) = P(S=i,R=j,W=k|evidence)
wolffd@0 510 should have T(i,1,k) = 0 for all i,k, since R=1 is incompatible
wolffd@0 511 with the evidence that R=2.
wolffd@0 512 Instead of creating large tables with many 0s, BNT sets the effective
wolffd@0 513 size of observed (discrete) nodes to 1, as explained above.
wolffd@0 514 This is why m.T has size 2x1x2.
wolffd@0 515 To get a 2x2x2 table, type
wolffd@0 516 <pre>
wolffd@0 517 m = marginal_nodes(engine, [S R W], 1)
wolffd@0 518 m =
wolffd@0 519 domain: [2 3 4]
wolffd@0 520 T: [2x2x2 double]
wolffd@0 521 >> m.T
wolffd@0 522 m.T
wolffd@0 523 ans(:,:,1) =
wolffd@0 524 0 0.082
wolffd@0 525 0 0.0018
wolffd@0 526 ans(:,:,2) =
wolffd@0 527 0 0.738
wolffd@0 528 0 0.1782
wolffd@0 529 </pre>
wolffd@0 530
wolffd@0 531 <p>
wolffd@0 532 Note: It is not always possible to compute the joint on arbitrary
wolffd@0 533 sets of nodes: it depends on which inference engine you use, as discussed
wolffd@0 534 in more detail <a href="#engines">below</a>.
wolffd@0 535
wolffd@0 536
wolffd@0 537 <h2><a name="soft">Soft/virtual evidence</h2>
wolffd@0 538
wolffd@0 539 Sometimes a node is not observed, but we have some distribution over
wolffd@0 540 its possible values; this is often called "soft" or "virtual"
wolffd@0 541 evidence.
wolffd@0 542 One can use this as follows
wolffd@0 543 <pre>
wolffd@0 544 [engine, loglik] = enter_evidence(engine, evidence, 'soft', soft_evidence);
wolffd@0 545 </pre>
wolffd@0 546 where soft_evidence{i} is either [] (if node i has no soft evidence)
wolffd@0 547 or is a vector representing the probability distribution over i's
wolffd@0 548 possible values.
wolffd@0 549 For example, if we don't know i's exact value, but we know its
wolffd@0 550 likelihood ratio is 60/40, we can write evidence{i} = [] and
wolffd@0 551 soft_evidence{i} = [0.6 0.4].
wolffd@0 552 <p>
wolffd@0 553 Currently only jtree_inf_engine supports this option.
wolffd@0 554 It assumes that all hidden nodes, and all nodes for
wolffd@0 555 which we have soft evidence, are discrete.
wolffd@0 556 For a longer example, see BNT/examples/static/softev1.m.
wolffd@0 557
wolffd@0 558
wolffd@0 559 <h2><a name="mpe">Most probable explanation</h2>
wolffd@0 560
wolffd@0 561 To compute the most probable explanation (MPE) of the evidence (i.e.,
wolffd@0 562 the most probable assignment, or a mode of the joint), use
wolffd@0 563 <pre>
wolffd@0 564 [mpe, ll] = calc_mpe(engine, evidence);
wolffd@0 565 </pre>
wolffd@0 566 mpe{i} is the most likely value of node i.
wolffd@0 567 This calls enter_evidence with the 'maximize' flag set to 1, which
wolffd@0 568 causes the engine to do max-product instead of sum-product.
wolffd@0 569 The resulting max-marginals are then thresholded.
wolffd@0 570 If there is more than one maximum probability assignment, we must take
wolffd@0 571 care to break ties in a consistent manner (thresholding the
wolffd@0 572 max-marginals may give the wrong result). To force this behavior,
wolffd@0 573 type
wolffd@0 574 <pre>
wolffd@0 575 [mpe, ll] = calc_mpe(engine, evidence, 1);
wolffd@0 576 </pre>
wolffd@0 577 Note that computing the MPE is someties called abductive reasoning.
wolffd@0 578
wolffd@0 579 <p>
wolffd@0 580 You can also use <tt>calc_mpe_bucket</tt> written by Ron Zohar,
wolffd@0 581 that does a forwards max-product pass, and then a backwards traceback
wolffd@0 582 pass, which is how Viterbi is traditionally implemented.
wolffd@0 583
wolffd@0 584
wolffd@0 585
wolffd@0 586 <h1><a name="cpd">Conditional Probability Distributions</h1>
wolffd@0 587
wolffd@0 588 A Conditional Probability Distributions (CPD)
wolffd@0 589 defines P(X(i) | X(Pa(i))), where X(i) is the i'th node, and X(Pa(i))
wolffd@0 590 are the parents of node i. There are many ways to represent this
wolffd@0 591 distribution, which depend in part on whether X(i) and X(Pa(i)) are
wolffd@0 592 discrete, continuous, or a combination.
wolffd@0 593 We will discuss various representations below.
wolffd@0 594
wolffd@0 595
wolffd@0 596 <h2><a name="tabular">Tabular nodes</h2>
wolffd@0 597
wolffd@0 598 If the CPD is represented as a table (i.e., if it is a multinomial
wolffd@0 599 distribution), it has a number of parameters that is exponential in
wolffd@0 600 the number of parents. See the example <a href="#cpt">above</a>.
wolffd@0 601
wolffd@0 602
wolffd@0 603 <h2><a name="noisyor">Noisy-or nodes</h2>
wolffd@0 604
wolffd@0 605 A noisy-OR node is like a regular logical OR gate except that
wolffd@0 606 sometimes the effects of parents that are on get inhibited.
wolffd@0 607 Let the prob. that parent i gets inhibited be q(i).
wolffd@0 608 Then a node, C, with 2 parents, A and B, has the following CPD, where
wolffd@0 609 we use F and T to represent off and on (1 and 2 in BNT).
wolffd@0 610 <pre>
wolffd@0 611 A B P(C=off) P(C=on)
wolffd@0 612 ---------------------------
wolffd@0 613 F F 1.0 0.0
wolffd@0 614 T F q(A) 1-q(A)
wolffd@0 615 F T q(B) 1-q(B)
wolffd@0 616 T T q(A)q(B) 1-q(A)q(B)
wolffd@0 617 </pre>
wolffd@0 618 Thus we see that the causes get inhibited independently.
wolffd@0 619 It is common to associate a "leak" node with a noisy-or CPD, which is
wolffd@0 620 like a parent that is always on. This can account for all other unmodelled
wolffd@0 621 causes which might turn the node on.
wolffd@0 622 <p>
wolffd@0 623 The noisy-or distribution is similar to the logistic distribution.
wolffd@0 624 To see this, let the nodes, S(i), have values in {0,1}, and let q(i,j)
wolffd@0 625 be the prob. that j inhibits i. Then
wolffd@0 626 <pre>
wolffd@0 627 Pr(S(i)=1 | parents(S(i))) = 1 - prod_{j} q(i,j)^S(j)
wolffd@0 628 </pre>
wolffd@0 629 Now define w(i,j) = -ln q(i,j) and rho(x) = 1-exp(-x). Then
wolffd@0 630 <pre>
wolffd@0 631 Pr(S(i)=1 | parents(S(i))) = rho(sum_j w(i,j) S(j))
wolffd@0 632 </pre>
wolffd@0 633 For a sigmoid node, we have
wolffd@0 634 <pre>
wolffd@0 635 Pr(S(i)=1 | parents(S(i))) = sigma(-sum_j w(i,j) S(j))
wolffd@0 636 </pre>
wolffd@0 637 where sigma(x) = 1/(1+exp(-x)). Hence they differ in the choice of
wolffd@0 638 the activation function (although both are monotonically increasing).
wolffd@0 639 In addition, in the case of a noisy-or, the weights are constrained to be
wolffd@0 640 positive, since they derive from probabilities q(i,j).
wolffd@0 641 In both cases, the number of parameters is <em>linear</em> in the
wolffd@0 642 number of parents, unlike the case of a multinomial distribution,
wolffd@0 643 where the number of parameters is exponential in the number of parents.
wolffd@0 644 We will see an example of noisy-OR nodes <a href="#qmr">below</a>.
wolffd@0 645
wolffd@0 646
wolffd@0 647 <h2><a name="deterministic">Other (noisy) deterministic nodes</h2>
wolffd@0 648
wolffd@0 649 Deterministic CPDs for discrete random variables can be created using
wolffd@0 650 the deterministic_CPD class. It is also possible to 'flip' the output
wolffd@0 651 of the function with some probability, to simulate noise.
wolffd@0 652 The boolean_CPD class is just a special case of a
wolffd@0 653 deterministic CPD, where the parents and child are all binary.
wolffd@0 654 <p>
wolffd@0 655 Both of these classes are just "syntactic sugar" for the tabular_CPD
wolffd@0 656 class.
wolffd@0 657
wolffd@0 658
wolffd@0 659
wolffd@0 660 <h2><a name="softmax">Softmax nodes</h2>
wolffd@0 661
wolffd@0 662 If we have a discrete node with a continuous parent,
wolffd@0 663 we can define its CPD using a softmax function
wolffd@0 664 (also known as the multinomial logit function).
wolffd@0 665 This acts like a soft thresholding operator, and is defined as follows:
wolffd@0 666 <pre>
wolffd@0 667 exp(w(:,i)'*x + b(i))
wolffd@0 668 Pr(Q=i | X=x) = -----------------------------
wolffd@0 669 sum_j exp(w(:,j)'*x + b(j))
wolffd@0 670
wolffd@0 671 </pre>
wolffd@0 672 The parameters of a softmax node, w(:,i) and b(i), i=1..|Q|, have the
wolffd@0 673 following interpretation: w(:,i)-w(:,j) is the normal vector to the
wolffd@0 674 decision boundary between classes i and j,
wolffd@0 675 and b(i)-b(j) is its offset (bias). For example, suppose
wolffd@0 676 X is a 2-vector, and Q is binary. Then
wolffd@0 677 <pre>
wolffd@0 678 w = [1 -1;
wolffd@0 679 0 0];
wolffd@0 680
wolffd@0 681 b = [0 0];
wolffd@0 682 </pre>
wolffd@0 683 means class 1 are points in the 2D plane with positive x coordinate,
wolffd@0 684 and class 2 are points in the 2D plane with negative x coordinate.
wolffd@0 685 If w has large magnitude, the decision boundary is sharp, otherwise it
wolffd@0 686 is soft.
wolffd@0 687 In the special case that Q is binary (0/1), the softmax function reduces to the logistic
wolffd@0 688 (sigmoid) function.
wolffd@0 689 <p>
wolffd@0 690 Fitting a softmax function can be done using the iteratively reweighted
wolffd@0 691 least squares (IRLS) algorithm.
wolffd@0 692 We use the implementation from
wolffd@0 693 <a href="http://www.ncrg.aston.ac.uk/netlab/">Netlab</a>.
wolffd@0 694 Note that since
wolffd@0 695 the softmax distribution is not in the exponential family, it does not
wolffd@0 696 have finite sufficient statistics, and hence we must store all the
wolffd@0 697 training data in uncompressed form.
wolffd@0 698 If this takes too much space, one should use online (stochastic) gradient
wolffd@0 699 descent (not implemented in BNT).
wolffd@0 700 <p>
wolffd@0 701 If a softmax node also has discrete parents,
wolffd@0 702 we use a different set of w/b parameters for each combination of
wolffd@0 703 parent values, as in the <a href="#gaussian">conditional linear
wolffd@0 704 Gaussian CPD</a>.
wolffd@0 705 This feature was implemented by Pierpaolo Brutti.
wolffd@0 706 He is currently extending it so that discrete parents can be treated
wolffd@0 707 as if they were continuous, by adding indicator variables to the X
wolffd@0 708 vector.
wolffd@0 709 <p>
wolffd@0 710 We will see an example of softmax nodes <a href="#mixexp">below</a>.
wolffd@0 711
wolffd@0 712
wolffd@0 713 <h2><a name="mlp">Neural network nodes</h2>
wolffd@0 714
wolffd@0 715 Pierpaolo Brutti has implemented the mlp_CPD class, which uses a multi layer perceptron
wolffd@0 716 to implement a mapping from continuous parents to discrete children,
wolffd@0 717 similar to the softmax function.
wolffd@0 718 (If there are also discrete parents, it creates a mixture of MLPs.)
wolffd@0 719 It uses code from <a
wolffd@0 720 href="http://www.ncrg.aston.ac.uk/netlab/">Netlab</a>.
wolffd@0 721 This is work in progress.
wolffd@0 722
wolffd@0 723 <h2><a name="root">Root nodes</h2>
wolffd@0 724
wolffd@0 725 A root node has no parents and no parameters; it can be used to model
wolffd@0 726 an observed, exogeneous input variable, i.e., one which is "outside"
wolffd@0 727 the model.
wolffd@0 728 This is useful for conditional density models.
wolffd@0 729 We will see an example of root nodes <a href="#mixexp">below</a>.
wolffd@0 730
wolffd@0 731
wolffd@0 732 <h2><a name="gaussian">Gaussian nodes</h2>
wolffd@0 733
wolffd@0 734 We now consider a distribution suitable for the continuous-valued nodes.
wolffd@0 735 Suppose the node is called Y, its continuous parents (if any) are
wolffd@0 736 called X, and its discrete parents (if any) are called Q.
wolffd@0 737 The distribution on Y is defined as follows:
wolffd@0 738 <pre>
wolffd@0 739 - no parents: Y ~ N(mu, Sigma)
wolffd@0 740 - cts parents : Y|X=x ~ N(mu + W x, Sigma)
wolffd@0 741 - discrete parents: Y|Q=i ~ N(mu(:,i), Sigma(:,:,i))
wolffd@0 742 - cts and discrete parents: Y|X=x,Q=i ~ N(mu(:,i) + W(:,:,i) * x, Sigma(:,:,i))
wolffd@0 743 </pre>
wolffd@0 744 where N(mu, Sigma) denotes a Normal distribution with mean mu and
wolffd@0 745 covariance Sigma. Let |X|, |Y| and |Q| denote the sizes of X, Y and Q
wolffd@0 746 respectively.
wolffd@0 747 If there are no discrete parents, |Q|=1; if there is
wolffd@0 748 more than one, then |Q| = a vector of the sizes of each discrete parent.
wolffd@0 749 If there are no continuous parents, |X|=0; if there is more than one,
wolffd@0 750 then |X| = the sum of their sizes.
wolffd@0 751 Then mu is a |Y|*|Q| vector, Sigma is a |Y|*|Y|*|Q| positive
wolffd@0 752 semi-definite matrix, and W is a |Y|*|X|*|Q| regression (weight)
wolffd@0 753 matrix.
wolffd@0 754 <p>
wolffd@0 755 We can create a Gaussian node with random parameters as follows.
wolffd@0 756 <pre>
wolffd@0 757 bnet.CPD{i} = gaussian_CPD(bnet, i);
wolffd@0 758 </pre>
wolffd@0 759 We can specify the value of one or more of the parameters as in the
wolffd@0 760 following example, in which |Y|=2, and |Q|=1.
wolffd@0 761 <pre>
wolffd@0 762 bnet.CPD{i} = gaussian_CPD(bnet, i, 'mean', [0; 0], 'weights', randn(Y,X), 'cov', eye(Y));
wolffd@0 763 </pre>
wolffd@0 764 <p>
wolffd@0 765 We will see an example of conditional linear Gaussian nodes <a
wolffd@0 766 href="#cg_model">below</a>.
wolffd@0 767 <p>
wolffd@0 768 <b>When learning Gaussians from data</b>, it is helpful to ensure the
wolffd@0 769 data has a small magnitde
wolffd@0 770 (see e.g., KPMstats/standardize) to prevent numerical problems.
wolffd@0 771 Unless you have a lot of data, it is also a very good idea to use
wolffd@0 772 diagonal instead of full covariance matrices.
wolffd@0 773 (BNT does not currently support spherical covariances, although it
wolffd@0 774 would be easy to add, since KPMstats/clg_Mstep supports this option;
wolffd@0 775 you would just need to modify gaussian_CPD/update_ess to accumulate
wolffd@0 776 weighted inner products.)
wolffd@0 777
wolffd@0 778
wolffd@0 779
wolffd@0 780 <h2><a name="nongauss">Other continuous distributions</h2>
wolffd@0 781
wolffd@0 782 Currently BNT does not support any CPDs for continuous nodes other
wolffd@0 783 than the Gaussian.
wolffd@0 784 However, you can use a mixture of Gaussians to
wolffd@0 785 approximate other continuous distributions. We will see some an example
wolffd@0 786 of this with the IFA model <a href="#pca">below</a>.
wolffd@0 787
wolffd@0 788
wolffd@0 789 <h2><a name="glm">Generalized linear model nodes</h2>
wolffd@0 790
wolffd@0 791 In the future, we may incorporate some of the functionality of
wolffd@0 792 <a href =
wolffd@0 793 "http://www.sci.usq.edu.au/staff/dunn/glmlab/glmlab.html">glmlab</a>
wolffd@0 794 into BNT.
wolffd@0 795
wolffd@0 796
wolffd@0 797 <h2><a name="dtree">Classification/regression tree nodes</h2>
wolffd@0 798
wolffd@0 799 We plan to add classification and regression trees to define CPDs for
wolffd@0 800 discrete and continuous nodes, respectively.
wolffd@0 801 Trees have many advantages: they are easy to interpret, they can do
wolffd@0 802 feature selection, they can
wolffd@0 803 handle discrete and continuous inputs, they do not make strong
wolffd@0 804 assumptions about the form of the distribution, the number of
wolffd@0 805 parameters can grow in a data-dependent way (i.e., they are
wolffd@0 806 semi-parametric), they can handle missing data, etc.
wolffd@0 807 However, they are not yet implemented.
wolffd@0 808 <!--
wolffd@0 809 Yimin Zhang is currently (Feb '02) implementing this.
wolffd@0 810 -->
wolffd@0 811
wolffd@0 812
wolffd@0 813 <h2><a name="cpd_summary">Summary of CPD types</h2>
wolffd@0 814
wolffd@0 815 We list all the different types of CPDs supported by BNT.
wolffd@0 816 For each CPD, we specify if the child and parents can be discrete (D) or
wolffd@0 817 continuous (C) (Binary (B) nodes are a special case).
wolffd@0 818 We also specify which methods each class supports.
wolffd@0 819 If a method is inherited, the name of the parent class is mentioned.
wolffd@0 820 If a parent class calls a child method, this is mentioned.
wolffd@0 821 <p>
wolffd@0 822 The <tt>CPD_to_CPT</tt> method converts a CPD to a table; this
wolffd@0 823 requires that the child and all parents are discrete.
wolffd@0 824 The CPT might be exponentially big...
wolffd@0 825 <tt>convert_to_table</tt> evaluates a CPD with evidence, and
wolffd@0 826 represents the the resulting potential as an array.
wolffd@0 827 This requires that the child is discrete, and any continuous parents
wolffd@0 828 are observed.
wolffd@0 829 <tt>convert_to_pot</tt> evaluates a CPD with evidence, and
wolffd@0 830 represents the resulting potential as a dpot, gpot, cgpot or upot, as
wolffd@0 831 requested. (d=discrete, g=Gaussian, cg = conditional Gaussian, u =
wolffd@0 832 utility).
wolffd@0 833
wolffd@0 834 <p>
wolffd@0 835 When we sample a node, all the parents are observed.
wolffd@0 836 When we compute the (log) probability of a node, all the parents and
wolffd@0 837 the child are observed.
wolffd@0 838 <p>
wolffd@0 839 We also specify if the parameters are learnable.
wolffd@0 840 For learning with EM, we require
wolffd@0 841 the methods <tt>reset_ess</tt>, <tt>update_ess</tt> and
wolffd@0 842 <tt>maximize_params</tt>.
wolffd@0 843 For learning from fully observed data, we require
wolffd@0 844 the method <tt>learn_params</tt>.
wolffd@0 845 By default, all classes inherit this from generic_CPD, which simply
wolffd@0 846 calls <tt>update_ess</tt> N times, once for each data case, followed
wolffd@0 847 by <tt>maximize_params</tt>, i.e., it is like EM, without the E step.
wolffd@0 848 Some classes implement a batch formula, which is quicker.
wolffd@0 849 <p>
wolffd@0 850 Bayesian learning means computing a posterior over the parameters
wolffd@0 851 given fully observed data.
wolffd@0 852 <p>
wolffd@0 853 Pearl means we implement the methods <tt>compute_pi</tt> and
wolffd@0 854 <tt>compute_lambda_msg</tt>, used by
wolffd@0 855 <tt>pearl_inf_engine</tt>, which runs on directed graphs.
wolffd@0 856 <tt>belprop_inf_engine</tt> only needs <tt>convert_to_pot</tt>.H
wolffd@0 857 The pearl methods can exploit special properties of the CPDs for
wolffd@0 858 computing the messages efficiently, whereas belprop does not.
wolffd@0 859 <p>
wolffd@0 860 The only method implemented by generic_CPD is <tt>adjustable_CPD</tt>,
wolffd@0 861 which is not shown, since it is not very interesting.
wolffd@0 862
wolffd@0 863
wolffd@0 864 <p>
wolffd@0 865
wolffd@0 866
wolffd@0 867 <table>
wolffd@0 868 <table border units = pixels><tr>
wolffd@0 869 <td align=center>Name
wolffd@0 870 <td align=center>Child
wolffd@0 871 <td align=center>Parents
wolffd@0 872 <td align=center>Comments
wolffd@0 873 <td align=center>CPD_to_CPT
wolffd@0 874 <td align=center>conv_to_table
wolffd@0 875 <td align=center>conv_to_pot
wolffd@0 876 <td align=center>sample
wolffd@0 877 <td align=center>prob
wolffd@0 878 <td align=center>learn
wolffd@0 879 <td align=center>Bayes
wolffd@0 880 <td align=center>Pearl
wolffd@0 881
wolffd@0 882
wolffd@0 883 <tr>
wolffd@0 884 <!-- Name--><td>
wolffd@0 885 <!-- Child--><td>
wolffd@0 886 <!-- Parents--><td>
wolffd@0 887 <!-- Comments--><td>
wolffd@0 888 <!-- CPD_to_CPT--><td>
wolffd@0 889 <!-- conv_to_table--><td>
wolffd@0 890 <!-- conv_to_pot--><td>
wolffd@0 891 <!-- sample--><td>
wolffd@0 892 <!-- prob--><td>
wolffd@0 893 <!-- learn--><td>
wolffd@0 894 <!-- Bayes--><td>
wolffd@0 895 <!-- Pearl--><td>
wolffd@0 896
wolffd@0 897 <tr>
wolffd@0 898 <!-- Name--><td>boolean
wolffd@0 899 <!-- Child--><td>B
wolffd@0 900 <!-- Parents--><td>B
wolffd@0 901 <!-- Comments--><td>Syntactic sugar for tabular
wolffd@0 902 <!-- CPD_to_CPT--><td>-
wolffd@0 903 <!-- conv_to_table--><td>-
wolffd@0 904 <!-- conv_to_pot--><td>-
wolffd@0 905 <!-- sample--><td>-
wolffd@0 906 <!-- prob--><td>-
wolffd@0 907 <!-- learn--><td>-
wolffd@0 908 <!-- Bayes--><td>-
wolffd@0 909 <!-- Pearl--><td>-
wolffd@0 910
wolffd@0 911 <tr>
wolffd@0 912 <!-- Name--><td>deterministic
wolffd@0 913 <!-- Child--><td>D
wolffd@0 914 <!-- Parents--><td>D
wolffd@0 915 <!-- Comments--><td>Syntactic sugar for tabular
wolffd@0 916 <!-- CPD_to_CPT--><td>-
wolffd@0 917 <!-- conv_to_table--><td>-
wolffd@0 918 <!-- conv_to_pot--><td>-
wolffd@0 919 <!-- sample--><td>-
wolffd@0 920 <!-- prob--><td>-
wolffd@0 921 <!-- learn--><td>-
wolffd@0 922 <!-- Bayes--><td>-
wolffd@0 923 <!-- Pearl--><td>-
wolffd@0 924
wolffd@0 925 <tr>
wolffd@0 926 <!-- Name--><td>Discrete
wolffd@0 927 <!-- Child--><td>D
wolffd@0 928 <!-- Parents--><td>C/D
wolffd@0 929 <!-- Comments--><td>Virtual class
wolffd@0 930 <!-- CPD_to_CPT--><td>N
wolffd@0 931 <!-- conv_to_table--><td>Calls CPD_to_CPT
wolffd@0 932 <!-- conv_to_pot--><td>Calls conv_to_table
wolffd@0 933 <!-- sample--><td>Calls conv_to_table
wolffd@0 934 <!-- prob--><td>Calls conv_to_table
wolffd@0 935 <!-- learn--><td>N
wolffd@0 936 <!-- Bayes--><td>N
wolffd@0 937 <!-- Pearl--><td>N
wolffd@0 938
wolffd@0 939 <tr>
wolffd@0 940 <!-- Name--><td>Gaussian
wolffd@0 941 <!-- Child--><td>C
wolffd@0 942 <!-- Parents--><td>C/D
wolffd@0 943 <!-- Comments--><td>-
wolffd@0 944 <!-- CPD_to_CPT--><td>N
wolffd@0 945 <!-- conv_to_table--><td>N
wolffd@0 946 <!-- conv_to_pot--><td>Y
wolffd@0 947 <!-- sample--><td>Y
wolffd@0 948 <!-- prob--><td>Y
wolffd@0 949 <!-- learn--><td>Y
wolffd@0 950 <!-- Bayes--><td>N
wolffd@0 951 <!-- Pearl--><td>N
wolffd@0 952
wolffd@0 953 <tr>
wolffd@0 954 <!-- Name--><td>gmux
wolffd@0 955 <!-- Child--><td>C
wolffd@0 956 <!-- Parents--><td>C/D
wolffd@0 957 <!-- Comments--><td>multiplexer
wolffd@0 958 <!-- CPD_to_CPT--><td>N
wolffd@0 959 <!-- conv_to_table--><td>N
wolffd@0 960 <!-- conv_to_pot--><td>Y
wolffd@0 961 <!-- sample--><td>N
wolffd@0 962 <!-- prob--><td>N
wolffd@0 963 <!-- learn--><td>N
wolffd@0 964 <!-- Bayes--><td>N
wolffd@0 965 <!-- Pearl--><td>Y
wolffd@0 966
wolffd@0 967
wolffd@0 968 <tr>
wolffd@0 969 <!-- Name--><td>MLP
wolffd@0 970 <!-- Child--><td>D
wolffd@0 971 <!-- Parents--><td>C/D
wolffd@0 972 <!-- Comments--><td>multi layer perceptron
wolffd@0 973 <!-- CPD_to_CPT--><td>N
wolffd@0 974 <!-- conv_to_table--><td>Y
wolffd@0 975 <!-- conv_to_pot--><td>Inherits from discrete
wolffd@0 976 <!-- sample--><td>Inherits from discrete
wolffd@0 977 <!-- prob--><td>Inherits from discrete
wolffd@0 978 <!-- learn--><td>Y
wolffd@0 979 <!-- Bayes--><td>N
wolffd@0 980 <!-- Pearl--><td>N
wolffd@0 981
wolffd@0 982
wolffd@0 983 <tr>
wolffd@0 984 <!-- Name--><td>noisy-or
wolffd@0 985 <!-- Child--><td>B
wolffd@0 986 <!-- Parents--><td>B
wolffd@0 987 <!-- Comments--><td>-
wolffd@0 988 <!-- CPD_to_CPT--><td>Y
wolffd@0 989 <!-- conv_to_table--><td>Inherits from discrete
wolffd@0 990 <!-- conv_to_pot--><td>Inherits from discrete
wolffd@0 991 <!-- sample--><td>Inherits from discrete
wolffd@0 992 <!-- prob--><td>Inherits from discrete
wolffd@0 993 <!-- learn--><td>N
wolffd@0 994 <!-- Bayes--><td>N
wolffd@0 995 <!-- Pearl--><td>Y
wolffd@0 996
wolffd@0 997
wolffd@0 998 <tr>
wolffd@0 999 <!-- Name--><td>root
wolffd@0 1000 <!-- Child--><td>C/D
wolffd@0 1001 <!-- Parents--><td>none
wolffd@0 1002 <!-- Comments--><td>no params
wolffd@0 1003 <!-- CPD_to_CPT--><td>N
wolffd@0 1004 <!-- conv_to_table--><td>N
wolffd@0 1005 <!-- conv_to_pot--><td>Y
wolffd@0 1006 <!-- sample--><td>Y
wolffd@0 1007 <!-- prob--><td>Y
wolffd@0 1008 <!-- learn--><td>N
wolffd@0 1009 <!-- Bayes--><td>N
wolffd@0 1010 <!-- Pearl--><td>N
wolffd@0 1011
wolffd@0 1012
wolffd@0 1013 <tr>
wolffd@0 1014 <!-- Name--><td>softmax
wolffd@0 1015 <!-- Child--><td>D
wolffd@0 1016 <!-- Parents--><td>C/D
wolffd@0 1017 <!-- Comments--><td>-
wolffd@0 1018 <!-- CPD_to_CPT--><td>N
wolffd@0 1019 <!-- conv_to_table--><td>Y
wolffd@0 1020 <!-- conv_to_pot--><td>Inherits from discrete
wolffd@0 1021 <!-- sample--><td>Inherits from discrete
wolffd@0 1022 <!-- prob--><td>Inherits from discrete
wolffd@0 1023 <!-- learn--><td>Y
wolffd@0 1024 <!-- Bayes--><td>N
wolffd@0 1025 <!-- Pearl--><td>N
wolffd@0 1026
wolffd@0 1027
wolffd@0 1028 <tr>
wolffd@0 1029 <!-- Name--><td>generic
wolffd@0 1030 <!-- Child--><td>C/D
wolffd@0 1031 <!-- Parents--><td>C/D
wolffd@0 1032 <!-- Comments--><td>Virtual class
wolffd@0 1033 <!-- CPD_to_CPT--><td>N
wolffd@0 1034 <!-- conv_to_table--><td>N
wolffd@0 1035 <!-- conv_to_pot--><td>N
wolffd@0 1036 <!-- sample--><td>N
wolffd@0 1037 <!-- prob--><td>N
wolffd@0 1038 <!-- learn--><td>N
wolffd@0 1039 <!-- Bayes--><td>N
wolffd@0 1040 <!-- Pearl--><td>N
wolffd@0 1041
wolffd@0 1042
wolffd@0 1043 <tr>
wolffd@0 1044 <!-- Name--><td>Tabular
wolffd@0 1045 <!-- Child--><td>D
wolffd@0 1046 <!-- Parents--><td>D
wolffd@0 1047 <!-- Comments--><td>-
wolffd@0 1048 <!-- CPD_to_CPT--><td>Y
wolffd@0 1049 <!-- conv_to_table--><td>Inherits from discrete
wolffd@0 1050 <!-- conv_to_pot--><td>Inherits from discrete
wolffd@0 1051 <!-- sample--><td>Inherits from discrete
wolffd@0 1052 <!-- prob--><td>Inherits from discrete
wolffd@0 1053 <!-- learn--><td>Y
wolffd@0 1054 <!-- Bayes--><td>Y
wolffd@0 1055 <!-- Pearl--><td>Y
wolffd@0 1056
wolffd@0 1057 </table>
wolffd@0 1058
wolffd@0 1059
wolffd@0 1060
wolffd@0 1061 <h1><a name="examples">Example models</h1>
wolffd@0 1062
wolffd@0 1063
wolffd@0 1064 <h2>Gaussian mixture models</h2>
wolffd@0 1065
wolffd@0 1066 Richard W. DeVaul has made a detailed tutorial on how to fit mixtures
wolffd@0 1067 of Gaussians using BNT. Available
wolffd@0 1068 <a href="http://www.media.mit.edu/wearables/mithril/BNT/mixtureBNT.txt">here</a>.
wolffd@0 1069
wolffd@0 1070
wolffd@0 1071 <h2><a name="pca">PCA, ICA, and all that </h2>
wolffd@0 1072
wolffd@0 1073 In Figure (a) below, we show how Factor Analysis can be thought of as a
wolffd@0 1074 graphical model. Here, X has an N(0,I) prior, and
wolffd@0 1075 Y|X=x ~ N(mu + Wx, Psi),
wolffd@0 1076 where Psi is diagonal and W is called the "factor loading matrix".
wolffd@0 1077 Since the noise on both X and Y is diagonal, the components of these
wolffd@0 1078 vectors are uncorrelated, and hence can be represented as individual
wolffd@0 1079 scalar nodes, as we show in (b).
wolffd@0 1080 (This is useful if parts of the observations on the Y vector are occasionally missing.)
wolffd@0 1081 We usually take k=|X| << |Y|=D, so the model tries to explain
wolffd@0 1082 many observations using a low-dimensional subspace.
wolffd@0 1083
wolffd@0 1084
wolffd@0 1085 <center>
wolffd@0 1086 <table>
wolffd@0 1087 <tr>
wolffd@0 1088 <td><img src="Figures/fa.gif">
wolffd@0 1089 <td><img src="Figures/fa_scalar.gif">
wolffd@0 1090 <td><img src="Figures/mfa.gif">
wolffd@0 1091 <td><img src="Figures/ifa.gif">
wolffd@0 1092 <tr>
wolffd@0 1093 <td align=center> (a)
wolffd@0 1094 <td align=center> (b)
wolffd@0 1095 <td align=center> (c)
wolffd@0 1096 <td align=center> (d)
wolffd@0 1097 </table>
wolffd@0 1098 </center>
wolffd@0 1099
wolffd@0 1100 <p>
wolffd@0 1101 We can create this model in BNT as follows.
wolffd@0 1102 <pre>
wolffd@0 1103 ns = [k D];
wolffd@0 1104 dag = zeros(2,2);
wolffd@0 1105 dag(1,2) = 1;
wolffd@0 1106 bnet = mk_bnet(dag, ns, 'discrete', []);
wolffd@0 1107 bnet.CPD{1} = gaussian_CPD(bnet, 1, 'mean', zeros(k,1), 'cov', eye(k), ...
wolffd@0 1108 'cov_type', 'diag', 'clamp_mean', 1, 'clamp_cov', 1);
wolffd@0 1109 bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', zeros(D,1), 'cov', diag(Psi0), 'weights', W0, ...
wolffd@0 1110 'cov_type', 'diag', 'clamp_mean', 1);
wolffd@0 1111 </pre>
wolffd@0 1112
wolffd@0 1113 The root node is clamped to the N(0,I) distribution, so that we will
wolffd@0 1114 not update these parameters during learning.
wolffd@0 1115 The mean of the leaf node is clamped to 0,
wolffd@0 1116 since we assume the data has been centered (had its mean subtracted
wolffd@0 1117 off); this is just for simplicity.
wolffd@0 1118 Finally, the covariance of the leaf node is constrained to be
wolffd@0 1119 diagonal. W0 and Psi0 are the initial parameter guesses.
wolffd@0 1120
wolffd@0 1121 <p>
wolffd@0 1122 We can fit this model (i.e., estimate its parameters in a maximum
wolffd@0 1123 likelihood (ML) sense) using EM, as we
wolffd@0 1124 explain <a href="#em">below</a>.
wolffd@0 1125 Not surprisingly, the ML estimates for mu and Psi turn out to be
wolffd@0 1126 identical to the
wolffd@0 1127 sample mean and variance, which can be computed directly as
wolffd@0 1128 <pre>
wolffd@0 1129 mu_ML = mean(data);
wolffd@0 1130 Psi_ML = diag(cov(data));
wolffd@0 1131 </pre>
wolffd@0 1132 Note that W can only be identified up to a rotation matrix, because of
wolffd@0 1133 the spherical symmetry of the source.
wolffd@0 1134
wolffd@0 1135 <p>
wolffd@0 1136 If we restrict Psi to be spherical, i.e., Psi = sigma*I,
wolffd@0 1137 there is a closed-form solution for W as well,
wolffd@0 1138 i.e., we do not need to use EM.
wolffd@0 1139 In particular, W contains the first |X| eigenvectors of the sample covariance
wolffd@0 1140 matrix, with scalings determined by the eigenvalues and sigma.
wolffd@0 1141 Classical PCA can be obtained by taking the sigma->0 limit.
wolffd@0 1142 For details, see
wolffd@0 1143
wolffd@0 1144 <ul>
wolffd@0 1145 <li> <a href="ftp://hope.caltech.edu/pub/roweis/Empca/empca.ps">
wolffd@0 1146 "EM algorithms for PCA and SPCA"</a>, Sam Roweis, NIPS 97.
wolffd@0 1147 (<a href="ftp://hope.caltech.edu/pub/roweis/Code/empca.tar.gz">
wolffd@0 1148 Matlab software</a>)
wolffd@0 1149
wolffd@0 1150 <p>
wolffd@0 1151 <li>
wolffd@0 1152 <a
wolffd@0 1153 href=http://neural-server.aston.ac.uk/cgi-bin/tr_avail.pl?trnumber=NCRG/97/003>
wolffd@0 1154 "Mixtures of probabilistic principal component analyzers"</a>,
wolffd@0 1155 Tipping and Bishop, Neural Computation 11(2):443--482, 1999.
wolffd@0 1156 </ul>
wolffd@0 1157
wolffd@0 1158 <p>
wolffd@0 1159 By adding a hidden discrete variable, we can create mixtures of FA
wolffd@0 1160 models, as shown in (c).
wolffd@0 1161 Now we can explain the data using a set of subspaces.
wolffd@0 1162 We can create this model in BNT as follows.
wolffd@0 1163 <pre>
wolffd@0 1164 ns = [M k D];
wolffd@0 1165 dag = zeros(3);
wolffd@0 1166 dag(1,3) = 1;
wolffd@0 1167 dag(2,3) = 1;
wolffd@0 1168 bnet = mk_bnet(dag, ns, 'discrete', 1);
wolffd@0 1169 bnet.CPD{1} = tabular_CPD(bnet, 1, Pi0);
wolffd@0 1170 bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', zeros(k, 1), 'cov', eye(k), 'cov_type', 'diag', ...
wolffd@0 1171 'clamp_mean', 1, 'clamp_cov', 1);
wolffd@0 1172 bnet.CPD{3} = gaussian_CPD(bnet, 3, 'mean', Mu0', 'cov', repmat(diag(Psi0), [1 1 M]), ...
wolffd@0 1173 'weights', W0, 'cov_type', 'diag', 'tied_cov', 1);
wolffd@0 1174 </pre>
wolffd@0 1175 Notice how the covariance matrix for Y is the same for all values of
wolffd@0 1176 Q; that is, the noise level in each sub-space is assumed the same.
wolffd@0 1177 However, we allow the offset, mu, to vary.
wolffd@0 1178 For details, see
wolffd@0 1179 <ul>
wolffd@0 1180
wolffd@0 1181 <LI>
wolffd@0 1182 <a HREF="ftp://ftp.cs.toronto.edu/pub/zoubin/tr-96-1.ps.gz"> The EM
wolffd@0 1183 Algorithm for Mixtures of Factor Analyzers </A>,
wolffd@0 1184 Ghahramani, Z. and Hinton, G.E. (1996),
wolffd@0 1185 University of Toronto
wolffd@0 1186 Technical Report CRG-TR-96-1.
wolffd@0 1187 (<A HREF="ftp://ftp.cs.toronto.edu/pub/zoubin/mfa.tar.gz">Matlab software</A>)
wolffd@0 1188
wolffd@0 1189 <p>
wolffd@0 1190 <li>
wolffd@0 1191 <a
wolffd@0 1192 href=http://neural-server.aston.ac.uk/cgi-bin/tr_avail.pl?trnumber=NCRG/97/003>
wolffd@0 1193 "Mixtures of probabilistic principal component analyzers"</a>,
wolffd@0 1194 Tipping and Bishop, Neural Computation 11(2):443--482, 1999.
wolffd@0 1195 </ul>
wolffd@0 1196
wolffd@0 1197 <p>
wolffd@0 1198 I have included Zoubin's specialized MFA code (with his permission)
wolffd@0 1199 with the toolbox, so you can check that BNT gives the same results:
wolffd@0 1200 see 'BNT/examples/static/mfa1.m'.
wolffd@0 1201
wolffd@0 1202 <p>
wolffd@0 1203 Independent Factor Analysis (IFA) generalizes FA by allowing a
wolffd@0 1204 non-Gaussian prior on each component of X.
wolffd@0 1205 (Note that we can approximate a non-Gaussian prior using a mixture of
wolffd@0 1206 Gaussians.)
wolffd@0 1207 This means that the likelihood function is no longer rotationally
wolffd@0 1208 invariant, so we can uniquely identify W and the hidden
wolffd@0 1209 sources X.
wolffd@0 1210 IFA also allows a non-diagonal Psi (i.e. correlations between the components of Y).
wolffd@0 1211 We recover classical Independent Components Analysis (ICA)
wolffd@0 1212 in the Psi -> 0 limit, and by assuming that |X|=|Y|, so that the
wolffd@0 1213 weight matrix W is square and invertible.
wolffd@0 1214 For details, see
wolffd@0 1215 <ul>
wolffd@0 1216 <li>
wolffd@0 1217 <a href="http://www.gatsby.ucl.ac.uk/~hagai/ifa.ps">Independent Factor
wolffd@0 1218 Analysis</a>, H. Attias, Neural Computation 11: 803--851, 1998.
wolffd@0 1219 </ul>
wolffd@0 1220
wolffd@0 1221
wolffd@0 1222
wolffd@0 1223 <h2><a name="mixexp">Mixtures of experts</h2>
wolffd@0 1224
wolffd@0 1225 As an example of the use of the softmax function,
wolffd@0 1226 we introduce the Mixture of Experts model.
wolffd@0 1227 <!--
wolffd@0 1228 We also show
wolffd@0 1229 the Hierarchical Mixture of Experts model, where the hierarchy has two
wolffd@0 1230 levels.
wolffd@0 1231 (This is essentially a probabilistic decision tree of height two.)
wolffd@0 1232 -->
wolffd@0 1233 As before,
wolffd@0 1234 circles denote continuous-valued nodes,
wolffd@0 1235 squares denote discrete nodes, clear
wolffd@0 1236 means hidden, and shaded means observed.
wolffd@0 1237 <p>
wolffd@0 1238 <center>
wolffd@0 1239 <table>
wolffd@0 1240 <tr>
wolffd@0 1241 <td><img src="Figures/mixexp.gif">
wolffd@0 1242 <!--
wolffd@0 1243 <td><img src="Figures/hme.gif">
wolffd@0 1244 -->
wolffd@0 1245 </table>
wolffd@0 1246 </center>
wolffd@0 1247 <p>
wolffd@0 1248 X is the observed
wolffd@0 1249 input, Y is the output, and
wolffd@0 1250 the Q nodes are hidden "gating" nodes, which select the appropriate
wolffd@0 1251 set of parameters for Y. During training, Y is assumed observed,
wolffd@0 1252 but for testing, the goal is to predict Y given X.
wolffd@0 1253 Note that this is a <em>conditional</em> density model, so we don't
wolffd@0 1254 associate any parameters with X.
wolffd@0 1255 Hence X's CPD will be a root CPD, which is a way of modelling
wolffd@0 1256 exogenous nodes.
wolffd@0 1257 If the output is a continuous-valued quantity,
wolffd@0 1258 we assume the "experts" are linear-regression units,
wolffd@0 1259 and set Y's CPD to linear-Gaussian.
wolffd@0 1260 If the output is discrete, we set Y's CPD to a softmax function.
wolffd@0 1261 The Q CPDs will always be softmax functions.
wolffd@0 1262
wolffd@0 1263 <p>
wolffd@0 1264 As a concrete example, consider the mixture of experts model where X and Y are
wolffd@0 1265 scalars, and Q is binary.
wolffd@0 1266 This is just piecewise linear regression, where
wolffd@0 1267 we have two line segments, i.e.,
wolffd@0 1268 <P>
wolffd@0 1269 <IMG ALIGN=BOTTOM SRC="Eqns/lin_reg_eqn.gif">
wolffd@0 1270 <P>
wolffd@0 1271 We can create this model with random parameters as follows.
wolffd@0 1272 (This code is bundled in BNT/examples/static/mixexp2.m.)
wolffd@0 1273 <PRE>
wolffd@0 1274 X = 1;
wolffd@0 1275 Q = 2;
wolffd@0 1276 Y = 3;
wolffd@0 1277 dag = zeros(3,3);
wolffd@0 1278 dag(X,[Q Y]) = 1
wolffd@0 1279 dag(Q,Y) = 1;
wolffd@0 1280 ns = [1 2 1]; % make X and Y scalars, and have 2 experts
wolffd@0 1281 onodes = [1 3];
wolffd@0 1282 bnet = mk_bnet(dag, ns, 'discrete', 2, 'observed', onodes);
wolffd@0 1283
wolffd@0 1284 rand('state', 0);
wolffd@0 1285 randn('state', 0);
wolffd@0 1286 bnet.CPD{1} = root_CPD(bnet, 1);
wolffd@0 1287 bnet.CPD{2} = softmax_CPD(bnet, 2);
wolffd@0 1288 bnet.CPD{3} = gaussian_CPD(bnet, 3);
wolffd@0 1289 </PRE>
wolffd@0 1290 Now let us fit this model using <a href="#em">EM</a>.
wolffd@0 1291 First we <a href="#load_data">load the data</a> (1000 training cases) and plot them.
wolffd@0 1292 <P>
wolffd@0 1293 <PRE>
wolffd@0 1294 data = load('/examples/static/Misc/mixexp_data.txt', '-ascii');
wolffd@0 1295 plot(data(:,1), data(:,2), '.');
wolffd@0 1296 </PRE>
wolffd@0 1297 <p>
wolffd@0 1298 <center>
wolffd@0 1299 <IMG SRC="Figures/mixexp_data.gif">
wolffd@0 1300 </center>
wolffd@0 1301 <p>
wolffd@0 1302 This is what the model looks like before training.
wolffd@0 1303 (Thanks to Thomas Hofman for writing this plotting routine.)
wolffd@0 1304 <p>
wolffd@0 1305 <center>
wolffd@0 1306 <IMG SRC="Figures/mixexp_before.gif">
wolffd@0 1307 </center>
wolffd@0 1308 <p>
wolffd@0 1309 Now let's train the model, and plot the final performance.
wolffd@0 1310 (We will discuss how to train models in more detail <a href="#param_learning">below</a>.)
wolffd@0 1311 <P>
wolffd@0 1312 <PRE>
wolffd@0 1313 ncases = size(data, 1); % each row of data is a training case
wolffd@0 1314 cases = cell(3, ncases);
wolffd@0 1315 cases([1 3], :) = num2cell(data'); % each column of cases is a training case
wolffd@0 1316 engine = jtree_inf_engine(bnet);
wolffd@0 1317 max_iter = 20;
wolffd@0 1318 [bnet2, LLtrace] = learn_params_em(engine, cases, max_iter);
wolffd@0 1319 </PRE>
wolffd@0 1320 (We specify which nodes will be observed when we create the engine.
wolffd@0 1321 Hence BNT knows that the hidden nodes are all discrete.
wolffd@0 1322 For complex models, this can lead to a significant speedup.)
wolffd@0 1323 Below we show what the model looks like after 16 iterations of EM
wolffd@0 1324 (with 100 IRLS iterations per M step), when it converged
wolffd@0 1325 using the default convergence tolerance (that the
wolffd@0 1326 fractional change in the log-likelihood be less than 1e-3).
wolffd@0 1327 Before learning, the log-likelihood was
wolffd@0 1328 -322.927442; afterwards, it was -13.728778.
wolffd@0 1329 <p>
wolffd@0 1330 <center>
wolffd@0 1331 <IMG SRC="Figures/mixexp_after.gif">
wolffd@0 1332 </center>
wolffd@0 1333 (See BNT/examples/static/mixexp2.m for details of the code.)
wolffd@0 1334
wolffd@0 1335
wolffd@0 1336
wolffd@0 1337 <h2><a name="hme">Hierarchical mixtures of experts</h2>
wolffd@0 1338
wolffd@0 1339 A hierarchical mixture of experts (HME) extends the mixture of experts
wolffd@0 1340 model by having more than one hidden node. A two-level example is shown below, along
wolffd@0 1341 with its more traditional representation as a neural network.
wolffd@0 1342 This is like a (balanced) probabilistic decision tree of height 2.
wolffd@0 1343 <p>
wolffd@0 1344 <center>
wolffd@0 1345 <IMG SRC="Figures/HMEforMatlab.jpg">
wolffd@0 1346 </center>
wolffd@0 1347 <p>
wolffd@0 1348 <a href="mailto:pbrutti@stat.cmu.edu">Pierpaolo Brutti</a>
wolffd@0 1349 has written an extensive set of routines for HMEs,
wolffd@0 1350 which are bundled with BNT: see the examples/static/HME directory.
wolffd@0 1351 These routines allow you to choose the number of hidden (gating)
wolffd@0 1352 layers, and the form of the experts (softmax or MLP).
wolffd@0 1353 See the file hmemenu, which provides a demo.
wolffd@0 1354 For example, the figure below shows the decision boundaries learned
wolffd@0 1355 for a ternary classification problem, using a 2 level HME with softmax
wolffd@0 1356 gates and softmax experts; the training set is on the left, the
wolffd@0 1357 testing set on the right.
wolffd@0 1358 <p>
wolffd@0 1359 <center>
wolffd@0 1360 <!--<IMG SRC="Figures/hme_dec_boundary.gif">-->
wolffd@0 1361 <IMG SRC="Figures/hme_dec_boundary.png">
wolffd@0 1362 </center>
wolffd@0 1363 <p>
wolffd@0 1364
wolffd@0 1365
wolffd@0 1366 <p>
wolffd@0 1367 For more details, see the following:
wolffd@0 1368 <ul>
wolffd@0 1369
wolffd@0 1370 <li> <a href="http://www.cs.berkeley.edu/~jordan/papers/hierarchies.ps.Z">
wolffd@0 1371 Hierarchical mixtures of experts and the EM algorithm</a>
wolffd@0 1372 M. I. Jordan and R. A. Jacobs. Neural Computation, 6, 181-214, 1994.
wolffd@0 1373
wolffd@0 1374 <li> <a href =
wolffd@0 1375 "http://www.cs.berkeley.edu/~dmartin/software">David Martin's
wolffd@0 1376 matlab code for HME</a>
wolffd@0 1377
wolffd@0 1378 <li> <a
wolffd@0 1379 href="http://www.cs.berkeley.edu/~jordan/papers/uai.ps.Z">Why the
wolffd@0 1380 logistic function? A tutorial discussion on
wolffd@0 1381 probabilities and neural networks.</a> M. I. Jordan. MIT Computational
wolffd@0 1382 Cognitive Science Report 9503, August 1995.
wolffd@0 1383
wolffd@0 1384 <li> "Generalized Linear Models", McCullagh and Nelder, Chapman and
wolffd@0 1385 Halll, 1983.
wolffd@0 1386
wolffd@0 1387 <li>
wolffd@0 1388 "Improved learning algorithms for mixtures of experts in multiclass
wolffd@0 1389 classification".
wolffd@0 1390 K. Chen, L. Xu, H. Chi.
wolffd@0 1391 Neural Networks (1999) 12: 1229-1252.
wolffd@0 1392
wolffd@0 1393 <li> <a href="http://www.oigeeza.com/steve/">
wolffd@0 1394 Classification Using Hierarchical Mixtures of Experts</a>
wolffd@0 1395 S.R. Waterhouse and A.J. Robinson.
wolffd@0 1396 In Proc. IEEE Workshop on Neural Network for Signal Processing IV (1994), pp. 177-186
wolffd@0 1397
wolffd@0 1398 <li> <a href="http://www.idiap.ch/~perry/">
wolffd@0 1399 Localized mixtures of experts</a>,
wolffd@0 1400 P. Moerland, 1998.
wolffd@0 1401
wolffd@0 1402 <li> "Nonlinear gated experts for time series",
wolffd@0 1403 A.S. Weigend and M. Mangeas, 1995.
wolffd@0 1404
wolffd@0 1405 </ul>
wolffd@0 1406
wolffd@0 1407
wolffd@0 1408 <h2><a name="qmr">QMR</h2>
wolffd@0 1409
wolffd@0 1410 Bayes nets originally arose out of an attempt to add probabilities to
wolffd@0 1411 expert systems, and this is still the most common use for BNs.
wolffd@0 1412 A famous example is
wolffd@0 1413 QMR-DT, a decision-theoretic reformulation of the Quick Medical
wolffd@0 1414 Reference (QMR) model.
wolffd@0 1415 <p>
wolffd@0 1416 <center>
wolffd@0 1417 <IMG ALIGN=BOTTOM SRC="Figures/qmr.gif">
wolffd@0 1418 </center>
wolffd@0 1419 Here, the top layer represents hidden disease nodes, and the bottom
wolffd@0 1420 layer represents observed symptom nodes.
wolffd@0 1421 The goal is to infer the posterior probability of each disease given
wolffd@0 1422 all the symptoms (which can be present, absent or unknown).
wolffd@0 1423 Each node in the top layer has a Bernoulli prior (with a low prior
wolffd@0 1424 probability that the disease is present).
wolffd@0 1425 Since each node in the bottom layer has a high fan-in, we use a
wolffd@0 1426 noisy-OR parameterization; each disease has an independent chance of
wolffd@0 1427 causing each symptom.
wolffd@0 1428 The real QMR-DT model is copyright, but
wolffd@0 1429 we can create a random QMR-like model as follows.
wolffd@0 1430 <pre>
wolffd@0 1431 function bnet = mk_qmr_bnet(G, inhibit, leak, prior)
wolffd@0 1432 % MK_QMR_BNET Make a QMR model
wolffd@0 1433 % bnet = mk_qmr_bnet(G, inhibit, leak, prior)
wolffd@0 1434 %
wolffd@0 1435 % G(i,j) = 1 iff there is an arc from disease i to finding j
wolffd@0 1436 % inhibit(i,j) = inhibition probability on i->j arc
wolffd@0 1437 % leak(j) = inhibition prob. on leak->j arc
wolffd@0 1438 % prior(i) = prob. disease i is on
wolffd@0 1439
wolffd@0 1440 [Ndiseases Nfindings] = size(inhibit);
wolffd@0 1441 N = Ndiseases + Nfindings;
wolffd@0 1442 finding_node = Ndiseases+1:N;
wolffd@0 1443 ns = 2*ones(1,N);
wolffd@0 1444 dag = zeros(N,N);
wolffd@0 1445 dag(1:Ndiseases, finding_node) = G;
wolffd@0 1446 bnet = mk_bnet(dag, ns, 'observed', finding_node);
wolffd@0 1447
wolffd@0 1448 for d=1:Ndiseases
wolffd@0 1449 CPT = [1-prior(d) prior(d)];
wolffd@0 1450 bnet.CPD{d} = tabular_CPD(bnet, d, CPT');
wolffd@0 1451 end
wolffd@0 1452
wolffd@0 1453 for i=1:Nfindings
wolffd@0 1454 fnode = finding_node(i);
wolffd@0 1455 ps = parents(G, i);
wolffd@0 1456 bnet.CPD{fnode} = noisyor_CPD(bnet, fnode, leak(i), inhibit(ps, i));
wolffd@0 1457 end
wolffd@0 1458 </pre>
wolffd@0 1459 In the file BNT/examples/static/qmr1, we create a random bipartite
wolffd@0 1460 graph G, with 5 diseases and 10 findings, and random parameters.
wolffd@0 1461 (In general, to create a random dag, use 'mk_random_dag'.)
wolffd@0 1462 We can visualize the resulting graph structure using
wolffd@0 1463 the methods discussed <a href="#graphdraw">below</a>, with the
wolffd@0 1464 following results:
wolffd@0 1465 <p>
wolffd@0 1466 <img src="Figures/qmr.rnd.jpg">
wolffd@0 1467
wolffd@0 1468 <p>
wolffd@0 1469 Now let us put some random evidence on all the leaves except the very
wolffd@0 1470 first and very last, and compute the disease posteriors.
wolffd@0 1471 <pre>
wolffd@0 1472 pos = 2:floor(Nfindings/2);
wolffd@0 1473 neg = (pos(end)+1):(Nfindings-1);
wolffd@0 1474 onodes = myunion(pos, neg);
wolffd@0 1475 evidence = cell(1, N);
wolffd@0 1476 evidence(findings(pos)) = num2cell(repmat(2, 1, length(pos)));
wolffd@0 1477 evidence(findings(neg)) = num2cell(repmat(1, 1, length(neg)));
wolffd@0 1478
wolffd@0 1479 engine = jtree_inf_engine(bnet);
wolffd@0 1480 [engine, ll] = enter_evidence(engine, evidence);
wolffd@0 1481 post = zeros(1, Ndiseases);
wolffd@0 1482 for i=diseases(:)'
wolffd@0 1483 m = marginal_nodes(engine, i);
wolffd@0 1484 post(i) = m.T(2);
wolffd@0 1485 end
wolffd@0 1486 </pre>
wolffd@0 1487 Junction tree can be quite slow on large QMR models.
wolffd@0 1488 Fortunately, it is possible to exploit properties of the noisy-OR
wolffd@0 1489 function to speed up exact inference using an algorithm called
wolffd@0 1490 <a href="#quickscore">quickscore</a>, discussed below.
wolffd@0 1491
wolffd@0 1492
wolffd@0 1493
wolffd@0 1494
wolffd@0 1495
wolffd@0 1496 <h2><a name="cg_model">Conditional Gaussian models</h2>
wolffd@0 1497
wolffd@0 1498 A conditional Gaussian model is one in which, conditioned on all the discrete
wolffd@0 1499 nodes, the distribution over the remaining (continuous) nodes is
wolffd@0 1500 multivariate Gaussian. This means we can have arcs from discrete (D)
wolffd@0 1501 to continuous (C) nodes, but not vice versa.
wolffd@0 1502 (We <em>are</em> allowed C->D arcs if the continuous nodes are observed,
wolffd@0 1503 as in the <a href="#mixexp">mixture of experts</a> model,
wolffd@0 1504 since this distribution can be represented with a discrete potential.)
wolffd@0 1505 <p>
wolffd@0 1506 We now give an example of a CG model, from
wolffd@0 1507 the paper "Propagation of Probabilities, Means amd
wolffd@0 1508 Variances in Mixed Graphical Association Models", Steffen Lauritzen,
wolffd@0 1509 JASA 87(420):1098--1108, 1992 (reprinted in the book "Probabilistic Networks and Expert
wolffd@0 1510 Systems", R. G. Cowell, A. P. Dawid, S. L. Lauritzen and
wolffd@0 1511 D. J. Spiegelhalter, Springer, 1999.)
wolffd@0 1512
wolffd@0 1513 <h3>Specifying the graph</h3>
wolffd@0 1514
wolffd@0 1515 Consider the model of waste emissions from an incinerator plant shown below.
wolffd@0 1516 We follow the standard convention that shaded nodes are observed,
wolffd@0 1517 clear nodes are hidden.
wolffd@0 1518 We also use the non-standard convention that
wolffd@0 1519 square nodes are discrete (tabular) and round nodes are
wolffd@0 1520 Gaussian.
wolffd@0 1521
wolffd@0 1522 <p>
wolffd@0 1523 <center>
wolffd@0 1524 <IMG SRC="Figures/cg1.gif">
wolffd@0 1525 </center>
wolffd@0 1526 <p>
wolffd@0 1527
wolffd@0 1528 We can create this model as follows.
wolffd@0 1529 <pre>
wolffd@0 1530 F = 1; W = 2; E = 3; B = 4; C = 5; D = 6; Min = 7; Mout = 8; L = 9;
wolffd@0 1531 n = 9;
wolffd@0 1532
wolffd@0 1533 dag = zeros(n);
wolffd@0 1534 dag(F,E)=1;
wolffd@0 1535 dag(W,[E Min D]) = 1;
wolffd@0 1536 dag(E,D)=1;
wolffd@0 1537 dag(B,[C D])=1;
wolffd@0 1538 dag(D,[L Mout])=1;
wolffd@0 1539 dag(Min,Mout)=1;
wolffd@0 1540
wolffd@0 1541 % node sizes - all cts nodes are scalar, all discrete nodes are binary
wolffd@0 1542 ns = ones(1, n);
wolffd@0 1543 dnodes = [F W B];
wolffd@0 1544 cnodes = mysetdiff(1:n, dnodes);
wolffd@0 1545 ns(dnodes) = 2;
wolffd@0 1546
wolffd@0 1547 bnet = mk_bnet(dag, ns, 'discrete', dnodes);
wolffd@0 1548 </pre>
wolffd@0 1549 'dnodes' is a list of the discrete nodes; 'cnodes' is the continuous
wolffd@0 1550 nodes. 'mysetdiff' is a faster version of the built-in 'setdiff'.
wolffd@0 1551 <p>
wolffd@0 1552
wolffd@0 1553
wolffd@0 1554 <h3>Specifying the parameters</h3>
wolffd@0 1555
wolffd@0 1556 The parameters of the discrete nodes can be specified as follows.
wolffd@0 1557 <pre>
wolffd@0 1558 bnet.CPD{B} = tabular_CPD(bnet, B, 'CPT', [0.85 0.15]); % 1=stable, 2=unstable
wolffd@0 1559 bnet.CPD{F} = tabular_CPD(bnet, F, 'CPT', [0.95 0.05]); % 1=intact, 2=defect
wolffd@0 1560 bnet.CPD{W} = tabular_CPD(bnet, W, 'CPT', [2/7 5/7]); % 1=industrial, 2=household
wolffd@0 1561 </pre>
wolffd@0 1562
wolffd@0 1563 <p>
wolffd@0 1564 The parameters of the continuous nodes can be specified as follows.
wolffd@0 1565 <pre>
wolffd@0 1566 bnet.CPD{E} = gaussian_CPD(bnet, E, 'mean', [-3.9 -0.4 -3.2 -0.5], ...
wolffd@0 1567 'cov', [0.00002 0.0001 0.00002 0.0001]);
wolffd@0 1568 bnet.CPD{D} = gaussian_CPD(bnet, D, 'mean', [6.5 6.0 7.5 7.0], ...
wolffd@0 1569 'cov', [0.03 0.04 0.1 0.1], 'weights', [1 1 1 1]);
wolffd@0 1570 bnet.CPD{C} = gaussian_CPD(bnet, C, 'mean', [-2 -1], 'cov', [0.1 0.3]);
wolffd@0 1571 bnet.CPD{L} = gaussian_CPD(bnet, L, 'mean', 3, 'cov', 0.25, 'weights', -0.5);
wolffd@0 1572 bnet.CPD{Min} = gaussian_CPD(bnet, Min, 'mean', [0.5 -0.5], 'cov', [0.01 0.005]);
wolffd@0 1573 bnet.CPD{Mout} = gaussian_CPD(bnet, Mout, 'mean', 0, 'cov', 0.002, 'weights', [1 1]);
wolffd@0 1574 </pre>
wolffd@0 1575
wolffd@0 1576
wolffd@0 1577 <h3><a name="cg_infer">Inference</h3>
wolffd@0 1578
wolffd@0 1579 <!--Let us perform inference in the <a href="#cg_model">waste incinerator example</a>.-->
wolffd@0 1580 First we compute the unconditional marginals.
wolffd@0 1581 <pre>
wolffd@0 1582 engine = jtree_inf_engine(bnet);
wolffd@0 1583 evidence = cell(1,n);
wolffd@0 1584 [engine, ll] = enter_evidence(engine, evidence);
wolffd@0 1585 marg = marginal_nodes(engine, E);
wolffd@0 1586 </pre>
wolffd@0 1587 <!--(Of course, we could use <tt>cond_gauss_inf_engine</tt> instead of jtree.)-->
wolffd@0 1588 'marg' is a structure that contains the fields 'mu' and 'Sigma', which
wolffd@0 1589 contain the mean and (co)variance of the marginal on E.
wolffd@0 1590 In this case, they are both scalars.
wolffd@0 1591 Let us check they match the published figures (to 2 decimal places).
wolffd@0 1592 <!--(We can't expect
wolffd@0 1593 more precision than this in general because I have implemented the algorithm of
wolffd@0 1594 Lauritzen (1992), which can be numerically unstable.)-->
wolffd@0 1595 <pre>
wolffd@0 1596 tol = 1e-2;
wolffd@0 1597 assert(approxeq(marg.mu, -3.25, tol));
wolffd@0 1598 assert(approxeq(sqrt(marg.Sigma), 0.709, tol));
wolffd@0 1599 </pre>
wolffd@0 1600 We can compute the other posteriors similarly.
wolffd@0 1601 Now let us add some evidence.
wolffd@0 1602 <pre>
wolffd@0 1603 evidence = cell(1,n);
wolffd@0 1604 evidence{W} = 1; % industrial
wolffd@0 1605 evidence{L} = 1.1;
wolffd@0 1606 evidence{C} = -0.9;
wolffd@0 1607 [engine, ll] = enter_evidence(engine, evidence);
wolffd@0 1608 </pre>
wolffd@0 1609 Now we find
wolffd@0 1610 <pre>
wolffd@0 1611 marg = marginal_nodes(engine, E);
wolffd@0 1612 assert(approxeq(marg.mu, -3.8983, tol));
wolffd@0 1613 assert(approxeq(sqrt(marg.Sigma), 0.0763, tol));
wolffd@0 1614 </pre>
wolffd@0 1615
wolffd@0 1616
wolffd@0 1617 We can also compute the joint probability on a set of nodes.
wolffd@0 1618 For example, P(D, Mout | evidence) is a 2D Gaussian:
wolffd@0 1619 <pre>
wolffd@0 1620 marg = marginal_nodes(engine, [D Mout])
wolffd@0 1621 marg =
wolffd@0 1622 domain: [6 8]
wolffd@0 1623 mu: [2x1 double]
wolffd@0 1624 Sigma: [2x2 double]
wolffd@0 1625 T: 1.0000
wolffd@0 1626 </pre>
wolffd@0 1627 The mean is
wolffd@0 1628 <pre>
wolffd@0 1629 marg.mu
wolffd@0 1630 ans =
wolffd@0 1631 3.6077
wolffd@0 1632 4.1077
wolffd@0 1633 </pre>
wolffd@0 1634 and the covariance matrix is
wolffd@0 1635 <pre>
wolffd@0 1636 marg.Sigma
wolffd@0 1637 ans =
wolffd@0 1638 0.1062 0.1062
wolffd@0 1639 0.1062 0.1182
wolffd@0 1640 </pre>
wolffd@0 1641 It is easy to visualize this posterior using standard Matlab plotting
wolffd@0 1642 functions, e.g.,
wolffd@0 1643 <pre>
wolffd@0 1644 gaussplot2d(marg.mu, marg.Sigma);
wolffd@0 1645 </pre>
wolffd@0 1646 produces the following picture.
wolffd@0 1647
wolffd@0 1648 <p>
wolffd@0 1649 <center>
wolffd@0 1650 <IMG SRC="Figures/gaussplot.png">
wolffd@0 1651 </center>
wolffd@0 1652 <p>
wolffd@0 1653
wolffd@0 1654
wolffd@0 1655 The T field indicates that the mixing weight of this Gaussian
wolffd@0 1656 component is 1.0.
wolffd@0 1657 If the joint contains discrete and continuous variables, the result
wolffd@0 1658 will be a mixture of Gaussians, e.g.,
wolffd@0 1659 <pre>
wolffd@0 1660 marg = marginal_nodes(engine, [F E])
wolffd@0 1661 domain: [1 3]
wolffd@0 1662 mu: [-3.9000 -0.4003]
wolffd@0 1663 Sigma: [1x1x2 double]
wolffd@0 1664 T: [0.9995 4.7373e-04]
wolffd@0 1665 </pre>
wolffd@0 1666 The interpretation is
wolffd@0 1667 Sigma(i,j,k) = Cov[ E(i) E(j) | F=k ].
wolffd@0 1668 In this case, E is a scalar, so i=j=1; k specifies the mixture component.
wolffd@0 1669 <p>
wolffd@0 1670 We saw in the sprinkler network that BNT sets the effective size of
wolffd@0 1671 observed discrete nodes to 1, since they only have one legal value.
wolffd@0 1672 For continuous nodes, BNT sets their length to 0,
wolffd@0 1673 since they have been reduced to a point.
wolffd@0 1674 For example,
wolffd@0 1675 <pre>
wolffd@0 1676 marg = marginal_nodes(engine, [B C])
wolffd@0 1677 domain: [4 5]
wolffd@0 1678 mu: []
wolffd@0 1679 Sigma: []
wolffd@0 1680 T: [0.0123 0.9877]
wolffd@0 1681 </pre>
wolffd@0 1682 It is simple to post-process the output of marginal_nodes.
wolffd@0 1683 For example, the file BNT/examples/static/cg1 sets the mu term of
wolffd@0 1684 observed nodes to their observed value, and the Sigma term to 0 (since
wolffd@0 1685 observed nodes have no variance).
wolffd@0 1686
wolffd@0 1687 <p>
wolffd@0 1688 Note that the implemented version of the junction tree is numerically
wolffd@0 1689 unstable when using CG potentials
wolffd@0 1690 (which is why, in the example above, we only required our answers to agree with
wolffd@0 1691 the published ones to 2dp.)
wolffd@0 1692 This is why you might want to use <tt>stab_cond_gauss_inf_engine</tt>,
wolffd@0 1693 implemented by Shan Huang. This is described in
wolffd@0 1694
wolffd@0 1695 <ul>
wolffd@0 1696 <li> "Stable Local Computation with Conditional Gaussian Distributions",
wolffd@0 1697 S. Lauritzen and F. Jensen, Tech Report R-99-2014,
wolffd@0 1698 Dept. Math. Sciences, Allborg Univ., 1999.
wolffd@0 1699 </ul>
wolffd@0 1700
wolffd@0 1701 However, even the numerically stable version
wolffd@0 1702 can be computationally intractable if there are many hidden discrete
wolffd@0 1703 nodes, because the number of mixture components grows exponentially e.g., in a
wolffd@0 1704 <a href="usage_dbn.html#lds">switching linear dynamical system</a>.
wolffd@0 1705 In general, one must resort to approximate inference techniques: see
wolffd@0 1706 the discussion on <a href="#engines">inference engines</a> below.
wolffd@0 1707
wolffd@0 1708
wolffd@0 1709 <h2><a name="hybrid">Other hybrid models</h2>
wolffd@0 1710
wolffd@0 1711 When we have C->D arcs, where C is hidden, we need to use
wolffd@0 1712 approximate inference.
wolffd@0 1713 One approach (not implemented in BNT) is described in
wolffd@0 1714 <ul>
wolffd@0 1715 <li> <a
wolffd@0 1716 href="http://www.cs.berkeley.edu/~murphyk/Papers/hybrid_uai99.ps.gz">A
wolffd@0 1717 Variational Approximation for Bayesian Networks with
wolffd@0 1718 Discrete and Continuous Latent Variables</a>,
wolffd@0 1719 K. Murphy, UAI 99.
wolffd@0 1720 </ul>
wolffd@0 1721 Of course, one can always use <a href="#sampling">sampling</a> methods
wolffd@0 1722 for approximate inference in such models.
wolffd@0 1723
wolffd@0 1724
wolffd@0 1725
wolffd@0 1726 <h1><a name="param_learning">Parameter Learning</h1>
wolffd@0 1727
wolffd@0 1728 The parameter estimation routines in BNT can be classified into 4
wolffd@0 1729 types, depending on whether the goal is to compute
wolffd@0 1730 a full (Bayesian) posterior over the parameters or just a point
wolffd@0 1731 estimate (e.g., Maximum Likelihood or Maximum A Posteriori),
wolffd@0 1732 and whether all the variables are fully observed or there is missing
wolffd@0 1733 data/ hidden variables (partial observability).
wolffd@0 1734 <p>
wolffd@0 1735
wolffd@0 1736 <TABLE BORDER>
wolffd@0 1737 <tr>
wolffd@0 1738 <TH></TH>
wolffd@0 1739 <th>Full obs</th>
wolffd@0 1740 <th>Partial obs</th>
wolffd@0 1741 </tr>
wolffd@0 1742 <tr>
wolffd@0 1743 <th>Point</th>
wolffd@0 1744 <td><tt>learn_params</tt></td>
wolffd@0 1745 <td><tt>learn_params_em</tt></td>
wolffd@0 1746 </tr>
wolffd@0 1747 <tr>
wolffd@0 1748 <th>Bayes</th>
wolffd@0 1749 <td><tt>bayes_update_params</tt></td>
wolffd@0 1750 <td>not yet supported</td>
wolffd@0 1751 </tr>
wolffd@0 1752 </table>
wolffd@0 1753
wolffd@0 1754
wolffd@0 1755 <h2><a name="load_data">Loading data from a file</h2>
wolffd@0 1756
wolffd@0 1757 To load numeric data from an ASCII text file called 'dat.txt', where each row is a
wolffd@0 1758 case and columns are separated by white-space, such as
wolffd@0 1759 <pre>
wolffd@0 1760 011979 1626.5 0.0
wolffd@0 1761 021979 1367.0 0.0
wolffd@0 1762 ...
wolffd@0 1763 </pre>
wolffd@0 1764 you can use
wolffd@0 1765 <pre>
wolffd@0 1766 data = load('dat.txt');
wolffd@0 1767 </pre>
wolffd@0 1768 or
wolffd@0 1769 <pre>
wolffd@0 1770 load dat.txt -ascii
wolffd@0 1771 </pre>
wolffd@0 1772 In the latter case, the data is stored in a variable called 'dat' (the
wolffd@0 1773 filename minus the extension).
wolffd@0 1774 Alternatively, suppose the data is stored in a .csv file (has commas
wolffd@0 1775 separating the columns, and contains a header line), such as
wolffd@0 1776 <pre>
wolffd@0 1777 header info goes here
wolffd@0 1778 ORD,011979,1626.5,0.0
wolffd@0 1779 DSM,021979,1367.0,0.0
wolffd@0 1780 ...
wolffd@0 1781 </pre>
wolffd@0 1782 You can load this using
wolffd@0 1783 <pre>
wolffd@0 1784 [a,b,c,d] = textread('dat.txt', '%s %d %f %f', 'delimiter', ',', 'headerlines', 1);
wolffd@0 1785 </pre>
wolffd@0 1786 If your file is not in either of these formats, you can either use Perl to convert
wolffd@0 1787 it to this format, or use the Matlab scanf command.
wolffd@0 1788 Type
wolffd@0 1789 <tt>
wolffd@0 1790 help iofun
wolffd@0 1791 </tt>
wolffd@0 1792 for more information on Matlab's file functions.
wolffd@0 1793 <!--
wolffd@0 1794 <p>
wolffd@0 1795 To load data directly from Excel,
wolffd@0 1796 you should buy the
wolffd@0 1797 <a href="http://www.mathworks.com/products/excellink/">Excel Link</a>.
wolffd@0 1798 To load data directly from a relational database,
wolffd@0 1799 you should buy the
wolffd@0 1800 <a href="http://www.mathworks.com/products/database">Database
wolffd@0 1801 toolbox</a>.
wolffd@0 1802 -->
wolffd@0 1803 <p>
wolffd@0 1804 BNT learning routines require data to be stored in a cell array.
wolffd@0 1805 data{i,m} is the value of node i in case (example) m, i.e., each
wolffd@0 1806 <em>column</em> is a case.
wolffd@0 1807 If node i is not observed in case m (missing value), set
wolffd@0 1808 data{i,m} = [].
wolffd@0 1809 (Not all the learning routines can cope with such missing values, however.)
wolffd@0 1810 In the special case that all the nodes are observed and are
wolffd@0 1811 scalar-valued (as opposed to vector-valued), the data can be
wolffd@0 1812 stored in a matrix (as opposed to a cell-array).
wolffd@0 1813 <p>
wolffd@0 1814 Suppose, as in the <a href="#mixexp">mixture of experts example</a>,
wolffd@0 1815 that we have 3 nodes in the graph: X(1) is the observed input, X(3) is
wolffd@0 1816 the observed output, and X(2) is a hidden (gating) node. We can
wolffd@0 1817 create the dataset as follows.
wolffd@0 1818 <pre>
wolffd@0 1819 data = load('dat.txt');
wolffd@0 1820 ncases = size(data, 1);
wolffd@0 1821 cases = cell(3, ncases);
wolffd@0 1822 cases([1 3], :) = num2cell(data');
wolffd@0 1823 </pre>
wolffd@0 1824 Notice how we transposed the data, to convert rows into columns.
wolffd@0 1825 Also, cases{2,m} = [] for all m, since X(2) is always hidden.
wolffd@0 1826
wolffd@0 1827
wolffd@0 1828 <h2><a name="mle_complete">Maximum likelihood parameter estimation from complete data</h2>
wolffd@0 1829
wolffd@0 1830 As an example, let's generate some data from the sprinkler network, randomize the parameters,
wolffd@0 1831 and then try to recover the original model.
wolffd@0 1832 First we create some training data using forwards sampling.
wolffd@0 1833 <pre>
wolffd@0 1834 samples = cell(N, nsamples);
wolffd@0 1835 for i=1:nsamples
wolffd@0 1836 samples(:,i) = sample_bnet(bnet);
wolffd@0 1837 end
wolffd@0 1838 </pre>
wolffd@0 1839 samples{j,i} contains the value of the j'th node in case i.
wolffd@0 1840 sample_bnet returns a cell array because, in general, each node might
wolffd@0 1841 be a vector of different length.
wolffd@0 1842 In this case, all nodes are discrete (and hence scalars), so we
wolffd@0 1843 could have used a regular array instead (which can be quicker):
wolffd@0 1844 <pre>
wolffd@0 1845 data = cell2num(samples);
wolffd@0 1846 </pre
wolffd@0 1847 So now data(j,i) = samples{j,i}.
wolffd@0 1848 <p>
wolffd@0 1849 Now we create a network with random parameters.
wolffd@0 1850 (The initial values of bnet2 don't matter in this case, since we can find the
wolffd@0 1851 globally optimal MLE independent of where we start.)
wolffd@0 1852 <pre>
wolffd@0 1853 % Make a tabula rasa
wolffd@0 1854 bnet2 = mk_bnet(dag, node_sizes);
wolffd@0 1855 seed = 0;
wolffd@0 1856 rand('state', seed);
wolffd@0 1857 bnet2.CPD{C} = tabular_CPD(bnet2, C);
wolffd@0 1858 bnet2.CPD{R} = tabular_CPD(bnet2, R);
wolffd@0 1859 bnet2.CPD{S} = tabular_CPD(bnet2, S);
wolffd@0 1860 bnet2.CPD{W} = tabular_CPD(bnet2, W);
wolffd@0 1861 </pre>
wolffd@0 1862 Finally, we find the maximum likelihood estimates of the parameters.
wolffd@0 1863 <pre>
wolffd@0 1864 bnet3 = learn_params(bnet2, samples);
wolffd@0 1865 </pre>
wolffd@0 1866 To view the learned parameters, we use a little Matlab hackery.
wolffd@0 1867 <pre>
wolffd@0 1868 CPT3 = cell(1,N);
wolffd@0 1869 for i=1:N
wolffd@0 1870 s=struct(bnet3.CPD{i}); % violate object privacy
wolffd@0 1871 CPT3{i}=s.CPT;
wolffd@0 1872 end
wolffd@0 1873 </pre>
wolffd@0 1874 Here are the parameters learned for node 4.
wolffd@0 1875 <pre>
wolffd@0 1876 dispcpt(CPT3{4})
wolffd@0 1877 1 1 : 1.0000 0.0000
wolffd@0 1878 2 1 : 0.2000 0.8000
wolffd@0 1879 1 2 : 0.2273 0.7727
wolffd@0 1880 2 2 : 0.0000 1.0000
wolffd@0 1881 </pre>
wolffd@0 1882 So we see that the learned parameters are fairly close to the "true"
wolffd@0 1883 ones, which we display below.
wolffd@0 1884 <pre>
wolffd@0 1885 dispcpt(CPT{4})
wolffd@0 1886 1 1 : 1.0000 0.0000
wolffd@0 1887 2 1 : 0.1000 0.9000
wolffd@0 1888 1 2 : 0.1000 0.9000
wolffd@0 1889 2 2 : 0.0100 0.9900
wolffd@0 1890 </pre>
wolffd@0 1891 We can get better results by using a larger training set, or using
wolffd@0 1892 informative priors (see <a href="#prior">below</a>).
wolffd@0 1893
wolffd@0 1894
wolffd@0 1895
wolffd@0 1896 <h2><a name="prior">Parameter priors</h2>
wolffd@0 1897
wolffd@0 1898 Currently, only tabular CPDs can have priors on their parameters.
wolffd@0 1899 The conjugate prior for a multinomial is the Dirichlet.
wolffd@0 1900 (For binary random variables, the multinomial is the same as the
wolffd@0 1901 Bernoulli, and the Dirichlet is the same as the Beta.)
wolffd@0 1902 <p>
wolffd@0 1903 The Dirichlet has a simple interpretation in terms of pseudo counts.
wolffd@0 1904 If we let N_ijk = the num. times X_i=k and Pa_i=j occurs in the
wolffd@0 1905 training set, where Pa_i are the parents of X_i,
wolffd@0 1906 then the maximum likelihood (ML) estimate is
wolffd@0 1907 T_ijk = N_ijk / N_ij (where N_ij = sum_k' N_ijk'), which will be 0 if N_ijk=0.
wolffd@0 1908 To prevent us from declaring that (X_i=k, Pa_i=j) is impossible just because this
wolffd@0 1909 event was not seen in the training set,
wolffd@0 1910 we can pretend we saw value k of X_i, for each value j of Pa_i some number (alpha_ijk)
wolffd@0 1911 of times in the past.
wolffd@0 1912 The MAP (maximum a posterior) estimate is then
wolffd@0 1913 <pre>
wolffd@0 1914 T_ijk = (N_ijk + alpha_ijk) / (N_ij + alpha_ij)
wolffd@0 1915 </pre>
wolffd@0 1916 and is never 0 if all alpha_ijk > 0.
wolffd@0 1917 For example, consider the network A->B, where A is binary and B has 3
wolffd@0 1918 values.
wolffd@0 1919 A uniform prior for B has the form
wolffd@0 1920 <pre>
wolffd@0 1921 B=1 B=2 B=3
wolffd@0 1922 A=1 1 1 1
wolffd@0 1923 A=2 1 1 1
wolffd@0 1924 </pre>
wolffd@0 1925 which can be created using
wolffd@0 1926 <pre>
wolffd@0 1927 tabular_CPD(bnet, i, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif');
wolffd@0 1928 </pre>
wolffd@0 1929 This prior does not satisfy the likelihood equivalence principle,
wolffd@0 1930 which says that <a href="#markov_equiv">Markov equivalent</a> models
wolffd@0 1931 should have the same marginal likelihood.
wolffd@0 1932 A prior that does satisfy this principle is shown below.
wolffd@0 1933 Heckerman (1995) calls this the
wolffd@0 1934 BDeu prior (likelihood equivalent uniform Bayesian Dirichlet).
wolffd@0 1935 <pre>
wolffd@0 1936 B=1 B=2 B=3
wolffd@0 1937 A=1 1/6 1/6 1/6
wolffd@0 1938 A=2 1/6 1/6 1/6
wolffd@0 1939 </pre>
wolffd@0 1940 where we put N/(q*r) in each bin; N is the equivalent sample size,
wolffd@0 1941 r=|A|, q = |B|.
wolffd@0 1942 This can be created as follows
wolffd@0 1943 <pre>
wolffd@0 1944 tabular_CPD(bnet, i, 'prior_type', 'dirichlet', 'dirichlet_type', 'BDeu');
wolffd@0 1945 </pre>
wolffd@0 1946 Here, 1 is the equivalent sample size, and is the strength of the
wolffd@0 1947 prior.
wolffd@0 1948 You can change this using
wolffd@0 1949 <pre>
wolffd@0 1950 tabular_CPD(bnet, i, 'prior_type', 'dirichlet', 'dirichlet_type', ...
wolffd@0 1951 'BDeu', 'dirichlet_weight', 10);
wolffd@0 1952 </pre>
wolffd@0 1953 <!--where counts is an array of pseudo-counts of the same size as the
wolffd@0 1954 CPT.-->
wolffd@0 1955 <!--
wolffd@0 1956 <p>
wolffd@0 1957 When you specify a prior, you should set row i of the CPT to the
wolffd@0 1958 normalized version of row i of the pseudo-count matrix, i.e., to the
wolffd@0 1959 expected values of the parameters. This will ensure that computing the
wolffd@0 1960 marginal likelihood sequentially (see <a
wolffd@0 1961 href="#bayes_learn">below</a>) and in batch form gives the same
wolffd@0 1962 results.
wolffd@0 1963 To do this, proceed as follows.
wolffd@0 1964 <pre>
wolffd@0 1965 tabular_CPD(bnet, i, 'prior', counts, 'CPT', mk_stochastic(counts));
wolffd@0 1966 </pre>
wolffd@0 1967 For a non-informative prior, you can just write
wolffd@0 1968 <pre>
wolffd@0 1969 tabular_CPD(bnet, i, 'prior', 'unif', 'CPT', 'unif');
wolffd@0 1970 </pre>
wolffd@0 1971 -->
wolffd@0 1972
wolffd@0 1973
wolffd@0 1974 <h2><a name="bayes_learn">(Sequential) Bayesian parameter updating from complete data</h2>
wolffd@0 1975
wolffd@0 1976 If we use conjugate priors and have fully observed data, we can
wolffd@0 1977 compute the posterior over the parameters in batch form as follows.
wolffd@0 1978 <pre>
wolffd@0 1979 cases = sample_bnet(bnet, nsamples);
wolffd@0 1980 bnet = bayes_update_params(bnet, cases);
wolffd@0 1981 LL = log_marg_lik_complete(bnet, cases);
wolffd@0 1982 </pre>
wolffd@0 1983 bnet.CPD{i}.prior contains the new Dirichlet pseudocounts,
wolffd@0 1984 and bnet.CPD{i}.CPT is set to the mean of the posterior (the
wolffd@0 1985 normalized counts).
wolffd@0 1986 (Hence if the initial pseudo counts are 0,
wolffd@0 1987 <tt>bayes_update_params</tt> and <tt>learn_params</tt> will give the
wolffd@0 1988 same result.)
wolffd@0 1989
wolffd@0 1990
wolffd@0 1991
wolffd@0 1992
wolffd@0 1993 <p>
wolffd@0 1994 We can compute the same result sequentially (on-line) as follows.
wolffd@0 1995 <pre>
wolffd@0 1996 LL = 0;
wolffd@0 1997 for m=1:nsamples
wolffd@0 1998 LL = LL + log_marg_lik_complete(bnet, cases(:,m));
wolffd@0 1999 bnet = bayes_update_params(bnet, cases(:,m));
wolffd@0 2000 end
wolffd@0 2001 </pre>
wolffd@0 2002
wolffd@0 2003 The file <tt>BNT/examples/static/StructLearn/model_select1</tt> has an example of
wolffd@0 2004 sequential model selection which uses the same idea.
wolffd@0 2005 We generate data from the model A->B
wolffd@0 2006 and compute the posterior prob of all 3 dags on 2 nodes:
wolffd@0 2007 (1) A B, (2) A <- B , (3) A -> B
wolffd@0 2008 Models 2 and 3 are <a href="#markov_equiv">Markov equivalent</a>, and therefore indistinguishable from
wolffd@0 2009 observational data alone, so we expect their posteriors to be the same
wolffd@0 2010 (assuming a prior which satisfies likelihood equivalence).
wolffd@0 2011 If we use random parameters, the "true" model only gets a higher posterior after 2000 trials!
wolffd@0 2012 However, if we make B a noisy NOT gate, the true model "wins" after 12
wolffd@0 2013 trials, as shown below (red = model 1, blue/green (superimposed)
wolffd@0 2014 represents models 2/3).
wolffd@0 2015 <p>
wolffd@0 2016 <img src="Figures/model_select.png">
wolffd@0 2017 <p>
wolffd@0 2018 The use of marginal likelihood for model selection is discussed in
wolffd@0 2019 greater detail in the
wolffd@0 2020 section on <a href="structure_learning">structure learning</a>.
wolffd@0 2021
wolffd@0 2022
wolffd@0 2023
wolffd@0 2024
wolffd@0 2025 <h2><a name="em">Maximum likelihood parameter estimation with missing values</h2>
wolffd@0 2026
wolffd@0 2027 Now we consider learning when some values are not observed.
wolffd@0 2028 Let us randomly hide half the values generated from the water
wolffd@0 2029 sprinkler example.
wolffd@0 2030 <pre>
wolffd@0 2031 samples2 = samples;
wolffd@0 2032 hide = rand(N, nsamples) > 0.5;
wolffd@0 2033 [I,J]=find(hide);
wolffd@0 2034 for k=1:length(I)
wolffd@0 2035 samples2{I(k), J(k)} = [];
wolffd@0 2036 end
wolffd@0 2037 </pre>
wolffd@0 2038 samples2{i,l} is the value of node i in training case l, or [] if unobserved.
wolffd@0 2039 <p>
wolffd@0 2040 Now we will compute the MLEs using the EM algorithm.
wolffd@0 2041 We need to use an inference algorithm to compute the expected
wolffd@0 2042 sufficient statistics in the E step; the M (maximization) step is as
wolffd@0 2043 above.
wolffd@0 2044 <pre>
wolffd@0 2045 engine2 = jtree_inf_engine(bnet2);
wolffd@0 2046 max_iter = 10;
wolffd@0 2047 [bnet4, LLtrace] = learn_params_em(engine2, samples2, max_iter);
wolffd@0 2048 </pre>
wolffd@0 2049 LLtrace(i) is the log-likelihood at iteration i. We can plot this as
wolffd@0 2050 follows:
wolffd@0 2051 <pre>
wolffd@0 2052 plot(LLtrace, 'x-')
wolffd@0 2053 </pre>
wolffd@0 2054 Let's display the results after 10 iterations of EM.
wolffd@0 2055 <pre>
wolffd@0 2056 celldisp(CPT4)
wolffd@0 2057 CPT4{1} =
wolffd@0 2058 0.6616
wolffd@0 2059 0.3384
wolffd@0 2060 CPT4{2} =
wolffd@0 2061 0.6510 0.3490
wolffd@0 2062 0.8751 0.1249
wolffd@0 2063 CPT4{3} =
wolffd@0 2064 0.8366 0.1634
wolffd@0 2065 0.0197 0.9803
wolffd@0 2066 CPT4{4} =
wolffd@0 2067 (:,:,1) =
wolffd@0 2068 0.8276 0.0546
wolffd@0 2069 0.5452 0.1658
wolffd@0 2070 (:,:,2) =
wolffd@0 2071 0.1724 0.9454
wolffd@0 2072 0.4548 0.8342
wolffd@0 2073 </pre>
wolffd@0 2074 We can get improved performance by using one or more of the following
wolffd@0 2075 methods:
wolffd@0 2076 <ul>
wolffd@0 2077 <li> Increasing the size of the training set.
wolffd@0 2078 <li> Decreasing the amount of hidden data.
wolffd@0 2079 <li> Running EM for longer.
wolffd@0 2080 <li> Using informative priors.
wolffd@0 2081 <li> Initialising EM from multiple starting points.
wolffd@0 2082 </ul>
wolffd@0 2083
wolffd@0 2084 Click <a href="#gaussian">here</a> for a discussion of learning
wolffd@0 2085 Gaussians, which can cause numerical problems.
wolffd@0 2086 <p>
wolffd@0 2087 For a more complete example of learning with EM,
wolffd@0 2088 see the script BNT/examples/static/learn1.m.
wolffd@0 2089
wolffd@0 2090 <h2><a name="tying">Parameter tying</h2>
wolffd@0 2091
wolffd@0 2092 In networks with repeated structure (e.g., chains and grids), it is
wolffd@0 2093 common to assume that the parameters are the same at every node. This
wolffd@0 2094 is called parameter tying, and reduces the amount of data needed for
wolffd@0 2095 learning.
wolffd@0 2096 <p>
wolffd@0 2097 When we have tied parameters, there is no longer a one-to-one
wolffd@0 2098 correspondence between nodes and CPDs.
wolffd@0 2099 Rather, each CPD species the parameters for a whole equivalence class
wolffd@0 2100 of nodes.
wolffd@0 2101 It is easiest to see this by example.
wolffd@0 2102 Consider the following <a href="usage_dbn.html#hmm">hidden Markov
wolffd@0 2103 model (HMM)</a>
wolffd@0 2104 <p>
wolffd@0 2105 <img src="Figures/hmm3.gif">
wolffd@0 2106 <p>
wolffd@0 2107 <!--
wolffd@0 2108 We can create this graph structure, assuming we have T time-slices,
wolffd@0 2109 as follows.
wolffd@0 2110 (We number the nodes as shown in the figure, but we could equally well
wolffd@0 2111 number the hidden nodes 1:T, and the observed nodes T+1:2T.)
wolffd@0 2112 <pre>
wolffd@0 2113 N = 2*T;
wolffd@0 2114 dag = zeros(N);
wolffd@0 2115 hnodes = 1:2:2*T;
wolffd@0 2116 for i=1:T-1
wolffd@0 2117 dag(hnodes(i), hnodes(i+1))=1;
wolffd@0 2118 end
wolffd@0 2119 onodes = 2:2:2*T;
wolffd@0 2120 for i=1:T
wolffd@0 2121 dag(hnodes(i), onodes(i)) = 1;
wolffd@0 2122 end
wolffd@0 2123 </pre>
wolffd@0 2124 <p>
wolffd@0 2125 The hidden nodes are always discrete, and have Q possible values each,
wolffd@0 2126 but the observed nodes can be discrete or continuous, and have O possible values/length.
wolffd@0 2127 <pre>
wolffd@0 2128 if cts_obs
wolffd@0 2129 dnodes = hnodes;
wolffd@0 2130 else
wolffd@0 2131 dnodes = 1:N;
wolffd@0 2132 end
wolffd@0 2133 ns = ones(1,N);
wolffd@0 2134 ns(hnodes) = Q;
wolffd@0 2135 ns(onodes) = O;
wolffd@0 2136 </pre>
wolffd@0 2137 -->
wolffd@0 2138 When HMMs are used for semi-infinite processes like speech recognition,
wolffd@0 2139 we assume the transition matrix
wolffd@0 2140 P(H(t+1)|H(t)) is the same for all t; this is called a time-invariant
wolffd@0 2141 or homogenous Markov chain.
wolffd@0 2142 Hence hidden nodes 2, 3, ..., T
wolffd@0 2143 are all in the same equivalence class, say class Hclass.
wolffd@0 2144 Similarly, the observation matrix P(O(t)|H(t)) is assumed to be the
wolffd@0 2145 same for all t, so the observed nodes are all in the same equivalence
wolffd@0 2146 class, say class Oclass.
wolffd@0 2147 Finally, the prior term P(H(1)) is in a class all by itself, say class
wolffd@0 2148 H1class.
wolffd@0 2149 This is illustrated below, where we explicitly represent the
wolffd@0 2150 parameters as random variables (dotted nodes).
wolffd@0 2151 <p>
wolffd@0 2152 <img src="Figures/hmm4_params.gif">
wolffd@0 2153 <p>
wolffd@0 2154 In BNT, we cannot represent parameters as random variables (nodes).
wolffd@0 2155 Instead, we "hide" the
wolffd@0 2156 parameters inside one CPD for each equivalence class,
wolffd@0 2157 and then specify that the other CPDs should share these parameters, as
wolffd@0 2158 follows.
wolffd@0 2159 <pre>
wolffd@0 2160 hnodes = 1:2:2*T;
wolffd@0 2161 onodes = 2:2:2*T;
wolffd@0 2162 H1class = 1; Hclass = 2; Oclass = 3;
wolffd@0 2163 eclass = ones(1,N);
wolffd@0 2164 eclass(hnodes(2:end)) = Hclass;
wolffd@0 2165 eclass(hnodes(1)) = H1class;
wolffd@0 2166 eclass(onodes) = Oclass;
wolffd@0 2167 % create dag and ns in the usual way
wolffd@0 2168 bnet = mk_bnet(dag, ns, 'discrete', dnodes, 'equiv_class', eclass);
wolffd@0 2169 </pre>
wolffd@0 2170 Finally, we define the parameters for each equivalence class:
wolffd@0 2171 <pre>
wolffd@0 2172 bnet.CPD{H1class} = tabular_CPD(bnet, hnodes(1)); % prior
wolffd@0 2173 bnet.CPD{Hclass} = tabular_CPD(bnet, hnodes(2)); % transition matrix
wolffd@0 2174 if cts_obs
wolffd@0 2175 bnet.CPD{Oclass} = gaussian_CPD(bnet, onodes(1));
wolffd@0 2176 else
wolffd@0 2177 bnet.CPD{Oclass} = tabular_CPD(bnet, onodes(1));
wolffd@0 2178 end
wolffd@0 2179 </pre>
wolffd@0 2180 In general, if bnet.CPD{e} = xxx_CPD(bnet, j), then j should be a
wolffd@0 2181 member of e's equivalence class; that is, it is not always the case
wolffd@0 2182 that e == j. You can use bnet.rep_of_eclass(e) to return the
wolffd@0 2183 representative of equivalence class e.
wolffd@0 2184 BNT will look up the parents of j to determine the size
wolffd@0 2185 of the CPT to use. It assumes that this is the same for all members of
wolffd@0 2186 the equivalence class.
wolffd@0 2187 Click <a href="param_tieing.html">here</a> for
wolffd@0 2188 a more complex example of parameter tying.
wolffd@0 2189 <p>
wolffd@0 2190 Note:
wolffd@0 2191 Normally one would define an HMM as a
wolffd@0 2192 <a href = "usage_dbn.html">Dynamic Bayes Net</a>
wolffd@0 2193 (see the function BNT/examples/dynamic/mk_chmm.m).
wolffd@0 2194 However, one can define an HMM as a static BN using the function
wolffd@0 2195 BNT/examples/static/Models/mk_hmm_bnet.m.
wolffd@0 2196
wolffd@0 2197
wolffd@0 2198
wolffd@0 2199 <h1><a name="structure_learning">Structure learning</h1>
wolffd@0 2200
wolffd@0 2201 Update (9/29/03):
wolffd@0 2202 Phillipe LeRay is developing some additional structure learning code
wolffd@0 2203 on top of BNT. Click
wolffd@0 2204 <a href="http://banquiseasi.insa-rouen.fr/projects/bnt-slp/">
wolffd@0 2205 here</a>
wolffd@0 2206 for details.
wolffd@0 2207
wolffd@0 2208 <p>
wolffd@0 2209
wolffd@0 2210 There are two very different approaches to structure learning:
wolffd@0 2211 constraint-based and search-and-score.
wolffd@0 2212 In the <a href="#constraint">constraint-based approach</a>,
wolffd@0 2213 we start with a fully connected graph, and remove edges if certain
wolffd@0 2214 conditional independencies are measured in the data.
wolffd@0 2215 This has the disadvantage that repeated independence tests lose
wolffd@0 2216 statistical power.
wolffd@0 2217 <p>
wolffd@0 2218 In the more popular search-and-score approach,
wolffd@0 2219 we perform a search through the space of possible DAGs, and either
wolffd@0 2220 return the best one found (a point estimate), or return a sample of the
wolffd@0 2221 models found (an approximation to the Bayesian posterior).
wolffd@0 2222 <p>
wolffd@0 2223 The number of DAGs as a function of the number of
wolffd@0 2224 nodes, G(n), is super-exponential in n,
wolffd@0 2225 and is given by the following recurrence
wolffd@0 2226 <!--(where R(i)=G(n)):-->
wolffd@0 2227 <p>
wolffd@0 2228 <center>
wolffd@0 2229 <IMG SRC="numDAGsEqn2.png">
wolffd@0 2230 </center>
wolffd@0 2231 <p>
wolffd@0 2232 The first few values
wolffd@0 2233 are shown below.
wolffd@0 2234
wolffd@0 2235 <table>
wolffd@0 2236 <tr> <th>n</th> <th align=left>G(n)</th> </tr>
wolffd@0 2237 <tr> <td>1</td> <td>1</td> </tr>
wolffd@0 2238 <tr> <td>2</td> <td>3</td> </tr>
wolffd@0 2239 <tr> <td>3</td> <td>25</td> </tr>
wolffd@0 2240 <tr> <td>4</td> <td>543</td> </tr>
wolffd@0 2241 <tr> <td>5</td> <td>29,281</td> </tr>
wolffd@0 2242 <tr> <td>6</td> <td>3,781,503</td> </tr>
wolffd@0 2243 <tr> <td>7</td> <td>1.1 x 10^9</td> </tr>
wolffd@0 2244 <tr> <td>8</td> <td>7.8 x 10^11</td> </tr>
wolffd@0 2245 <tr> <td>9</td> <td>1.2 x 10^15</td> </tr>
wolffd@0 2246 <tr> <td>10</td> <td>4.2 x 10^18</td> </tr>
wolffd@0 2247 </table>
wolffd@0 2248
wolffd@0 2249 Since the number of DAGs is super-exponential in the number of nodes,
wolffd@0 2250 we cannot exhaustively search the space, so we either use a local
wolffd@0 2251 search algorithm (e.g., greedy hill climbining, perhaps with multiple
wolffd@0 2252 restarts) or a global search algorithm (e.g., Markov Chain Monte
wolffd@0 2253 Carlo).
wolffd@0 2254 <p>
wolffd@0 2255 If we know a total ordering on the nodes,
wolffd@0 2256 finding the best structure amounts to picking the best set of parents
wolffd@0 2257 for each node independently.
wolffd@0 2258 This is what the K2 algorithm does.
wolffd@0 2259 If the ordering is unknown, we can search over orderings,
wolffd@0 2260 which is more efficient than searching over DAGs (Koller and Friedman, 2000).
wolffd@0 2261 <p>
wolffd@0 2262 In addition to the search procedure, we must specify the scoring
wolffd@0 2263 function. There are two popular choices. The Bayesian score integrates
wolffd@0 2264 out the parameters, i.e., it is the marginal likelihood of the model.
wolffd@0 2265 The BIC (Bayesian Information Criterion) is defined as
wolffd@0 2266 log P(D|theta_hat) - 0.5*d*log(N), where D is the data, theta_hat is
wolffd@0 2267 the ML estimate of the parameters, d is the number of parameters, and
wolffd@0 2268 N is the number of data cases.
wolffd@0 2269 The BIC method has the advantage of not requiring a prior.
wolffd@0 2270 <p>
wolffd@0 2271 BIC can be derived as a large sample
wolffd@0 2272 approximation to the marginal likelihood.
wolffd@0 2273 (It is also equal to the Minimum Description Length of a model.)
wolffd@0 2274 However, in practice, the sample size does not need to be very large
wolffd@0 2275 for the approximation to be good.
wolffd@0 2276 For example, in the figure below, we plot the ratio between the log marginal likelihood
wolffd@0 2277 and the BIC score against data-set size; we see that the ratio rapidly
wolffd@0 2278 approaches 1, especially for non-informative priors.
wolffd@0 2279 (This plot was generated by the file BNT/examples/static/bic1.m. It
wolffd@0 2280 uses the water sprinkler BN with BDeu Dirichlet priors with different
wolffd@0 2281 equivalent sample sizes.)
wolffd@0 2282
wolffd@0 2283 <p>
wolffd@0 2284 <center>
wolffd@0 2285 <IMG SRC="Figures/bic.png">
wolffd@0 2286 </center>
wolffd@0 2287 <p>
wolffd@0 2288
wolffd@0 2289 <p>
wolffd@0 2290 As with parameter learning, handling missing data/ hidden variables is
wolffd@0 2291 much harder than the fully observed case.
wolffd@0 2292 The structure learning routines in BNT can therefore be classified into 4
wolffd@0 2293 types, analogously to the parameter learning case.
wolffd@0 2294 <p>
wolffd@0 2295
wolffd@0 2296 <TABLE BORDER>
wolffd@0 2297 <tr>
wolffd@0 2298 <TH></TH>
wolffd@0 2299 <th>Full obs</th>
wolffd@0 2300 <th>Partial obs</th>
wolffd@0 2301 </tr>
wolffd@0 2302 <tr>
wolffd@0 2303 <th>Point</th>
wolffd@0 2304 <td><tt>learn_struct_K2</tt> <br>
wolffd@0 2305 <!-- <tt>learn_struct_hill_climb</tt></td> -->
wolffd@0 2306 <td><tt>not yet supported</tt></td>
wolffd@0 2307 </tr>
wolffd@0 2308 <tr>
wolffd@0 2309 <th>Bayes</th>
wolffd@0 2310 <td><tt>learn_struct_mcmc</tt></td>
wolffd@0 2311 <td>not yet supported</td>
wolffd@0 2312 </tr>
wolffd@0 2313 </table>
wolffd@0 2314
wolffd@0 2315
wolffd@0 2316 <h2><a name="markov_equiv">Markov equivalence</h2>
wolffd@0 2317
wolffd@0 2318 If two DAGs encode the same conditional independencies, they are
wolffd@0 2319 called Markov equivalent. The set of all DAGs can be paritioned into
wolffd@0 2320 Markov equivalence classes. Graphs within the same class can
wolffd@0 2321 have
wolffd@0 2322 the direction of some of their arcs reversed without changing any of
wolffd@0 2323 the CI relationships.
wolffd@0 2324 Each class can be represented by a PDAG
wolffd@0 2325 (partially directed acyclic graph) called an essential graph or
wolffd@0 2326 pattern. This specifies which edges must be oriented in a certain
wolffd@0 2327 direction, and which may be reversed.
wolffd@0 2328
wolffd@0 2329 <p>
wolffd@0 2330 When learning graph structure from observational data,
wolffd@0 2331 the best one can hope to do is to identify the model up to Markov
wolffd@0 2332 equivalence. To distinguish amongst graphs within the same equivalence
wolffd@0 2333 class, one needs interventional data: see the discussion on <a
wolffd@0 2334 href="#active">active learning</a> below.
wolffd@0 2335
wolffd@0 2336
wolffd@0 2337
wolffd@0 2338 <h2><a name="enumerate">Exhaustive search</h2>
wolffd@0 2339
wolffd@0 2340 The brute-force approach to structure learning is to enumerate all
wolffd@0 2341 possible DAGs, and score each one. This provides a "gold standard"
wolffd@0 2342 with which to compare other algorithms. We can do this as follows.
wolffd@0 2343 <pre>
wolffd@0 2344 dags = mk_all_dags(N);
wolffd@0 2345 score = score_dags(data, ns, dags);
wolffd@0 2346 </pre>
wolffd@0 2347 where data(i,m) is the value of node i in case m,
wolffd@0 2348 and ns(i) is the size of node i.
wolffd@0 2349 If the DAGs have a lot of families in common, we can cache the sufficient statistics,
wolffd@0 2350 making this potentially more efficient than scoring the DAGs one at a time.
wolffd@0 2351 (Caching is not currently implemented, however.)
wolffd@0 2352 <p>
wolffd@0 2353 By default, we use the Bayesian scoring metric, and assume CPDs are
wolffd@0 2354 represented by tables with BDeu(1) priors.
wolffd@0 2355 We can override these defaults as follows.
wolffd@0 2356 If we want to use uniform priors, we can say
wolffd@0 2357 <pre>
wolffd@0 2358 params = cell(1,N);
wolffd@0 2359 for i=1:N
wolffd@0 2360 params{i} = {'prior', 'unif'};
wolffd@0 2361 end
wolffd@0 2362 score = score_dags(data, ns, dags, 'params', params);
wolffd@0 2363 </pre>
wolffd@0 2364 params{i} is a cell-array, containing optional arguments that are
wolffd@0 2365 passed to the constructor for CPD i.
wolffd@0 2366 <p>
wolffd@0 2367 Now suppose we want to use different node types, e.g.,
wolffd@0 2368 Suppose nodes 1 and 2 are Gaussian, and nodes 3 and 4 softmax (both
wolffd@0 2369 these CPDs can support discrete and continuous parents, which is
wolffd@0 2370 necessary since all other nodes will be considered as parents).
wolffd@0 2371 The Bayesian scoring metric currently only works for tabular CPDs, so
wolffd@0 2372 we will use BIC:
wolffd@0 2373 <pre>
wolffd@0 2374 score = score_dags(data, ns, dags, 'discrete', [3 4], 'params', [],
wolffd@0 2375 'type', {'gaussian', 'gaussian', 'softmax', softmax'}, 'scoring_fn', 'bic')
wolffd@0 2376 </pre>
wolffd@0 2377 In practice, one can't enumerate all possible DAGs for N > 5,
wolffd@0 2378 but one can evaluate any reasonably-sized set of hypotheses in this
wolffd@0 2379 way (e.g., nearest neighbors of your current best guess).
wolffd@0 2380 Think of this as "computer assisted model refinement" as opposed to de
wolffd@0 2381 novo learning.
wolffd@0 2382
wolffd@0 2383
wolffd@0 2384 <h2><a name="K2">K2</h2>
wolffd@0 2385
wolffd@0 2386 The K2 algorithm (Cooper and Herskovits, 1992) is a greedy search algorithm that works as follows.
wolffd@0 2387 Initially each node has no parents. It then adds incrementally that parent whose addition most
wolffd@0 2388 increases the score of the resulting structure. When the addition of no single
wolffd@0 2389 parent can increase the score, it stops adding parents to the node.
wolffd@0 2390 Since we are using a fixed ordering, we do not need to check for
wolffd@0 2391 cycles, and can choose the parents for each node independently.
wolffd@0 2392 <p>
wolffd@0 2393 The original paper used the Bayesian scoring
wolffd@0 2394 metric with tabular CPDs and Dirichlet priors.
wolffd@0 2395 BNT generalizes this to allow any kind of CPD, and either the Bayesian
wolffd@0 2396 scoring metric or BIC, as in the example <a href="#enumerate">above</a>.
wolffd@0 2397 In addition, you can specify
wolffd@0 2398 an optional upper bound on the number of parents for each node.
wolffd@0 2399 The file BNT/examples/static/k2demo1.m gives an example of how to use K2.
wolffd@0 2400 We use the water sprinkler network and sample 100 cases from it as before.
wolffd@0 2401 Then we see how much data it takes to recover the generating structure:
wolffd@0 2402 <pre>
wolffd@0 2403 order = [C S R W];
wolffd@0 2404 max_fan_in = 2;
wolffd@0 2405 sz = 5:5:100;
wolffd@0 2406 for i=1:length(sz)
wolffd@0 2407 dag2 = learn_struct_K2(data(:,1:sz(i)), node_sizes, order, 'max_fan_in', max_fan_in);
wolffd@0 2408 correct(i) = isequal(dag, dag2);
wolffd@0 2409 end
wolffd@0 2410 </pre>
wolffd@0 2411 Here are the results.
wolffd@0 2412 <pre>
wolffd@0 2413 correct =
wolffd@0 2414 Columns 1 through 12
wolffd@0 2415 0 0 0 0 0 0 0 1 0 1 1 1
wolffd@0 2416 Columns 13 through 20
wolffd@0 2417 1 1 1 1 1 1 1 1
wolffd@0 2418 </pre>
wolffd@0 2419 So we see it takes about sz(10)=50 cases. (BIC behaves similarly,
wolffd@0 2420 showing that the prior doesn't matter too much.)
wolffd@0 2421 In general, we cannot hope to recover the "true" generating structure,
wolffd@0 2422 only one that is in its <a href="#markov_equiv">Markov equivalence
wolffd@0 2423 class</a>.
wolffd@0 2424
wolffd@0 2425
wolffd@0 2426 <h2><a name="hill_climb">Hill-climbing</h2>
wolffd@0 2427
wolffd@0 2428 Hill-climbing starts at a specific point in space,
wolffd@0 2429 considers all nearest neighbors, and moves to the neighbor
wolffd@0 2430 that has the highest score; if no neighbors have higher
wolffd@0 2431 score than the current point (i.e., we have reached a local maximum),
wolffd@0 2432 the algorithm stops. One can then restart in another part of the space.
wolffd@0 2433 <p>
wolffd@0 2434 A common definition of "neighbor" is all graphs that can be
wolffd@0 2435 generated from the current graph by adding, deleting or reversing a
wolffd@0 2436 single arc, subject to the acyclicity constraint.
wolffd@0 2437 Other neighborhoods are possible: see
wolffd@0 2438 <a href="http://research.microsoft.com/~dmax/publications/jmlr02.pdf">
wolffd@0 2439 Optimal Structure Identification with Greedy Search</a>, Max
wolffd@0 2440 Chickering, JMLR 2002.
wolffd@0 2441
wolffd@0 2442 <!--
wolffd@0 2443 Note: This algorithm is currently (Feb '02) being implemented by Qian
wolffd@0 2444 Diao.
wolffd@0 2445 -->
wolffd@0 2446
wolffd@0 2447
wolffd@0 2448 <h2><a name="mcmc">MCMC</h2>
wolffd@0 2449
wolffd@0 2450 We can use a Markov Chain Monte Carlo (MCMC) algorithm called
wolffd@0 2451 Metropolis-Hastings (MH) to search the space of all
wolffd@0 2452 DAGs.
wolffd@0 2453 The standard proposal distribution is to consider moving to all
wolffd@0 2454 nearest neighbors in the sense defined <a href="#hill_climb">above</a>.
wolffd@0 2455 <p>
wolffd@0 2456 The function can be called
wolffd@0 2457 as in the following example.
wolffd@0 2458 <pre>
wolffd@0 2459 [sampled_graphs, accept_ratio] = learn_struct_mcmc(data, ns, 'nsamples', 100, 'burnin', 10);
wolffd@0 2460 </pre>
wolffd@0 2461 We can convert our set of sampled graphs to a histogram
wolffd@0 2462 (empirical posterior over all the DAGs) thus
wolffd@0 2463 <pre>
wolffd@0 2464 all_dags = mk_all_dags(N);
wolffd@0 2465 mcmc_post = mcmc_sample_to_hist(sampled_graphs, all_dags);
wolffd@0 2466 </pre>
wolffd@0 2467 To see how well this performs, let us compute the exact posterior exhaustively.
wolffd@0 2468 <p>
wolffd@0 2469 <pre>
wolffd@0 2470 score = score_dags(data, ns, all_dags);
wolffd@0 2471 post = normalise(exp(score)); % assuming uniform structural prior
wolffd@0 2472 </pre>
wolffd@0 2473 We plot the results below.
wolffd@0 2474 (The data set was 100 samples drawn from a random 4 node bnet; see the
wolffd@0 2475 file BNT/examples/static/mcmc1.)
wolffd@0 2476 <pre>
wolffd@0 2477 subplot(2,1,1)
wolffd@0 2478 bar(post)
wolffd@0 2479 subplot(2,1,2)
wolffd@0 2480 bar(mcmc_post)
wolffd@0 2481 </pre>
wolffd@0 2482 <img src="Figures/mcmc_post.jpg" width="800" height="500">
wolffd@0 2483 <p>
wolffd@0 2484 We can also plot the acceptance ratio versus number of MCMC steps,
wolffd@0 2485 as a crude convergence diagnostic.
wolffd@0 2486 <pre>
wolffd@0 2487 clf
wolffd@0 2488 plot(accept_ratio)
wolffd@0 2489 </pre>
wolffd@0 2490 <img src="Figures/mcmc_accept.jpg" width="800" height="300">
wolffd@0 2491 <p>
wolffd@0 2492 Even though the number of samples needed by MCMC is theoretically
wolffd@0 2493 polynomial (not exponential) in the dimensionality of the search space, in practice it has been
wolffd@0 2494 found that MCMC does not converge in reasonable time for graphs with
wolffd@0 2495 more than about 10 nodes.
wolffd@0 2496
wolffd@0 2497
wolffd@0 2498
wolffd@0 2499
wolffd@0 2500 <h2><a name="active">Active structure learning</h2>
wolffd@0 2501
wolffd@0 2502 As was mentioned <a href="#markov_equiv">above</a>,
wolffd@0 2503 one can only learn a DAG up to Markov equivalence, even given infinite data.
wolffd@0 2504 If one is interested in learning the structure of a causal network,
wolffd@0 2505 one needs interventional data.
wolffd@0 2506 (By "intervention" we mean forcing a node to take on a specific value,
wolffd@0 2507 thereby effectively severing its incoming arcs.)
wolffd@0 2508 <p>
wolffd@0 2509 Most of the scoring functions accept an optional argument
wolffd@0 2510 that specifies whether a node was observed to have a certain value, or
wolffd@0 2511 was forced to have that value: we set clamped(i,m)=1 if node i was
wolffd@0 2512 forced in training case m. e.g., see the file
wolffd@0 2513 BNT/examples/static/cooper_yoo.
wolffd@0 2514 <p>
wolffd@0 2515 An interesting question is to decide which interventions to perform
wolffd@0 2516 (c.f., design of experiments). For details, see the following tech
wolffd@0 2517 report
wolffd@0 2518 <ul>
wolffd@0 2519 <li> <a href = "../../Papers/alearn.ps.gz">
wolffd@0 2520 Active learning of causal Bayes net structure</a>, Kevin Murphy, March
wolffd@0 2521 2001.
wolffd@0 2522 </ul>
wolffd@0 2523
wolffd@0 2524
wolffd@0 2525 <h2><a name="struct_em">Structural EM</h2>
wolffd@0 2526
wolffd@0 2527 Computing the Bayesian score when there is partial observability is
wolffd@0 2528 computationally challenging, because the parameter posterior becomes
wolffd@0 2529 multimodal (the hidden nodes induce a mixture distribution).
wolffd@0 2530 One therefore needs to use approximations such as BIC.
wolffd@0 2531 Unfortunately, search algorithms are still expensive, because we need
wolffd@0 2532 to run EM at each step to compute the MLE, which is needed to compute
wolffd@0 2533 the score of each model. An alternative approach is
wolffd@0 2534 to do the local search steps inside of the M step of EM, which is more
wolffd@0 2535 efficient since the data has been "filled in" - this is
wolffd@0 2536 called the structural EM algorithm (Friedman 1997), and provably
wolffd@0 2537 converges to a local maximum of the BIC score.
wolffd@0 2538 <p>
wolffd@0 2539 Wei Hu has implemented SEM for discrete nodes.
wolffd@0 2540 You can download his package from
wolffd@0 2541 <a href="../SEM.zip">here</a>.
wolffd@0 2542 Please address all questions about this code to
wolffd@0 2543 wei.hu@intel.com.
wolffd@0 2544 See also <a href="#phl">Phl's implementation of SEM</a>.
wolffd@0 2545
wolffd@0 2546 <!--
wolffd@0 2547 <h2><a name="reveal">REVEAL algorithm</h2>
wolffd@0 2548
wolffd@0 2549 A simple way to learn the structure of a fully observed, discrete,
wolffd@0 2550 factored DBN from a time series is described <a
wolffd@0 2551 href="usage_dbn.html#struct_learn">here</a>.
wolffd@0 2552 -->
wolffd@0 2553
wolffd@0 2554
wolffd@0 2555 <h2><a name="graphdraw">Visualizing the graph</h2>
wolffd@0 2556
wolffd@0 2557 Click <a href="graphviz.html">here</a> for more information
wolffd@0 2558 on graph visualization.
wolffd@0 2559
wolffd@0 2560 <h2><a name = "constraint">Constraint-based methods</h2>
wolffd@0 2561
wolffd@0 2562 The IC algorithm (Pearl and Verma, 1991),
wolffd@0 2563 and the faster, but otherwise equivalent, PC algorithm (Spirtes, Glymour, and Scheines 1993),
wolffd@0 2564 computes many conditional independence tests,
wolffd@0 2565 and combines these constraints into a
wolffd@0 2566 PDAG to represent the whole
wolffd@0 2567 <a href="#markov_equiv">Markov equivalence class</a>.
wolffd@0 2568 <p>
wolffd@0 2569 IC*/FCI extend IC/PC to handle latent variables: see <a href="#ic_star">below</a>.
wolffd@0 2570 (IC stands for inductive causation; PC stands for Peter and Clark,
wolffd@0 2571 the first names of Spirtes and Glymour; FCI stands for fast causal
wolffd@0 2572 inference.
wolffd@0 2573 What we, following Pearl (2000), call IC* was called
wolffd@0 2574 IC in the original Pearl and Verma paper.)
wolffd@0 2575 For details, see
wolffd@0 2576 <ul>
wolffd@0 2577 <li>
wolffd@0 2578 <a href="http://hss.cmu.edu/html/departments/philosophy/TETRAD/tetrad.html">Causation,
wolffd@0 2579 Prediction, and Search</a>, Spirtes, Glymour and
wolffd@0 2580 Scheines (SGS), 2001 (2nd edition), MIT Press.
wolffd@0 2581 <li>
wolffd@0 2582 <a href="http://bayes.cs.ucla.edu/BOOK-2K/index.html">Causality: Models, Reasoning and Inference</a>, J. Pearl,
wolffd@0 2583 2000, Cambridge University Press.
wolffd@0 2584 </ul>
wolffd@0 2585
wolffd@0 2586 <p>
wolffd@0 2587
wolffd@0 2588 The PC algorithm takes as arguments a function f, the number of nodes N,
wolffd@0 2589 the maximum fan in K, and additional arguments A which are passed to f.
wolffd@0 2590 The function f(X,Y,S,A) returns 1 if X is conditionally independent of Y given S, and 0
wolffd@0 2591 otherwise.
wolffd@0 2592 For example, suppose we cheat by
wolffd@0 2593 passing in a CI "oracle" which has access to the true DAG; the oracle
wolffd@0 2594 tests for d-separation in this DAG, i.e.,
wolffd@0 2595 f(X,Y,S) calls dsep(X,Y,S,dag). We can to this as follows.
wolffd@0 2596 <pre>
wolffd@0 2597 pdag = learn_struct_pdag_pc('dsep', N, max_fan_in, dag);
wolffd@0 2598 </pre>
wolffd@0 2599 pdag(i,j) = -1 if there is definitely an i->j arc,
wolffd@0 2600 and pdag(i,j) = 1 if there is either an i->j or and i<-j arc.
wolffd@0 2601 <p>
wolffd@0 2602 Applied to the sprinkler network, this returns
wolffd@0 2603 <pre>
wolffd@0 2604 pdag =
wolffd@0 2605 0 1 1 0
wolffd@0 2606 1 0 0 -1
wolffd@0 2607 1 0 0 -1
wolffd@0 2608 0 0 0 0
wolffd@0 2609 </pre>
wolffd@0 2610 So as expected, we see that the V-structure at the W node is uniquely identified,
wolffd@0 2611 but the other arcs have ambiguous orientation.
wolffd@0 2612 <p>
wolffd@0 2613 We now give an example from p141 (1st edn) / p103 (2nd end) of the SGS
wolffd@0 2614 book.
wolffd@0 2615 This example concerns the female orgasm.
wolffd@0 2616 We are given a correlation matrix C between 7 measured factors (such
wolffd@0 2617 as subjective experiences of coital and masturbatory experiences),
wolffd@0 2618 derived from 281 samples, and want to learn a causal model of the
wolffd@0 2619 data. We will not discuss the merits of this type of work here, but
wolffd@0 2620 merely show how to reproduce the results in the SGS book.
wolffd@0 2621 Their program,
wolffd@0 2622 <a href="http://hss.cmu.edu/html/departments/philosophy/TETRAD/tetrad.html">Tetrad</a>,
wolffd@0 2623 makes use of the Fisher Z-test for conditional
wolffd@0 2624 independence, so we do the same:
wolffd@0 2625 <pre>
wolffd@0 2626 max_fan_in = 4;
wolffd@0 2627 nsamples = 281;
wolffd@0 2628 alpha = 0.05;
wolffd@0 2629 pdag = learn_struct_pdag_pc('cond_indep_fisher_z', n, max_fan_in, C, nsamples, alpha);
wolffd@0 2630 </pre>
wolffd@0 2631 In this case, the CI test is
wolffd@0 2632 <pre>
wolffd@0 2633 f(X,Y,S) = cond_indep_fisher_z(X,Y,S, C,nsamples,alpha)
wolffd@0 2634 </pre>
wolffd@0 2635 The results match those of Fig 12a of SGS apart from two edge
wolffd@0 2636 differences; presumably this is due to rounding error (although it
wolffd@0 2637 could be a bug, either in BNT or in Tetrad).
wolffd@0 2638 This example can be found in the file BNT/examples/static/pc2.m.
wolffd@0 2639
wolffd@0 2640 <p>
wolffd@0 2641
wolffd@0 2642 The IC* algorithm (Pearl and Verma, 1991),
wolffd@0 2643 and the faster FCI algorithm (Spirtes, Glymour, and Scheines 1993),
wolffd@0 2644 are like the IC/PC algorithm, except that they can detect the presence
wolffd@0 2645 of latent variables.
wolffd@0 2646 See the file <tt>learn_struct_pdag_ic_star</tt> written by Tamar
wolffd@0 2647 Kushnir. The output is a matrix P, defined as follows
wolffd@0 2648 (see Pearl (2000), p52 for details):
wolffd@0 2649 <pre>
wolffd@0 2650 % P(i,j) = -1 if there is either a latent variable L such that i <-L->j OR there is a directed edge from i->j.
wolffd@0 2651 % P(i,j) = -2 if there is a marked directed i-*>j edge.
wolffd@0 2652 % P(i,j) = P(j,i) = 1 if there is and undirected edge i--j
wolffd@0 2653 % P(i,j) = P(j,i) = 2 if there is a latent variable L such that i<-L->j.
wolffd@0 2654 </pre>
wolffd@0 2655
wolffd@0 2656
wolffd@0 2657 <h2><a name="phl">Philippe Leray's structure learning package</h2>
wolffd@0 2658
wolffd@0 2659 Philippe Leray has written a
wolffd@0 2660 <a href="http://bnt.insa-rouen.fr/ajouts.html">
wolffd@0 2661 structure learning package</a> that uses BNT.
wolffd@0 2662
wolffd@0 2663 It currently (Juen 2003) has the following features:
wolffd@0 2664 <ul>
wolffd@0 2665 <li>PC with Chi2 statistical test
wolffd@0 2666 <li> MWST : Maximum weighted Spanning Tree
wolffd@0 2667 <li> Hill Climbing
wolffd@0 2668 <li> Greedy Search
wolffd@0 2669 <li> Structural EM
wolffd@0 2670 <li> hist_ic : optimal Histogram based on IC information criterion
wolffd@0 2671 <li> cpdag_to_dag
wolffd@0 2672 <li> dag_to_cpdag
wolffd@0 2673 <li> ...
wolffd@0 2674 </ul>
wolffd@0 2675
wolffd@0 2676
wolffd@0 2677 </a>
wolffd@0 2678
wolffd@0 2679
wolffd@0 2680 <!--
wolffd@0 2681 <h2><a name="read_learning">Further reading on learning</h2>
wolffd@0 2682
wolffd@0 2683 I recommend the following tutorials for more details on learning.
wolffd@0 2684 <ul>
wolffd@0 2685 <li> <a
wolffd@0 2686 href="http://www.cs.berkeley.edu/~murphyk/Papers/intel.ps.gz">My short
wolffd@0 2687 tutorial</a> on graphical models, which contains an overview of learning.
wolffd@0 2688
wolffd@0 2689 <li>
wolffd@0 2690 <A HREF="ftp://ftp.research.microsoft.com/pub/tr/TR-95-06.PS">
wolffd@0 2691 A tutorial on learning with Bayesian networks</a>, D. Heckerman,
wolffd@0 2692 Microsoft Research Tech Report, 1995.
wolffd@0 2693
wolffd@0 2694 <li> <A HREF="http://www-cad.eecs.berkeley.edu/~wray/Mirror/lwgmja">
wolffd@0 2695 Operations for Learning with Graphical Models</a>,
wolffd@0 2696 W. L. Buntine, JAIR'94, 159--225.
wolffd@0 2697 </ul>
wolffd@0 2698 <p>
wolffd@0 2699 -->
wolffd@0 2700
wolffd@0 2701
wolffd@0 2702
wolffd@0 2703
wolffd@0 2704
wolffd@0 2705 <h1><a name="engines">Inference engines</h1>
wolffd@0 2706
wolffd@0 2707 Up until now, we have used the junction tree algorithm for inference.
wolffd@0 2708 However, sometimes this is too slow, or not even applicable.
wolffd@0 2709 In general, there are many inference algorithms each of which make
wolffd@0 2710 different tradeoffs between speed, accuracy, complexity and
wolffd@0 2711 generality. Furthermore, there might be many implementations of the
wolffd@0 2712 same algorithm; for instance, a general purpose, readable version,
wolffd@0 2713 and a highly-optimized, specialized one.
wolffd@0 2714 To cope with this variety, we treat each inference algorithm as an
wolffd@0 2715 object, which we call an inference engine.
wolffd@0 2716
wolffd@0 2717 <p>
wolffd@0 2718 An inference engine is an object that contains a bnet and supports the
wolffd@0 2719 'enter_evidence' and 'marginal_nodes' methods. The engine constructor
wolffd@0 2720 takes the bnet as argument and may do some model-specific processing.
wolffd@0 2721 When 'enter_evidence' is called, the engine may do some
wolffd@0 2722 evidence-specific processing. Finally, when 'marginal_nodes' is
wolffd@0 2723 called, the engine may do some query-specific processing.
wolffd@0 2724
wolffd@0 2725 <p>
wolffd@0 2726 The amount of work done when each stage is specified -- structure,
wolffd@0 2727 parameters, evidence, and query -- depends on the engine. The cost of
wolffd@0 2728 work done early in this sequence can be amortized. On the other hand,
wolffd@0 2729 one can make better optimizations if one waits until later in the
wolffd@0 2730 sequence.
wolffd@0 2731 For example, the parameters might imply
wolffd@0 2732 conditional indpendencies that are not evident in the graph structure,
wolffd@0 2733 but can nevertheless be exploited; the evidence indicates which nodes
wolffd@0 2734 are observed and hence can effectively be disconnected from the
wolffd@0 2735 graph; and the query might indicate that large parts of the network
wolffd@0 2736 are d-separated from the query nodes. (Since it is not the actual
wolffd@0 2737 <em>values</em> of the evidence that matters, just which nodes are observed,
wolffd@0 2738 many engines allow you to specify which nodes will be observed when they are constructed,
wolffd@0 2739 i.e., before calling 'enter_evidence'. Some engines can still cope if
wolffd@0 2740 the actual pattern of evidence is different, e.g., if there is missing
wolffd@0 2741 data.)
wolffd@0 2742 <p>
wolffd@0 2743
wolffd@0 2744 Although being maximally lazy (i.e., only doing work when a query is
wolffd@0 2745 issued) may seem desirable,
wolffd@0 2746 this is not always the most efficient.
wolffd@0 2747 For example,
wolffd@0 2748 when learning using EM, we need to call marginal_nodes N times, where N is the
wolffd@0 2749 number of nodes. <a href="varelim">Variable elimination</a> would end
wolffd@0 2750 up repeating a lot of work
wolffd@0 2751 each time marginal_nodes is called, making it inefficient for
wolffd@0 2752 learning. The junction tree algorithm, by contrast, uses dynamic
wolffd@0 2753 programming to avoid this redundant computation --- it calculates all
wolffd@0 2754 marginals in two passes during 'enter_evidence', so calling
wolffd@0 2755 'marginal_nodes' takes constant time.
wolffd@0 2756 <p>
wolffd@0 2757 We will discuss some of the inference algorithms implemented in BNT
wolffd@0 2758 below, and finish with a <a href="#engine_summary">summary</a> of all
wolffd@0 2759 of them.
wolffd@0 2760
wolffd@0 2761
wolffd@0 2762
wolffd@0 2763
wolffd@0 2764
wolffd@0 2765
wolffd@0 2766
wolffd@0 2767 <h2><a name="varelim">Variable elimination</h2>
wolffd@0 2768
wolffd@0 2769 The variable elimination algorithm, also known as bucket elimination
wolffd@0 2770 or peeling, is one of the simplest inference algorithms.
wolffd@0 2771 The basic idea is to "push sums inside of products"; this is explained
wolffd@0 2772 in more detail
wolffd@0 2773 <a
wolffd@0 2774 href="http://HTTP.CS.Berkeley.EDU/~murphyk/Bayes/bayes.html#infer">here</a>.
wolffd@0 2775 <p>
wolffd@0 2776 The principle of distributing sums over products can be generalized
wolffd@0 2777 greatly to apply to any commutative semiring.
wolffd@0 2778 This forms the basis of many common algorithms, such as Viterbi
wolffd@0 2779 decoding and the Fast Fourier Transform. For details, see
wolffd@0 2780
wolffd@0 2781 <ul>
wolffd@0 2782 <li> R. McEliece and S. M. Aji, 2000.
wolffd@0 2783 <!--<a href="http://www.systems.caltech.edu/EE/Faculty/rjm/papers/GDL.ps">-->
wolffd@0 2784 <a href="GDL.pdf">
wolffd@0 2785 The Generalized Distributive Law</a>,
wolffd@0 2786 IEEE Trans. Inform. Theory, vol. 46, no. 2 (March 2000),
wolffd@0 2787 pp. 325--343.
wolffd@0 2788
wolffd@0 2789
wolffd@0 2790 <li>
wolffd@0 2791 F. R. Kschischang, B. J. Frey and H.-A. Loeliger, 2001.
wolffd@0 2792 <a href="http://www.cs.toronto.edu/~frey/papers/fgspa.abs.html">
wolffd@0 2793 Factor graphs and the sum-product algorithm</a>
wolffd@0 2794 IEEE Transactions on Information Theory, February, 2001.
wolffd@0 2795
wolffd@0 2796 </ul>
wolffd@0 2797
wolffd@0 2798 <p>
wolffd@0 2799 Choosing an order in which to sum out the variables so as to minimize
wolffd@0 2800 computational cost is known to be NP-hard.
wolffd@0 2801 The implementation of this algorithm in
wolffd@0 2802 <tt>var_elim_inf_engine</tt> makes no attempt to optimize this
wolffd@0 2803 ordering (in contrast, say, to <tt>jtree_inf_engine</tt>, which uses a
wolffd@0 2804 greedy search procedure to find a good ordering).
wolffd@0 2805 <p>
wolffd@0 2806 Note: unlike most algorithms, var_elim does all its computational work
wolffd@0 2807 inside of <tt>marginal_nodes</tt>, not inside of
wolffd@0 2808 <tt>enter_evidence</tt>.
wolffd@0 2809
wolffd@0 2810
wolffd@0 2811
wolffd@0 2812
wolffd@0 2813 <h2><a name="global">Global inference methods</h2>
wolffd@0 2814
wolffd@0 2815 The simplest inference algorithm of all is to explicitely construct
wolffd@0 2816 the joint distribution over all the nodes, and then to marginalize it.
wolffd@0 2817 This is implemented in <tt>global_joint_inf_engine</tt>.
wolffd@0 2818 Since the size of the joint is exponential in the
wolffd@0 2819 number of discrete (hidden) nodes, this is not a very practical algorithm.
wolffd@0 2820 It is included merely for pedagogical and debugging purposes.
wolffd@0 2821 <p>
wolffd@0 2822 Three specialized versions of this algorithm have also been implemented,
wolffd@0 2823 corresponding to the cases where all the nodes are discrete (D), all
wolffd@0 2824 are Gaussian (G), and some are discrete and some Gaussian (CG).
wolffd@0 2825 They are called <tt>enumerative_inf_engine</tt>,
wolffd@0 2826 <tt>gaussian_inf_engine</tt>,
wolffd@0 2827 and <tt>cond_gauss_inf_engine</tt> respectively.
wolffd@0 2828 <p>
wolffd@0 2829 Note: unlike most algorithms, these global inference algorithms do all their computational work
wolffd@0 2830 inside of <tt>marginal_nodes</tt>, not inside of
wolffd@0 2831 <tt>enter_evidence</tt>.
wolffd@0 2832
wolffd@0 2833
wolffd@0 2834 <h2><a name="quickscore">Quickscore</h2>
wolffd@0 2835
wolffd@0 2836 The junction tree algorithm is quite slow on the <a href="#qmr">QMR</a> network,
wolffd@0 2837 since the cliques are so big.
wolffd@0 2838 One simple trick we can use is to notice that hidden leaves do not
wolffd@0 2839 affect the posteriors on the roots, and hence do not need to be
wolffd@0 2840 included in the network.
wolffd@0 2841 A second trick is to notice that the negative findings can be
wolffd@0 2842 "absorbed" into the prior:
wolffd@0 2843 see the file
wolffd@0 2844 BNT/examples/static/mk_minimal_qmr_bnet for details.
wolffd@0 2845 <p>
wolffd@0 2846
wolffd@0 2847 A much more significant speedup is obtained by exploiting special
wolffd@0 2848 properties of the noisy-or node, as done by the quickscore
wolffd@0 2849 algorithm. For details, see
wolffd@0 2850 <ul>
wolffd@0 2851 <li> Heckerman, "A tractable inference algorithm for diagnosing multiple diseases", UAI 89.
wolffd@0 2852 <li> Rish and Dechter, "On the impact of causal independence", UCI
wolffd@0 2853 tech report, 1998.
wolffd@0 2854 </ul>
wolffd@0 2855
wolffd@0 2856 This has been implemented in BNT as a special-purpose inference
wolffd@0 2857 engine, which can be created and used as follows:
wolffd@0 2858 <pre>
wolffd@0 2859 engine = quickscore_inf_engine(inhibit, leak, prior);
wolffd@0 2860 engine = enter_evidence(engine, pos, neg);
wolffd@0 2861 m = marginal_nodes(engine, i);
wolffd@0 2862 </pre>
wolffd@0 2863
wolffd@0 2864
wolffd@0 2865 <h2><a name="belprop">Belief propagation</h2>
wolffd@0 2866
wolffd@0 2867 Even using quickscore, exact inference takes time that is exponential
wolffd@0 2868 in the number of positive findings.
wolffd@0 2869 Hence for large networks we need to resort to approximate inference techniques.
wolffd@0 2870 See for example
wolffd@0 2871 <ul>
wolffd@0 2872 <li> T. Jaakkola and M. Jordan, "Variational probabilistic inference and the
wolffd@0 2873 QMR-DT network", JAIR 10, 1999.
wolffd@0 2874
wolffd@0 2875 <li> K. Murphy, Y. Weiss and M. Jordan, "Loopy belief propagation for approximate inference: an empirical study",
wolffd@0 2876 UAI 99.
wolffd@0 2877 </ul>
wolffd@0 2878 The latter approximation
wolffd@0 2879 entails applying Pearl's belief propagation algorithm to a model even
wolffd@0 2880 if it has loops (hence the name loopy belief propagation).
wolffd@0 2881 Pearl's algorithm, implemented as <tt>pearl_inf_engine</tt>, gives
wolffd@0 2882 exact results when applied to singly-connected graphs
wolffd@0 2883 (a.k.a. polytrees, since
wolffd@0 2884 the underlying undirected topology is a tree, but a node may have
wolffd@0 2885 multiple parents).
wolffd@0 2886 To apply this algorithm to a graph with loops,
wolffd@0 2887 use <tt>pearl_inf_engine</tt>.
wolffd@0 2888 This can use a centralized or distributed message passing protocol.
wolffd@0 2889 You can use it as in the following example.
wolffd@0 2890 <pre>
wolffd@0 2891 engine = pearl_inf_engine(bnet, 'max_iter', 30);
wolffd@0 2892 engine = enter_evidence(engine, evidence);
wolffd@0 2893 m = marginal_nodes(engine, i);
wolffd@0 2894 </pre>
wolffd@0 2895 We found that this algorithm often converges, and when it does, often
wolffd@0 2896 is very accurate, but it depends on the precise setting of the
wolffd@0 2897 parameter values of the network.
wolffd@0 2898 (See the file BNT/examples/static/qmr1 to repeat the experiment for yourself.)
wolffd@0 2899 Understanding when and why belief propagation converges/ works
wolffd@0 2900 is a topic of ongoing research.
wolffd@0 2901 <p>
wolffd@0 2902 <tt>pearl_inf_engine</tt> can exploit special structure in noisy-or
wolffd@0 2903 and gmux nodes to compute messages efficiently.
wolffd@0 2904 <p>
wolffd@0 2905 <tt>belprop_inf_engine</tt> is like pearl, but uses potentials to
wolffd@0 2906 represent messages. Hence this is slower.
wolffd@0 2907 <p>
wolffd@0 2908 <tt>belprop_fg_inf_engine</tt> is like belprop,
wolffd@0 2909 but is designed for factor graphs.
wolffd@0 2910
wolffd@0 2911
wolffd@0 2912
wolffd@0 2913 <h2><a name="sampling">Sampling</h2>
wolffd@0 2914
wolffd@0 2915 BNT now (Mar '02) has two sampling (Monte Carlo) inference algorithms:
wolffd@0 2916 <ul>
wolffd@0 2917 <li> <tt>likelihood_weighting_inf_engine</tt> which does importance
wolffd@0 2918 sampling and can handle any node type.
wolffd@0 2919 <li> <tt>gibbs_sampling_inf_engine</tt>, written by Bhaskara Marthi.
wolffd@0 2920 Currently this can only handle tabular CPDs.
wolffd@0 2921 For a much faster and more powerful Gibbs sampling program, see
wolffd@0 2922 <a href="http://www.mrc-bsu.cam.ac.uk/bugs">BUGS</a>.
wolffd@0 2923 </ul>
wolffd@0 2924 Note: To generate samples from a network (which is not the same as inference!),
wolffd@0 2925 use <tt>sample_bnet</tt>.
wolffd@0 2926
wolffd@0 2927
wolffd@0 2928
wolffd@0 2929 <h2><a name="engine_summary">Summary of inference engines</h2>
wolffd@0 2930
wolffd@0 2931
wolffd@0 2932 The inference engines differ in many ways. Here are
wolffd@0 2933 some of the major "axes":
wolffd@0 2934 <ul>
wolffd@0 2935 <li> Works for all topologies or makes restrictions?
wolffd@0 2936 <li> Works for all node types or makes restrictions?
wolffd@0 2937 <li> Exact or approximate inference?
wolffd@0 2938 </ul>
wolffd@0 2939
wolffd@0 2940 <p>
wolffd@0 2941 In terms of topology, most engines handle any kind of DAG.
wolffd@0 2942 <tt>belprop_fg</tt> does approximate inference on factor graphs (FG), which
wolffd@0 2943 can be used to represent directed, undirected, and mixed (chain)
wolffd@0 2944 graphs.
wolffd@0 2945 (In the future, we plan to support exact inference on chain graphs.)
wolffd@0 2946 <tt>quickscore</tt> only works on QMR-like models.
wolffd@0 2947 <p>
wolffd@0 2948 In terms of node types: algorithms that use potentials can handle
wolffd@0 2949 discrete (D), Gaussian (G) or conditional Gaussian (CG) models.
wolffd@0 2950 Sampling algorithms can essentially handle any kind of node (distribution).
wolffd@0 2951 Other algorithms make more restrictive assumptions in exchange for
wolffd@0 2952 speed.
wolffd@0 2953 <p>
wolffd@0 2954 Finally, most algorithms are designed to give the exact answer.
wolffd@0 2955 The belief propagation algorithms are exact if applied to trees, and
wolffd@0 2956 in some other cases.
wolffd@0 2957 Sampling is considered approximate, even though, in the limit of an
wolffd@0 2958 infinite number of samples, it gives the exact answer.
wolffd@0 2959
wolffd@0 2960 <p>
wolffd@0 2961
wolffd@0 2962 Here is a summary of the properties
wolffd@0 2963 of all the engines in BNT which work on static networks.
wolffd@0 2964 <p>
wolffd@0 2965 <table>
wolffd@0 2966 <table border units = pixels><tr>
wolffd@0 2967 <td align=left width=0>Name
wolffd@0 2968 <td align=left width=0>Exact?
wolffd@0 2969 <td align=left width=0>Node type?
wolffd@0 2970 <td align=left width=0>topology
wolffd@0 2971 <tr>
wolffd@0 2972 <tr>
wolffd@0 2973 <td align=left> belprop
wolffd@0 2974 <td align=left> approx
wolffd@0 2975 <td align=left> D
wolffd@0 2976 <td align=left> DAG
wolffd@0 2977 <tr>
wolffd@0 2978 <td align=left> belprop_fg
wolffd@0 2979 <td align=left> approx
wolffd@0 2980 <td align=left> D
wolffd@0 2981 <td align=left> factor graph
wolffd@0 2982 <tr>
wolffd@0 2983 <td align=left> cond_gauss
wolffd@0 2984 <td align=left> exact
wolffd@0 2985 <td align=left> CG
wolffd@0 2986 <td align=left> DAG
wolffd@0 2987 <tr>
wolffd@0 2988 <td align=left> enumerative
wolffd@0 2989 <td align=left> exact
wolffd@0 2990 <td align=left> D
wolffd@0 2991 <td align=left> DAG
wolffd@0 2992 <tr>
wolffd@0 2993 <td align=left> gaussian
wolffd@0 2994 <td align=left> exact
wolffd@0 2995 <td align=left> G
wolffd@0 2996 <td align=left> DAG
wolffd@0 2997 <tr>
wolffd@0 2998 <td align=left> gibbs
wolffd@0 2999 <td align=left> approx
wolffd@0 3000 <td align=left> D
wolffd@0 3001 <td align=left> DAG
wolffd@0 3002 <tr>
wolffd@0 3003 <td align=left> global_joint
wolffd@0 3004 <td align=left> exact
wolffd@0 3005 <td align=left> D,G,CG
wolffd@0 3006 <td align=left> DAG
wolffd@0 3007 <tr>
wolffd@0 3008 <td align=left> jtree
wolffd@0 3009 <td align=left> exact
wolffd@0 3010 <td align=left> D,G,CG
wolffd@0 3011 <td align=left> DAG
wolffd@0 3012 b<tr>
wolffd@0 3013 <td align=left> likelihood_weighting
wolffd@0 3014 <td align=left> approx
wolffd@0 3015 <td align=left> any
wolffd@0 3016 <td align=left> DAG
wolffd@0 3017 <tr>
wolffd@0 3018 <td align=left> pearl
wolffd@0 3019 <td align=left> approx
wolffd@0 3020 <td align=left> D,G
wolffd@0 3021 <td align=left> DAG
wolffd@0 3022 <tr>
wolffd@0 3023 <td align=left> pearl
wolffd@0 3024 <td align=left> exact
wolffd@0 3025 <td align=left> D,G
wolffd@0 3026 <td align=left> polytree
wolffd@0 3027 <tr>
wolffd@0 3028 <td align=left> quickscore
wolffd@0 3029 <td align=left> exact
wolffd@0 3030 <td align=left> noisy-or
wolffd@0 3031 <td align=left> QMR
wolffd@0 3032 <tr>
wolffd@0 3033 <td align=left> stab_cond_gauss
wolffd@0 3034 <td align=left> exact
wolffd@0 3035 <td align=left> CG
wolffd@0 3036 <td align=left> DAG
wolffd@0 3037 <tr>
wolffd@0 3038 <td align=left> var_elim
wolffd@0 3039 <td align=left> exact
wolffd@0 3040 <td align=left> D,G,CG
wolffd@0 3041 <td align=left> DAG
wolffd@0 3042 </table>
wolffd@0 3043
wolffd@0 3044
wolffd@0 3045
wolffd@0 3046 <h1><a name="influence">Influence diagrams/ decision making</h1>
wolffd@0 3047
wolffd@0 3048 BNT implements an exact algorithm for solving LIMIDs (limited memory
wolffd@0 3049 influence diagrams), described in
wolffd@0 3050 <ul>
wolffd@0 3051 <li> S. L. Lauritzen and D. Nilsson.
wolffd@0 3052 <a href="http://www.math.auc.dk/~steffen/papers/limids.pdf">
wolffd@0 3053 Representing and solving decision problems with limited
wolffd@0 3054 information</a>
wolffd@0 3055 Management Science, 47, 1238 - 1251. September 2001.
wolffd@0 3056 </ul>
wolffd@0 3057 LIMIDs explicitely show all information arcs, rather than implicitely
wolffd@0 3058 assuming no forgetting. This allows them to model forgetful
wolffd@0 3059 controllers.
wolffd@0 3060 <p>
wolffd@0 3061 See the examples in <tt>BNT/examples/limids</tt> for details.
wolffd@0 3062
wolffd@0 3063
wolffd@0 3064
wolffd@0 3065
wolffd@0 3066 <h1>DBNs, HMMs, Kalman filters and all that</h1>
wolffd@0 3067
wolffd@0 3068 Click <a href="usage_dbn.html">here</a> for documentation about how to
wolffd@0 3069 use BNT for dynamical systems and sequence data.
wolffd@0 3070
wolffd@0 3071
wolffd@0 3072 </BODY>