Mercurial > hg > camir-aes2014

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