# HG changeset patch # User matthiasm # Date 1397228065 -3600 # Node ID b5b38998ef3b116a2806111cd93635e87fca6570 # Parent 12abff5474c8d05dee6b17f57b648ac66819e1be added all that other stuff diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/@assocarray/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/@assocarray/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,3 @@ +/assocarray.m/1.1.1.1/Wed May 29 15:59:52 2002// +/subsref.m/1.1.1.1/Wed Aug 4 19:36:30 2004// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/@assocarray/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/@assocarray/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/@assocarray diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/@assocarray/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/@assocarray/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/@assocarray/assocarray.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/@assocarray/assocarray.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,10 @@ +function A = assocarray(keys, vals) +% ASSOCARRAY Make an associative array +% function A = assocarray(keys, vals) +% +% keys{i} is the i'th string, vals{i} is the i'th value. +% After construction, A('foo') will return the value associated with foo. + +A.keys = keys; +A.vals = vals; +A = class(A, 'assocarray'); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/@assocarray/subsref.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/@assocarray/subsref.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,15 @@ +function val = subsref(A, S) +% SUBSREF Subscript reference for an associative array +% A('foo') will return the value associated with foo. +% If there are multiple identicaly keys, the first match is returned. +% Currently the search is sequential. + +i = 1; +while i <= length(A.keys) + if strcmp(S.subs{1}, A.keys{i}) + val = A.vals{i}; + return; + end + i = i + 1; +end +error(['can''t find ' S.subs{1}]) diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@boolean_CPD/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@boolean_CPD/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,2 @@ +/boolean_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@boolean_CPD/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@boolean_CPD/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@boolean_CPD diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@boolean_CPD/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@boolean_CPD/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@boolean_CPD/boolean_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@boolean_CPD/boolean_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,179 @@ +function CPD = boolean_CPD(bnet, self, ftype, fname, pfail) +% BOOLEAN_CPD Make a tabular CPD representing a (noisy) boolean function +% +% CPD = boolean_cpd(bnet, self, 'inline', f) uses the inline function f +% to specify the CPT. +% e.g., suppose X4 = X2 AND (NOT X3). Then we can write +% bnet.CPD{4} = boolean_CPD(bnet, 4, 'inline', inline('(x(1) & ~x(2)')); +% Note that x(1) refers pvals(1) = X2, and x(2) refers to pvals(2)=X3. +% +% CPD = boolean_cpd(bnet, self, 'named', f) assumes f is a function name. +% f can be built-in to matlab, or a file. +% e.g., If X4 = X2 AND X3, we can write +% bnet.CPD{4} = boolean_CPD(bnet, 4, 'named', 'and'); +% e.g., If X4 = X2 OR X3, we can write +% bnet.CPD{4} = boolean_CPD(bnet, 4, 'named', 'any'); +% +% CPD = boolean_cpd(bnet, self, 'rnd') makes a random non-redundant bool fn. +% +% CPD = boolean_CPD(bnet, self, 'inline'/'named', f, pfail) +% will put probability mass 1-pfail on f(parents), and put pfail on the other value. +% This is useful for simulating noisy boolean functions. +% If pfail is omitted, it is set to 0. +% (Note that adding noise to a random (non-redundant) boolean function just creates a different +% (potentially redundant) random boolean function.) +% +% Note: This cannot be used to simulate a noisy-OR gate. +% Example: suppose C has parents A and B, and the +% link of A->C fails with prob pA and the link B->C fails with pB. +% Then the noisy-OR gate defines the following distribution +% +% A B P(C=0) +% 0 0 1.0 +% 1 0 pA +% 0 1 pB +% 1 1 pA * PB +% +% By contrast, boolean_CPD(bnet, C, 'any', p) would define +% +% A B P(C=0) +% 0 0 1-p +% 1 0 p +% 0 1 p +% 1 1 p + + +if nargin==0 + % This occurs if we are trying to load an object from a file. + CPD = tabular_CPD(bnet, self); + return; +elseif isa(bnet, 'boolean_CPD') + % This might occur if we are copying an object. + CPD = bnet; + return; +end + +if nargin < 5, pfail = 0; end + +ps = parents(bnet.dag, self); +ns = bnet.node_sizes; +psizes = ns(ps); +self_size = ns(self); + +psucc = 1-pfail; + +k = length(ps); +switch ftype + case 'inline', f = eval_bool_fn(fname, k); + case 'named', f = eval_bool_fn(fname, k); + case 'rnd', f = mk_rnd_bool_fn(k); + otherwise, error(['unknown function type ' ftype]); +end + +CPT = zeros(prod(psizes), self_size); +ndx = find(f==0); +CPT(ndx, 1) = psucc; +CPT(ndx, 2) = pfail; +ndx = find(f==1); +CPT(ndx, 2) = psucc; +CPT(ndx, 1) = pfail; +if k > 0 + CPT = reshape(CPT, [psizes self_size]); +end + +clamp = 1; +CPD = tabular_CPD(bnet, self, CPT, [], clamp); + + + +%%%%%%%%%%%% + +function f = eval_bool_fn(fname, n) +% EVAL_BOOL_FN Evaluate a boolean function on all bit vectors of length n +% f = eval_bool_fn(fname, n) +% +% e.g. f = eval_bool_fn(inline('x(1) & x(3)'), 3) +% returns 0 0 0 0 0 1 0 1 + +ns = 2*ones(1, n); +f = zeros(1, 2^n); +bits = ind2subv(ns, 1:2^n); +for i=1:2^n + f(i) = feval(fname, bits(i,:)-1); +end + +%%%%%%%%%%%%%%% + +function f = mk_rnd_bool_fn(n) +% MK_RND_BOOL_FN Make a random bit vector of length n that encodes a non-redundant boolean function +% f = mk_rnd_bool_fn(n) + +red = 1; +while red + f = sample_discrete([0.5 0.5], 2^n, 1)-1; + red = redundant_bool_fn(f); +end + +%%%%%%%% + + +function red = redundant_bool_fn(f) +% REDUNDANT_BOOL_FN Does a boolean function depend on all its input values? +% r = redundant_bool_fn(f) +% +% f is a vector of length 2^n, representing the output for each bit vector. +% An input is redundant if there is no assignment to the other bits +% which changes the output e.g., input 1 is redundant if u(2:n) s.t., +% f([0 u(2:n)]) <> f([1 u(2:n)]). +% A function is redundant it it has any redundant inputs. + +n = log2(length(f)); +ns = 2*ones(1,n); +red = 0; +for i=1:n + ens = ns; + ens(i) = 1; + U = ind2subv(ens, 1:2^(n-1)); + U(:,i) = 1; + f1 = f(subv2ind(ns, U)); + U(:,i) = 2; + f2 = f(subv2ind(ns, U)); + if isequal(f1, f2) + red = 1; + return; + end +end + + +%%%%%%%%%% + +function [b, iter] = rnd_truth_table(N) +% RND_TRUTH_TABLE Construct the output of a random truth table s.t. each input is non-redundant +% b = rnd_truth_table(N) +% +% N is the number of inputs. +% b is a random bit string of length N, representing the output of the truth table. +% Non-redundant means that, for each input position k, +% there are at least two bit patterns, u and v, that differ only in the k'th position, +% s.t., f(u) ~= f(v), where f is the function represented by b. +% We use rejection sampling to ensure non-redundancy. +% +% Example: b = [0 0 0 1 0 0 0 1] is indep of 3rd input (AND of inputs 1 and 2) + +bits = ind2subv(2*ones(1,N), 1:2^N)-1; +redundant = 1; +iter = 0; +while redundant & (iter < 4) + iter = iter + 1; + b = sample_discrete([0.5 0.5], 1, 2^N)-1; + redundant = 0; + for i=1:N + on = find(bits(:,i)==1); + off = find(bits(:,i)==0); + if isequal(b(on), b(off)) + redundant = 1; + break; + end + end +end + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@deterministic_CPD/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@deterministic_CPD/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,2 @@ +/deterministic_CPD.m/1.1.1.1/Mon Oct 7 13:26:36 2002// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@deterministic_CPD/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@deterministic_CPD/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@deterministic_CPD diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@deterministic_CPD/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@deterministic_CPD/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@deterministic_CPD/deterministic_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@deterministic_CPD/deterministic_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,59 @@ +function CPD = deterministic_CPD(bnet, self, fname, pfail) +% DETERMINISTIC_CPD Make a tabular CPD representing a (noisy) deterministic function +% +% CPD = deterministic_CPD(bnet, self, fname) +% This calls feval(fname, pvals) for each possible vector of parent values. +% e.g., suppose there are 2 ternary parents, then pvals = +% [1 1], [2 1], [3 1], [1 2], [2 2], [3 2], [1 3], [2 3], [3 3] +% If v = feval(fname, pvals(i)), then +% CPD(x | parents=pvals(i)) = 1 if x==v, and = 0 if x<>v +% e.g., suppose X4 = X2 AND (NOT X3). Then +% bnet.CPD{4} = deterministic_CPD(bnet, 4, inline('((x(1)-1) & ~(x(2)-1)) + 1')); +% Note that x(1) refers pvals(1) = X2, and x(2) refers to pvals(2)=X3 +% See also boolean_CPD. +% +% CPD = deterministic_CPD(bnet, self, fname, pfail) +% will put probability mass 1-pfail on f(parents), and distribute pfail over the other values. +% This is useful for simulating noisy deterministic functions. +% If pfail is omitted, it is set to 0. +% + + +if nargin==0 + % This occurs if we are trying to load an object from a file. + CPD = tabular_CPD(bnet, self); + return; +elseif isa(bnet, 'deterministic_CPD') + % This might occur if we are copying an object. + CPD = bnet; + return; +end + +if nargin < 4, pfail = 0; end + +ps = parents(bnet.dag, self); +ns = bnet.node_sizes; +psizes = ns(ps); +self_size = ns(self); + +psucc = 1-pfail; + +CPT = zeros(prod(psizes), self_size); +pvals = zeros(1, length(ps)); +for i=1:prod(psizes) + pvals = ind2subv(psizes, i); + x = feval(fname, pvals); + %fprintf('%d ', [pvals x]); fprintf('\n'); + if psucc == 1 + CPT(i, x) = 1; + else + CPT(i, x) = psucc; + rest = mysetdiff(1:self_size, x); + CPT(i, rest) = pfail/length(rest); + end +end +CPT = reshape(CPT, [psizes self_size]); + +CPD = tabular_CPD(bnet, self, 'CPT',CPT, 'clamped',1); + + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/CPD_to_lambda_msg.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/CPD_to_lambda_msg.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,16 @@ +function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence) +% CPD_TO_LAMBDA_MSG Compute lambda message (discrete) +% lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence) +% Pearl p183 eq 4.52 + +switch msg_type + case 'd', + T = prod_CPT_and_pi_msgs(CPD, n, ps, msg, p); + mysize = length(msg{n}.lambda); + lambda = dpot(n, mysize, msg{n}.lambda); + T = multiply_by_pot(T, lambda); + lam_msg = pot_to_marginal(marginalize_pot(T, p)); + lam_msg = lam_msg.T; + case 'g', + error('discrete_CPD can''t create Gaussian msgs') +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/CPD_to_pi.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/CPD_to_pi.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,13 @@ +function pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence) +% COMPUTE_PI Compute pi vector (discrete) +% pi = compute_pi(CPD, msg_type, n, ps, msg, evidence) +% Pearl p183 eq 4.51 + +switch msg_type + case 'd', + T = prod_CPT_and_pi_msgs(CPD, n, ps, msg); + pi = pot_to_marginal(marginalize_pot(T, n)); + pi = pi.T(:); + case 'g', + error('can only convert discrete CPD to Gaussian pi if observed') +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/CPD_to_scgpot.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/CPD_to_scgpot.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,25 @@ +function pot = CPD_to_scgpot(CPD, domain, ns, cnodes, evidence) +% CPD_TO_SCGPOT Convert a CPD to a CG potential, incorporating any evidence (discrete) +% pot = CPD_to_scgpot(CPD, domain, ns, cnodes, evidence) +% +% domain is the domain of CPD. +% node_sizes(i) is the size of node i. +% cnodes +% evidence{i} is the evidence on the i'th node. + +%odom = domain(~isemptycell(evidence(domain))); + +%vals = cat(1, evidence{odom}); +%map = find_equiv_posns(odom, domain); +%index = mk_multi_index(length(domain), map, vals); +CPT = CPD_to_CPT(CPD); +%CPT = CPT(index{:}); +CPT = CPT(:); +%ns(odom) = 1; +potarray = cell(1, length(CPT)); +for i=1:length(CPT) + %p = CPT(i); + potarray{i} = scgcpot(0, 0, CPT(i)); + %scpot{i} = scpot(0, 0); +end +pot = scgpot(domain, [], [], ns, potarray); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,15 @@ +/CPD_to_lambda_msg.m/1.1.1.1/Wed May 29 15:59:52 2002// +/CPD_to_pi.m/1.1.1.1/Wed May 29 15:59:52 2002// +/CPD_to_scgpot.m/1.1.1.1/Wed May 29 15:59:52 2002// +/README/1.1.1.1/Wed May 29 15:59:52 2002// +/convert_CPD_to_table_hidden_ps.m/1.1.1.1/Wed May 29 15:59:52 2002// +/convert_obs_CPD_to_table.m/1.1.1.1/Wed May 29 15:59:52 2002// +/convert_to_pot.m/1.1.1.1/Fri Feb 20 22:00:38 2004// +/convert_to_sparse_table.c/1.1.1.1/Wed May 29 15:59:52 2002// +/convert_to_table.m/1.1.1.1/Wed May 29 15:59:52 2002// +/discrete_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002// +/dom_sizes.m/1.1.1.1/Wed May 29 15:59:52 2002// +/log_prob_node.m/1.1.1.1/Wed May 29 15:59:52 2002// +/prob_node.m/1.1.1.1/Wed May 29 15:59:52 2002// +/sample_node.m/1.1.1.1/Wed May 29 15:59:52 2002// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/CVS/Entries.Log --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/CVS/Entries.Log Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,2 @@ +A D/Old//// +A D/private//// diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@discrete_CPD diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/Old/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/Old/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,5 @@ +/convert_to_pot.m/1.1.1.1/Wed May 29 15:59:52 2002// +/convert_to_table.m/1.1.1.1/Wed May 29 15:59:52 2002// +/prob_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002// +/prob_node.m/1.1.1.1/Wed May 29 15:59:52 2002// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/Old/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/Old/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@discrete_CPD/Old diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/Old/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/Old/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/Old/convert_to_pot.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/Old/convert_to_pot.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,44 @@ +function pot = convert_to_pot(CPD, pot_type, domain, evidence) +% CONVERT_TO_POT Convert a tabular CPD to one or more potentials +% pots = convert_to_pot(CPD, pot_type, domain, evidence) +% +% pots{i} = CPD evaluated using evidence(domain(:,i)) +% If 'domains' is a single row vector, pots will be an object, not a cell array. + +ncases = size(domain,2); +assert(ncases==1); % not yet vectorized + +sz = dom_sizes(CPD); +ns = zeros(1, max(domain)); +ns(domain) = sz; + +local_ev = evidence(domain); +obs_bitv = ~isemptycell(local_ev); +odom = domain(obs_bitv); +T = convert_to_table(CPD, domain, local_ev, obs_bitv); + +switch pot_type + case 'u', + pot = upot(domain, sz, T, 0*myones(sz)); + case 'd', + ns(odom) = 1; + pot = dpot(domain, ns(domain), T); + case {'c','g'}, + % Since we want the output to be a Gaussian, the whole family must be observed. + % In other words, the potential is really just a constant. + p = T; + %p = prob_node(CPD, evidence(domain(end)), evidence(domain(1:end-1))); + ns(domain) = 0; + pot = cpot(domain, ns(domain), log(p)); + case 'cg', + T = T(:); + ns(odom) = 1; + can = cell(1, length(T)); + for i=1:length(T) + can{i} = cpot([], [], log(T(i))); + end + pot = cgpot(domain, [], ns, can); + otherwise, + error(['unrecognized pot type ' pot_type]) +end + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/Old/convert_to_table.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/Old/convert_to_table.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,23 @@ +function T = convert_to_table(CPD, domain, local_ev, obs_bitv) +% CONVERT_TO_TABLE Convert a discrete CPD to a table +% function T = convert_to_table(CPD, domain, local_ev, obs_bitv) +% +% We convert the CPD to a CPT, and then lookup the evidence on the discrete parents. +% The resulting table can easily be converted to a potential. + + +CPT = CPD_to_CPT(CPD); +obs_child_only = ~any(obs_bitv(1:end-1)) & obs_bitv(end); + +if obs_child_only + sz = size(CPT); + CPT = reshape(CPT, prod(sz(1:end-1)), sz(end)); + o = local_ev{end}; + T = CPT(:, o); +else + odom = domain(obs_bitv); + vals = cat(1, local_ev{find(obs_bitv)}); % undo cell array + map = find_equiv_posns(odom, domain); + index = mk_multi_index(length(domain), map, vals); + T = CPT(index{:}); +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/Old/prob_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/Old/prob_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,25 @@ +function p = prob_CPD(CPD, domain, ns, cnodes, evidence) +% PROB_CPD Compute prob of a node given evidence on the parents (discrete) +% p = prob_CPD(CPD, domain, ns, cnodes, evidence) +% +% domain is the domain of CPD. +% node_sizes(i) is the size of node i. +% cnodes = all the cts nodes +% evidence{i} is the evidence on the i'th node. + +ps = domain(1:end-1); +self = domain(end); +CPT = CPD_to_CPT(CPD); + +if isempty(ps) + T = CPT; +else + assert(~any(isemptycell(evidence(ps)))); + pvals = cat(1, evidence{ps}); + i = subv2ind(ns(ps), pvals(:)'); + T = reshape(CPT, [prod(ns(ps)) ns(self)]); + T = T(i,:); +end +p = T(evidence{self}); + + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/Old/prob_node.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/Old/prob_node.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,51 @@ +function [P, p] = prob_node(CPD, self_ev, pev) +% PROB_NODE Compute prod_m P(x(i,m)| x(pi_i,m), theta_i) for node i (discrete) +% [P, p] = prob_node(CPD, self_ev, pev) +% +% self_ev(m) is the evidence on this node in case m. +% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents). +% (These may also be cell arrays.) +% +% p(m) = P(x(i,m)| x(pi_i,m), theta_i) +% P = prod p(m) + +if iscell(self_ev), usecell = 1; else usecell = 0; end + +ncases = length(self_ev); +sz = dom_sizes(CPD); + +nparents = length(sz)-1; +if nparents == 0 + assert(isempty(pev)); +else + assert(isequal(size(pev), [nparents ncases])); +end + +n = length(sz); +dom = 1:n; +p = zeros(1, ncases); +if nparents == 0 + for m=1:ncases + if usecell + evidence = {self_ev{m}}; + else + evidence = num2cell(self_ev(m)); + end + T = convert_to_table(CPD, dom, evidence); + p(m) = T; + end +else + for m=1:ncases + if usecell + evidence = cell(1,n); + evidence(1:n-1) = pev(:,m); + evidence(n) = self_ev(m); + else + evidence = num2cell([pev(:,m)', self_ev(m)]); + end + T = convert_to_table(CPD, dom, evidence); + p(m) = T; + end +end +P = prod(p); + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/README --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/README Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,5 @@ +Any CPD on a discrete child with discrete parents +can be represented as a table (although this might be quite big). +discrete_CPD uses this tabular representation to implement various +functions. Subtypes are free to implement more efficient versions. + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/convert_CPD_to_table_hidden_ps.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/convert_CPD_to_table_hidden_ps.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,20 @@ +function T = convert_CPD_to_table_hidden_ps(CPD, child_obs) +% CONVERT_CPD_TO_TABLE_HIDDEN_PS Convert a discrete CPD to a table +% T = convert_CPD_to_table_hidden_ps(CPD, child_obs) +% +% This is like convert_to_table, except that we are guaranteed that +% none of the parents have evidence on them. +% child_obs may be an integer (1,2,...) or []. + +CPT = CPD_to_CPT(CPD); +if isempty(child_obs) + T = CPT(:); +else + sz = dom_sizes(CPD); + if length(sz)==1 % no parents + T = CPT(child_obs); + else + CPT = reshape(CPT, prod(sz(1:end-1)), sz(end)); + T = CPT(:, child_obs); + end +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/convert_obs_CPD_to_table.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/convert_obs_CPD_to_table.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,13 @@ +function T = convert_to_table(CPD, domain, evidence) +% CONVERT_TO_TABLE Convert a discrete CPD to a table +% T = convert_to_table(CPD, domain, evidence) +% +% We convert the CPD to a CPT, and then lookup the evidence on the discrete parents. +% The resulting table can easily be converted to a potential. + +CPT = CPD_to_CPT(CPD); +odom = domain(~isemptycell(evidence(domain))); +vals = cat(1, evidence{odom}); +map = find_equiv_posns(odom, domain); +index = mk_multi_index(length(domain), map, vals); +T = CPT(index{:}); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/convert_to_pot.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/convert_to_pot.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,62 @@ +function pot = convert_to_pot(CPD, pot_type, domain, evidence) +% CONVERT_TO_POT Convert a discrete CPD to a potential +% pot = convert_to_pot(CPD, pot_type, domain, evidence) +% +% pots = CPD evaluated using evidence(domain) + +ncases = size(domain,2); +assert(ncases==1); % not yet vectorized + +sz = dom_sizes(CPD); +ns = zeros(1, max(domain)); +ns(domain) = sz; + +CPT1 = CPD_to_CPT(CPD); +spar = issparse(CPT1); +odom = domain(~isemptycell(evidence(domain))); +if spar + T = convert_to_sparse_table(CPD, domain, evidence); +else + T = convert_to_table(CPD, domain, evidence); +end + +switch pot_type + case 'u', + pot = upot(domain, sz, T, 0*myones(sz)); + case 'd', + ns(odom) = 1; + pot = dpot(domain, ns(domain), T); + case {'c','g'}, + % Since we want the output to be a Gaussian, the whole family must be observed. + % In other words, the potential is really just a constant. + p = T; + %p = prob_node(CPD, evidence(domain(end)), evidence(domain(1:end-1))); + ns(domain) = 0; + pot = cpot(domain, ns(domain), log(p)); + + case 'cg', + T = T(:); + ns(odom) = 1; + can = cell(1, length(T)); + for i=1:length(T) + if T(i) == 0 + can{i} = cpot([], [], -Inf); % bug fix by Bob Welch 20/2/04 + else + can{i} = cpot([], [], log(T(i))); + end; + end + pot = cgpot(domain, [], ns, can); + + case 'scg' + T = T(:); + ns(odom) = 1; + pot_array = cell(1, length(T)); + for i=1:length(T) + pot_array{i} = scgcpot([], [], T(i)); + end + pot = scgpot(domain, [], [], ns, pot_array); + + otherwise, + error(['unrecognized pot type ' pot_type]) +end + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@discrete_CPD/convert_to_sparse_table.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@discrete_CPD/convert_to_sparse_table.c Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,154 @@ +/* convert_to_sparse_table.c convert a sparse discrete CPD with evidence into sparse table */ +/* convert_to_pot.m located in ../CPDs/discrete_CPD call it */ +/* 3 input */ +/* CPD prhs[0] with 1D sparse CPT */ +/* domain prhs[1] */ +/* evidence prhs[2] */ +/* 1 output */ +/* T plhs[0] sparse table */ + +#include +#include "mex.h" + +void ind_subv(int index, const int *cumprod, const int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i 1e-3) | isinf(P) + if isinf(P) % Y is observed + Sigma_lambda = zeros(self_size, self_size); % infinite precision => 0 variance + mu_lambda = msg{n}.lambda.mu; % observed_value; + else + Sigma_lambda = inv(P); + mu_lambda = Sigma_lambda * msg{n}.lambda.info_state; + end + C = inv(Sigma_lambda + BSigma); + lam_msg.precision = Bi' * C * Bi; + lam_msg.info_state = Bi' * C * (mu_lambda - Bmu); + else + % method that uses matrix inversion lemma to avoid inverting P + A = inv(P + inv(BSigma)); + C = P - P*A*P; + lam_msg.precision = Bi' * C * Bi; + D = eye(self_size) - P*A; + z = msg{n}.lambda.info_state; + lam_msg.info_state = Bi' * (D*z - D*P*Bmu); + end +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/CPD_to_pi.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/CPD_to_pi.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,22 @@ +function pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence) +% CPD_TO_PI Compute the pi vector (gaussian) +% function pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence) + +switch msg_type + case 'd', + error('gaussian_CPD can''t create discrete msgs') + case 'g', + [m, Q, W] = gaussian_CPD_params_given_dps(CPD, [ps n], evidence); + cps = ps(CPD.cps); + cpsizes = CPD.sizes(CPD.cps); + pi.mu = m; + pi.Sigma = Q; + for k=1:length(cps) % only get pi msgs from cts parents + %bk = block(k, cpsizes); + bk = CPD.cps_block_ndx{k}; + Bk = W(:, bk); + m = msg{n}.pi_from_parent{k}; + pi.Sigma = pi.Sigma + Bk * m.Sigma * Bk'; + pi.mu = pi.mu + Bk * m.mu; % m.mu = u(k) + end +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/CPD_to_scgpot.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/CPD_to_scgpot.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,58 @@ +function pot = CPD_to_scgpot(CPD, domain, ns, cnodes, evidence) +% CPD_TO_CGPOT Convert a Gaussian CPD to a CG potential, incorporating any evidence +% pot = CPD_to_cgpot(CPD, domain, ns, cnodes, evidence) + +self = CPD.self; +dnodes = mysetdiff(1:length(ns), cnodes); +odom = domain(~isemptycell(evidence(domain))); +cdom = myintersect(cnodes, domain); +cheaddom = myintersect(self, domain); +ctaildom = mysetdiff(cdom,cheaddom); +ddom = myintersect(dnodes, domain); +cobs = myintersect(cdom, odom); +dobs = myintersect(ddom, odom); +ens = ns; % effective node size +ens(cobs) = 0; +ens(dobs) = 1; + +% Extract the params compatible with the observations (if any) on the discrete parents (if any) +% parents are all but the last domain element +ps = domain(1:end-1); +dps = myintersect(ps, ddom); +dops = myintersect(dps, odom); + +map = find_equiv_posns(dops, dps); +dpvals = cat(1, evidence{dops}); +index = mk_multi_index(length(dps), map, dpvals); + +dpsize = prod(ens(dps)); +cpsize = size(CPD.weights(:,:,1), 2); % cts parents size +ss = size(CPD.mean, 1); % self size +% the reshape acts like a squeeze +m = reshape(CPD.mean(:, index{:}), [ss dpsize]); +C = reshape(CPD.cov(:, :, index{:}), [ss ss dpsize]); +W = reshape(CPD.weights(:, :, index{:}), [ss cpsize dpsize]); + + +% Convert each conditional Gaussian to a canonical potential +pot = cell(1, dpsize); +for i=1:dpsize + %pot{i} = linear_gaussian_to_scgcpot(m(:,i), C(:,:,i), W(:,:,i), cdom, ns, cnodes, evidence); + pot{i} = scgcpot(ss, cpsize, 1, m(:,i), W(:,:,i), C(:,:,i)); +end + +pot = scgpot(ddom, cheaddom, ctaildom, ens, pot); + + +function pot = linear_gaussian_to_scgcpot(mu, Sigma, W, domain, ns, cnodes, evidence) +% LINEAR_GAUSSIAN_TO_CPOT Convert a linear Gaussian CPD to a stable conditional potential element. +% pot = linear_gaussian_to_cpot(mu, Sigma, W, domain, ns, cnodes, evidence) + +p = 1; +A = mu; +B = W; +C = Sigma; +ns(odom) = 0; +%pot = scgcpot(, ns(domain), p, A, B, C); + + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,20 @@ +/CPD_to_lambda_msg.m/1.1.1.1/Wed May 29 15:59:52 2002// +/CPD_to_pi.m/1.1.1.1/Wed May 29 15:59:52 2002// +/CPD_to_scgpot.m/1.1.1.1/Wed May 29 15:59:52 2002// +/adjustable_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002// +/convert_CPD_to_table_hidden_ps.m/1.1.1.1/Wed May 29 15:59:52 2002// +/convert_to_pot.m/1.1.1.1/Sun Mar 9 23:03:16 2003// +/convert_to_table.m/1.1.1.1/Sun May 11 23:31:54 2003// +/display.m/1.1.1.1/Wed May 29 15:59:52 2002// +/gaussian_CPD.m/1.1.1.1/Wed Jun 15 21:13:06 2005// +/gaussian_CPD_params_given_dps.m/1.1.1.1/Sun May 11 23:13:40 2003// +/get_field.m/1.1.1.1/Wed May 29 15:59:52 2002// +/learn_params.m/1.1.1.1/Thu Jun 10 01:28:10 2004// +/log_prob_node.m/1.1.1.1/Tue Sep 10 17:44:00 2002// +/maximize_params.m/1.1.1.1/Tue May 20 14:10:06 2003// +/maximize_params_debug.m/1.1.1.1/Fri Jan 31 00:13:10 2003// +/reset_ess.m/1.1.1.1/Wed May 29 15:59:52 2002// +/sample_node.m/1.1.1.1/Wed May 29 15:59:52 2002// +/set_fields.m/1.1.1.1/Wed May 29 15:59:52 2002// +/update_ess.m/1.1.1.1/Tue Jul 22 22:55:46 2003// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/CVS/Entries.Log --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/CVS/Entries.Log Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,2 @@ +A D/Old//// +A D/private//// diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@gaussian_CPD diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/Old/CPD_to_lambda_msg.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/Old/CPD_to_lambda_msg.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,64 @@ +function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p) +% CPD_TO_LAMBDA_MSG Compute lambda message (gaussian) +% lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p) +% Pearl p183 eq 4.52 + +switch msg_type + case 'd', + error('gaussian_CPD can''t create discrete msgs') + case 'g', + self_size = CPD.sizes(end); + if all(msg{n}.lambda.precision == 0) % no info to send on + lam_msg.precision = zeros(self_size); + lam_msg.info_state = zeros(self_size, 1); + return; + end + cpsizes = CPD.sizes(CPD.cps); + dpval = 1; + Q = CPD.cov(:,:,dpval); + Sigmai = Q; + wmu = zeros(self_size, 1); + for k=1:length(ps) + pk = ps(k); + if pk ~= p + bk = block(k, cpsizes); + Bk = CPD.weights(:, bk, dpval); + m = msg{n}.pi_from_parent{k}; + Sigmai = Sigmai + Bk * m.Sigma * Bk'; + wmu = wmu + Bk * m.mu; % m.mu = u(k) + end + end + % Sigmai = Q + sum_{k \neq i} B_k Sigma_k B_k' + i = find_equiv_posns(p, ps); + bi = block(i, cpsizes); + Bi = CPD.weights(:,bi, dpval); + + if 0 + P = msg{n}.lambda.precision; + if isinf(P) % inv(P)=Sigma_lambda=0 + precision_temp = inv(Sigmai); + lam_msg.precision = Bi' * precision_temp * Bi; + lam_msg.info_state = precision_temp * (msg{n}.lambda.mu - wmu); + else + A = inv(P + inv(Sigmai)); + precision_temp = P + P*A*P; + lam_msg.precision = Bi' * precision_temp * Bi; + self_size = length(P); + C = eye(self_size) + P*A; + z = msg{n}.lambda.info_state; + lam_msg.info_state = C*z - C*P*wmu; + end + end + + if isinf(msg{n}.lambda.precision) + Sigma_lambda = zeros(self_size, self_size); % infinite precision => 0 variance + mu_lambda = msg{n}.lambda.mu; % observed_value; + else + Sigma_lambda = inv(msg{n}.lambda.precision); + mu_lambda = Sigma_lambda * msg{n}.lambda.info_state; + end + precision_temp = inv(Sigma_lambda + Sigmai); + lam_msg.precision = Bi' * precision_temp * Bi; + lam_msg.info_state = Bi' * precision_temp * (mu_lambda - wmu); +end + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/Old/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/Old/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,7 @@ +/CPD_to_lambda_msg.m/1.1.1.1/Wed May 29 15:59:52 2002// +/gaussian_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002// +/log_prob_node.m/1.1.1.1/Wed May 29 15:59:52 2002// +/maximize_params.m/1.1.1.1/Thu Jan 30 22:38:16 2003// +/update_ess.m/1.1.1.1/Wed May 29 15:59:52 2002// +/update_tied_ess.m/1.1.1.1/Wed May 29 15:59:52 2002// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/Old/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/Old/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@gaussian_CPD/Old diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/Old/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/Old/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/Old/gaussian_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/Old/gaussian_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,184 @@ +function CPD = gaussian_CPD(varargin) +% GAUSSIAN_CPD Make a conditional linear Gaussian distrib. +% +% To define this CPD precisely, call the continuous (cts) parents (if any) X, +% the discrete parents (if any) Q, and this node Y. Then the distribution on Y is: +% - no parents: Y ~ N(mu, Sigma) +% - cts parents : Y|X=x ~ N(mu + W x, Sigma) +% - discrete parents: Y|Q=i ~ N(mu(i), Sigma(i)) +% - cts and discrete parents: Y|X=x,Q=i ~ N(mu(i) + W(i) x, Sigma(i)) +% +% CPD = gaussian_CPD(bnet, node, ...) will create a CPD with random parameters, +% where node is the number of a node in this equivalence class. +% +% The list below gives optional arguments [default value in brackets]. +% (Let ns(i) be the size of node i, X = ns(X), Y = ns(Y) and Q = prod(ns(Q)).) +% +% mean - mu(:,i) is the mean given Q=i [ randn(Y,Q) ] +% cov - Sigma(:,:,i) is the covariance given Q=i [ repmat(eye(Y,Y), [1 1 Q]) ] +% weights - W(:,:,i) is the regression matrix given Q=i [ randn(Y,X,Q) ] +% cov_type - if 'diag', Sigma(:,:,i) is diagonal [ 'full' ] +% tied_cov - if 1, we constrain Sigma(:,:,i) to be the same for all i [0] +% clamp_mean - if 1, we do not adjust mu(:,i) during learning [0] +% clamp_cov - if 1, we do not adjust Sigma(:,:,i) during learning [0] +% clamp_weights - if 1, we do not adjust W(:,:,i) during learning [0] +% cov_prior_weight - weight given to I prior for estimating Sigma [0.01] +% +% e.g., CPD = gaussian_CPD(bnet, i, 'mean', [0; 0], 'clamp_mean', 'yes') +% +% For backwards compatibility with BNT2, you can also specify the parameters in the following order +% CPD = gaussian_CPD(bnet, self, mu, Sigma, W, cov_type, tied_cov, clamp_mean, clamp_cov, clamp_weight) +% +% Sometimes it is useful to create an "isolated" CPD, without needing to pass in a bnet. +% In this case, you must specify the discrete and cts parents (dps, cps) and the family sizes, followed +% by the optional arguments above: +% CPD = gaussian_CPD('self', i, 'dps', dps, 'cps', cps, 'sz', fam_size, ...) + + +if nargin==0 + % This occurs if we are trying to load an object from a file. + CPD = init_fields; + clamp = 0; + CPD = class(CPD, 'gaussian_CPD', generic_CPD(clamp)); + return; +elseif isa(varargin{1}, 'gaussian_CPD') + % This might occur if we are copying an object. + CPD = varargin{1}; + return; +end +CPD = init_fields; + +CPD = class(CPD, 'gaussian_CPD', generic_CPD(0)); + + +% parse mandatory arguments +if ~isstr(varargin{1}) % pass in bnet + bnet = varargin{1}; + self = varargin{2}; + args = varargin(3:end); + ns = bnet.node_sizes; + ps = parents(bnet.dag, self); + dps = myintersect(ps, bnet.dnodes); + cps = myintersect(ps, bnet.cnodes); + fam_sz = ns([ps self]); +else + disp('parsing new style') + for i=1:2:length(varargin) + switch varargin{i}, + case 'self', self = varargin{i+1}; + case 'dps', dps = varargin{i+1}; + case 'cps', cps = varargin{i+1}; + case 'sz', fam_sz = varargin{i+1}; + end + end + ps = myunion(dps, cps); + args = varargin; +end + +CPD.self = self; +CPD.sizes = fam_sz; + +% Figure out which (if any) of the parents are discrete, and which cts, and how big they are +% dps = discrete parents, cps = cts parents +CPD.cps = find_equiv_posns(cps, ps); % cts parent index +CPD.dps = find_equiv_posns(dps, ps); +ss = fam_sz(end); +psz = fam_sz(1:end-1); +dpsz = prod(psz(CPD.dps)); +cpsz = sum(psz(CPD.cps)); + +% set default params +CPD.mean = randn(ss, dpsz); +CPD.cov = 100*repmat(eye(ss), [1 1 dpsz]); +CPD.weights = randn(ss, cpsz, dpsz); +CPD.cov_type = 'full'; +CPD.tied_cov = 0; +CPD.clamped_mean = 0; +CPD.clamped_cov = 0; +CPD.clamped_weights = 0; +CPD.cov_prior_weight = 0.01; + +nargs = length(args); +if nargs > 0 + if ~isstr(args{1}) + % gaussian_CPD(bnet, self, mu, Sigma, W, cov_type, tied_cov, clamp_mean, clamp_cov, clamp_weights) + if nargs >= 1 & ~isempty(args{1}), CPD.mean = args{1}; end + if nargs >= 2 & ~isempty(args{2}), CPD.cov = args{2}; end + if nargs >= 3 & ~isempty(args{3}), CPD.weights = args{3}; end + if nargs >= 4 & ~isempty(args{4}), CPD.cov_type = args{4}; end + if nargs >= 5 & ~isempty(args{5}) & strcmp(args{5}, 'tied'), CPD.tied_cov = 1; end + if nargs >= 6 & ~isempty(args{6}), CPD.clamped_mean = 1; end + if nargs >= 7 & ~isempty(args{7}), CPD.clamped_cov = 1; end + if nargs >= 8 & ~isempty(args{8}), CPD.clamped_weights = 1; end + else + CPD = set_fields(CPD, args{:}); + end +end + +% Make sure the matrices have 1 dimension per discrete parent. +% Bug fix due to Xuejing Sun 3/6/01 +CPD.mean = myreshape(CPD.mean, [ss ns(dps)]); +CPD.cov = myreshape(CPD.cov, [ss ss ns(dps)]); +CPD.weights = myreshape(CPD.weights, [ss cpsz ns(dps)]); + +CPD.init_cov = CPD.cov; % we reset to this if things go wrong during learning + +% expected sufficient statistics +CPD.Wsum = zeros(dpsz,1); +CPD.WYsum = zeros(ss, dpsz); +CPD.WXsum = zeros(cpsz, dpsz); +CPD.WYYsum = zeros(ss, ss, dpsz); +CPD.WXXsum = zeros(cpsz, cpsz, dpsz); +CPD.WXYsum = zeros(cpsz, ss, dpsz); + +% For BIC +CPD.nsamples = 0; +switch CPD.cov_type + case 'full', + ncov_params = ss*(ss-1)/2; % since symmetric (and positive definite) + case 'diag', + ncov_params = ss; + otherwise + error(['unrecognized cov_type ' cov_type]); +end +% params = weights + mean + cov +if CPD.tied_cov + CPD.nparams = ss*cpsz*dpsz + ss*dpsz + ncov_params; +else + CPD.nparams = ss*cpsz*dpsz + ss*dpsz + dpsz*ncov_params; +end + + + +clamped = CPD.clamped_mean & CPD.clamped_cov & CPD.clamped_weights; +CPD = set_clamped(CPD, clamped); + +%%%%%%%%%%% + +function CPD = init_fields() +% This ensures we define the fields in the same order +% no matter whether we load an object from a file, +% or create it from scratch. (Matlab requires this.) + +CPD.self = []; +CPD.sizes = []; +CPD.cps = []; +CPD.dps = []; +CPD.mean = []; +CPD.cov = []; +CPD.weights = []; +CPD.clamped_mean = []; +CPD.clamped_cov = []; +CPD.clamped_weights = []; +CPD.init_cov = []; +CPD.cov_type = []; +CPD.tied_cov = []; +CPD.Wsum = []; +CPD.WYsum = []; +CPD.WXsum = []; +CPD.WYYsum = []; +CPD.WXXsum = []; +CPD.WXYsum = []; +CPD.nsamples = []; +CPD.nparams = []; +CPD.cov_prior_weight = []; diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/Old/log_prob_node.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/Old/log_prob_node.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,59 @@ +function L = log_prob_node(CPD, self_ev, pev) +% LOG_PROB_NODE Compute prod_m log P(x(i,m)| x(pi_i,m), theta_i) for node i (gaussian) +% L = log_prob_node(CPD, self_ev, pev) +% +% self_ev(m) is the evidence on this node in case m. +% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents). +% (These may also be cell arrays.) + +if iscell(self_ev), usecell = 1; else usecell = 0; end + +use_log = 1; +ncases = length(self_ev); +nparents = length(CPD.sizes)-1; +assert(ncases == size(pev, 2)); + +if ncases == 0 + L = 0; + return; +end + +if length(CPD.dps)==0 % no discrete parents, so we can vectorize + i = 1; + if usecell + Y = cell2num(self_ev); + else + Y = self_ev; + end + if length(CPD.cps) == 0 + L = gaussian_prob(Y, CPD.mean(:,i), CPD.cov(:,:,i), use_log); + else + if usecell + X = cell2num(pev); + else + X = pev; + end + L = gaussian_prob(Y, CPD.mean(:,i) + CPD.weights(:,:,i)*X, CPD.cov(:,:,i), use_log); + end +else % each case uses a (potentially) different set of parameters + L = 0; + for m=1:ncases + if usecell + dpvals = cat(1, pev{CPD.dps, m}); + else + dpvals = pev(CPD.dps, m); + end + i = subv2ind(CPD.sizes(CPD.dps), dpvals(:)'); + y = self_ev{m}; + if length(CPD.cps) == 0 + L = L + gaussian_prob(y, CPD.mean(:,i), CPD.cov(:,:,i), use_log); + else + if usecell + x = cat(1, pev{CPD.cps, m}); + else + x = pev(CPD.cps, m); + end + L = L + gaussian_prob(y, CPD.mean(:,i) + CPD.weights(:,:,i)*x, CPD.cov(:,:,i), use_log); + end + end +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/Old/maximize_params.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/Old/maximize_params.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,147 @@ +function CPD = maximize_params(CPD, temp) +% MAXIMIZE_PARAMS Set the params of a CPD to their ML values (Gaussian) +% CPD = maximize_params(CPD, temperature) +% +% Temperature is currently only used for entropic prior on Sigma + +% For details, see "Fitting a Conditional Gaussian Distribution", Kevin Murphy, tech. report, +% 1998, available at www.cs.berkeley.edu/~murphyk/papers.html +% Refering to table 2, we use equations 1/2 to estimate the covariance matrix in the untied/tied case, +% and equation 9 to estimate the weight matrix and mean. +% We do not implement spherical Gaussians - the code is already pretty complicated! + +if ~adjustable_CPD(CPD), return; end + +%assert(approxeq(CPD.nsamples, sum(CPD.Wsum))); +assert(~any(isnan(CPD.WXXsum))) +assert(~any(isnan(CPD.WXYsum))) +assert(~any(isnan(CPD.WYYsum))) + +[self_size cpsize dpsize] = size(CPD.weights); + +% Append 1s to the parents, and derive the corresponding cross products. +% This is used when estimate the means and weights simultaneosuly, +% and when estimatting Sigma. +% Let x2 = [x 1]' +XY = zeros(cpsize+1, self_size, dpsize); % XY(:,:,i) = sum_l w(l,i) x2(l) y(l)' +XX = zeros(cpsize+1, cpsize+1, dpsize); % XX(:,:,i) = sum_l w(l,i) x2(l) x2(l)' +YY = zeros(self_size, self_size, dpsize); % YY(:,:,i) = sum_l w(l,i) y(l) y(l)' +for i=1:dpsize + XY(:,:,i) = [CPD.WXYsum(:,:,i) % X*Y + CPD.WYsum(:,i)']; % 1*Y + % [x * [x' 1] = [xx' x + % 1] x' 1] + XX(:,:,i) = [CPD.WXXsum(:,:,i) CPD.WXsum(:,i); + CPD.WXsum(:,i)' CPD.Wsum(i)]; + YY(:,:,i) = CPD.WYYsum(:,:,i); +end + +w = CPD.Wsum(:); +% Set any zeros to one before dividing +% This is valid because w(i)=0 => WYsum(:,i)=0, etc +w = w + (w==0); + +if CPD.clamped_mean + % Estimating B2 and then setting the last column (the mean) to the clamped mean is *not* equivalent + % to estimating B and then adding the clamped_mean to the last column. + if ~CPD.clamped_weights + B = zeros(self_size, cpsize, dpsize); + for i=1:dpsize + if det(CPD.WXXsum(:,:,i))==0 + B(:,:,i) = 0; + else + % Eqn 9 in table 2 of TR + %B(:,:,i) = CPD.WXYsum(:,:,i)' * inv(CPD.WXXsum(:,:,i)); + B(:,:,i) = (CPD.WXXsum(:,:,i) \ CPD.WXYsum(:,:,i))'; + end + end + %CPD.weights = reshape(B, [self_size cpsize dpsize]); + CPD.weights = B; + end +elseif CPD.clamped_weights % KPM 1/25/02 + if ~CPD.clamped_mean % ML estimate is just sample mean of the residuals + for i=1:dpsize + CPD.mean(:,i) = (CPD.WYsum(:,i) - CPD.weights(:,:,i) * CPD.WXsum(:,i)) / w(i); + end + end +else % nothing is clamped, so estimate mean and weights simultaneously + B2 = zeros(self_size, cpsize+1, dpsize); + for i=1:dpsize + if det(XX(:,:,i))==0 % fix by U. Sondhauss 6/27/99 + B2(:,:,i)=0; + else + % Eqn 9 in table 2 of TR + %B2(:,:,i) = XY(:,:,i)' * inv(XX(:,:,i)); + B2(:,:,i) = (XX(:,:,i) \ XY(:,:,i))'; + end + CPD.mean(:,i) = B2(:,cpsize+1,i); + CPD.weights(:,:,i) = B2(:,1:cpsize,i); + end +end + +% Let B2 = [W mu] +if cpsize>0 + B2(:,1:cpsize,:) = reshape(CPD.weights, [self_size cpsize dpsize]); +end +B2(:,cpsize+1,:) = reshape(CPD.mean, [self_size dpsize]); + +% To avoid singular covariance matrices, +% we use the regularization method suggested in "A Quasi-Bayesian approach to estimating +% parameters for mixtures of normal distributions", Hamilton 91. +% If the ML estimate is Sigma = M/N, the MAP estimate is (M+gamma*I) / (N+gamma), +% where gamma >=0 is a smoothing parameter (equivalent sample size of I prior) + +gamma = CPD.cov_prior_weight; + +if ~CPD.clamped_cov + if CPD.cov_prior_entropic % eqn 12 of Brand AI/Stat 99 + Z = 1-temp; + % When temp > 1, Z is negative, so we are dividing by a smaller + % number, ie. increasing the variance. + else + Z = 0; + end + if CPD.tied_cov + S = zeros(self_size, self_size); + % Eqn 2 from table 2 in TR + for i=1:dpsize + S = S + (YY(:,:,i) - B2(:,:,i)*XY(:,:,i)); + end + %denom = max(1, CPD.nsamples + gamma + Z); + denom = CPD.nsamples + gamma + Z; + S = (S + gamma*eye(self_size)) / denom; + if strcmp(CPD.cov_type, 'diag') + S = diag(diag(S)); + end + CPD.cov = repmat(S, [1 1 dpsize]); + else + for i=1:dpsize + % Eqn 1 from table 2 in TR + S = YY(:,:,i) - B2(:,:,i)*XY(:,:,i); + %denom = max(1, w(i) + gamma + Z); % gives wrong answers on mhmm1 + denom = w(i) + gamma + Z; + S = (S + gamma*eye(self_size)) / denom; + CPD.cov(:,:,i) = S; + end + if strcmp(CPD.cov_type, 'diag') + for i=1:dpsize + CPD.cov(:,:,i) = diag(diag(CPD.cov(:,:,i))); + end + end + end +end + + +check_covars = 0; +min_covar = 1e-5; +if check_covars % prevent collapsing to a point + for i=1:dpsize + if min(svd(CPD.cov(:,:,i))) < min_covar + disp(['resetting singular covariance for node ' num2str(CPD.self)]); + CPD.cov(:,:,i) = CPD.init_cov(:,:,i); + end + end +end + + + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/Old/update_ess.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/Old/update_ess.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,85 @@ +function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv) +% UPDATE_ESS Update the Expected Sufficient Statistics of a Gaussian node +% function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv) + +%if nargin < 6 +% hidden_bitv = zeros(1, max(fmarginal.domain)); +% hidden_bitv(find(isempty(evidence)))=1; +%end + +dom = fmarginal.domain; +self = dom(end); +ps = dom(1:end-1); +hidden_self = hidden_bitv(self); +cps = myintersect(ps, cnodes); +dps = mysetdiff(ps, cps); +hidden_cps = all(hidden_bitv(cps)); +hidden_dps = all(hidden_bitv(dps)); + +CPD.nsamples = CPD.nsamples + 1; +[ss cpsz dpsz] = size(CPD.weights); % ss = self size + +% Let X be the cts parent (if any), Y be the cts child (self). + +if ~hidden_self & (isempty(cps) | ~hidden_cps) & hidden_dps % all cts nodes are observed, all discrete nodes are hidden + % Since X and Y are observed, SYY = 0, SXX = 0, SXY = 0 + % Since discrete parents are hidden, we do not need to add evidence to w. + w = fmarginal.T(:); + CPD.Wsum = CPD.Wsum + w; + y = evidence{self}; + Cyy = y*y'; + if ~CPD.useC + W = repmat(w(:)',ss,1); % W(y,i) = w(i) + W2 = repmat(reshape(W, [ss 1 dpsz]), [1 ss 1]); % W2(x,y,i) = w(i) + CPD.WYsum = CPD.WYsum + W .* repmat(y(:), 1, dpsz); + CPD.WYYsum = CPD.WYYsum + W2 .* repmat(reshape(Cyy, [ss ss 1]), [1 1 dpsz]); + else + W = w(:)'; + W2 = reshape(W, [1 1 dpsz]); + CPD.WYsum = CPD.WYsum + rep_mult(W, y(:), size(CPD.WYsum)); + CPD.WYYsum = CPD.WYYsum + rep_mult(W2, Cyy, size(CPD.WYYsum)); + end + if cpsz > 0 % X exists + x = cat(1, evidence{cps}); x = x(:); + Cxx = x*x'; + Cxy = x*y'; + if ~CPD.useC + CPD.WXsum = CPD.WXsum + W .* repmat(x(:), 1, dpsz); + CPD.WXXsum = CPD.WXXsum + W2 .* repmat(reshape(Cxx, [cpsz cpsz 1]), [1 1 dpsz]); + CPD.WXYsum = CPD.WXYsum + W2 .* repmat(reshape(Cxy, [cpsz ss 1]), [1 1 dpsz]); + else + CPD.WXsum = CPD.WXsum + rep_mult(W, x(:), size(CPD.WXsum)); + CPD.WXXsum = CPD.WXXsum + rep_mult(W2, Cxx, size(CPD.WXXsum)); + CPD.WXYsum = CPD.WXYsum + rep_mult(W2, Cxy, size(CPD.WXYsum)); + end + end + return; +end + +% general (non-vectorized) case +fullm = add_evidence_to_gmarginal(fmarginal, evidence, ns, cnodes); % slow! + +if dpsz == 1 % no discrete parents + w = 1; +else + w = fullm.T(:); +end + +CPD.Wsum = CPD.Wsum + w; +xi = 1:cpsz; +yi = (cpsz+1):(cpsz+ss); +for i=1:dpsz + muY = fullm.mu(yi, i); + SYY = fullm.Sigma(yi, yi, i); + CPD.WYsum(:,i) = CPD.WYsum(:,i) + w(i)*muY; + CPD.WYYsum(:,:,i) = CPD.WYYsum(:,:,i) + w(i)*(SYY + muY*muY'); % E[X Y] = Cov[X,Y] + E[X] E[Y] + if cpsz > 0 + muX = fullm.mu(xi, i); + SXX = fullm.Sigma(xi, xi, i); + SXY = fullm.Sigma(xi, yi, i); + CPD.WXsum(:,i) = CPD.WXsum(:,i) + w(i)*muX; + CPD.WXXsum(:,:,i) = CPD.WXXsum(:,:,i) + w(i)*(SXX + muX*muX'); + CPD.WXYsum(:,:,i) = CPD.WXYsum(:,:,i) + w(i)*(SXY + muX*muY'); + end +end + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/Old/update_tied_ess.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/Old/update_tied_ess.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,118 @@ +function CPD = update_tied_ess(CPD, domain, engine, evidence, ns, cnodes) + +if ~adjustable_CPD(CPD), return; end +nCPDs = size(domain, 2); +fmarginal = cell(1, nCPDs); +for l=1:nCPDs + fmarginal{l} = marginal_family(engine, nodes(l)); +end + +[ss cpsz dpsz] = size(CPD.weights); +if const_evidence_pattern(engine) + dom = domain(:,1); + dnodes = mysetdiff(1:length(ns), cnodes); + ddom = myintersect(dom, dnodes); + cdom = myintersect(dom, cnodes); + odom = dom(~isemptycell(evidence(dom))); + hdom = dom(isemptycell(evidence(dom))); + % If all hidden nodes are discrete and all cts nodes are observed + % (e.g., HMM with Gaussian output) + % we can add the observed evidence in parallel + if mysubset(ddom, hdom) & mysubset(cdom, odom) + [mu, Sigma, T] = add_cts_ev_to_marginals(fmarginal, evidence, ns, cnodes); + else + mu = zeros(ss, dpsz, nCPDs); + Sigma = zeros(ss, ss, dpsz, nCPDs); + T = zeros(dpsz, nCPDs); + for l=1:nCPDs + [mu(:,:,l), Sigma(:,:,:,l), T(:,l)] = add_ev_to_marginals(fmarginal{l}, evidence, ns, cnodes); + end + end +end +CPD.nsamples = CPD.nsamples + nCPDs; + + +if dpsz == 1 % no discrete parents + w = 1; +else + w = fullm.T(:); +end +CPD.Wsum = CPD.Wsum + w; +% Let X be the cts parent (if any), Y be the cts child (self). +xi = 1:cpsz; +yi = (cpsz+1):(cpsz+ss); +for i=1:dpsz + muY = fullm.mu(yi, i); + SYY = fullm.Sigma(yi, yi, i); + CPD.WYsum(:,i) = CPD.WYsum(:,i) + w(i)*muY; + CPD.WYYsum(:,:,i) = CPD.WYYsum(:,:,i) + w(i)*(SYY + muY*muY'); % E[X Y] = Cov[X,Y] + E[X] E[Y] + if cpsz > 0 + muX = fullm.mu(xi, i); + SXX = fullm.Sigma(xi, xi, i); + SXY = fullm.Sigma(xi, yi, i); + CPD.WXsum(:,i) = CPD.WXsum(:,i) + w(i)*muX; + CPD.WXYsum(:,:,i) = CPD.WXYsum(:,:,i) + w(i)*(SXY + muX*muY'); + CPD.WXXsum(:,:,i) = CPD.WXXsum(:,:,i) + w(i)*(SXX + muX*muX'); + end +end + + +%%%%%%%%%%%%% + +function fullm = add_evidence_to_marginal(fmarginal, evidence, ns, cnodes) + + +dom = fmarginal.domain; + +% Find out which values of the discrete parents (if any) are compatible with +% the discrete evidence (if any). +dnodes = mysetdiff(1:length(ns), cnodes); +ddom = myintersect(dom, dnodes); +cdom = myintersect(dom, cnodes); +odom = dom(~isemptycell(evidence(dom))); +hdom = dom(isemptycell(evidence(dom))); + +dobs = myintersect(ddom, odom); +dvals = cat(1, evidence{dobs}); +ens = ns; % effective node sizes +ens(dobs) = 1; +S = prod(ens(ddom)); +subs = ind2subv(ens(ddom), 1:S); +mask = find_equiv_posns(dobs, ddom); +subs(mask) = dvals; +supportedQs = subv2ind(ns(ddom), subs); + +if isempty(ddom) + Qarity = 1; +else + Qarity = prod(ns(ddom)); +end +fullm.T = zeros(Qarity, 1); +fullm.T(supportedQs) = fmarginal.T(:); + +% Now put the hidden cts parts into their right blocks, +% leaving the observed cts parts as 0. +cobs = myintersect(cdom, odom); +chid = myintersect(cdom, hdom); +cvals = cat(1, evidence{cobs}); +n = sum(ns(cdom)); +fullm.mu = zeros(n,Qarity); +fullm.Sigma = zeros(n,n,Qarity); + +if ~isempty(chid) + chid_blocks = block(find_equiv_posns(chid, cdom), ns(cdom)); +end +if ~isempty(cobs) + cobs_blocks = block(find_equiv_posns(cobs, cdom), ns(cdom)); +end + +for i=1:length(supportedQs) + Q = supportedQs(i); + if ~isempty(chid) + fullm.mu(chid_blocks, Q) = fmarginal.mu(:, i); + fullm.Sigma(chid_blocks, chid_blocks, Q) = fmarginal.Sigma(:,:,i); + end + if ~isempty(cobs) + fullm.mu(cobs_blocks, Q) = cvals(:); + end +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/adjustable_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/adjustable_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,5 @@ +function p = adjustable_CPD(CPD) +% ADJUSTABLE_CPD Does this CPD have any adjustable params? (gaussian) +% p = adjustable_CPD(CPD) + +p = ~CPD.clamped_mean | ~CPD.clamped_cov | ~CPD.clamped_weights; diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/convert_CPD_to_table_hidden_ps.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/convert_CPD_to_table_hidden_ps.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,20 @@ +function T = convert_CPD_to_table_hidden_ps(CPD, self_val) +% CONVERT_CPD_TO_TABLE_HIDDEN_PS Convert a Gaussian CPD to a table +% function T = convert_CPD_to_table_hidden_ps(CPD, self_val) +% +% self_val must be a non-empty vector. +% All the parents are hidden. +% +% This is used by misc/convert_dbn_CPDs_to_tables + +m = CPD.mean; +C = CPD.cov; +W = CPD.weights; + +[ssz dpsize] = size(m); + +T = zeros(dpsize, 1); +for i=1:dpsize + T(i) = gaussian_prob(self_val, m(:,i), C(:,:,i)); +end + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/convert_to_pot.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/convert_to_pot.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,71 @@ +function pot = convert_to_pot(CPD, pot_type, domain, evidence) +% CONVERT_TO_POT Convert a Gaussian CPD to one or more potentials +% pot = convert_to_pot(CPD, pot_type, domain, evidence) + +sz = CPD.sizes; +ns = zeros(1, max(domain)); +ns(domain) = sz; + +odom = domain(~isemptycell(evidence(domain))); +ps = domain(1:end-1); +cps = ps(CPD.cps); +dps = ps(CPD.dps); +self = domain(end); +cdom = [cps(:)' self]; +ddom = dps; +cnodes = cdom; + +switch pot_type + case 'u', + error('gaussian utility potentials not yet supported'); + + case 'd', + T = convert_to_table(CPD, domain, evidence); + ns(odom) = 1; + pot = dpot(domain, ns(domain), T); + + case {'c','g'}, + [m, C, W] = gaussian_CPD_params_given_dps(CPD, domain, evidence); + pot = linear_gaussian_to_cpot(m, C, W, domain, ns, cnodes, evidence); + + case 'cg', + [m, C, W] = gaussian_CPD_params_given_dps(CPD, domain, evidence); + % Convert each conditional Gaussian to a canonical potential + cobs = myintersect(cdom, odom); + dobs = myintersect(ddom, odom); + ens = ns; % effective node size + ens(cobs) = 0; + ens(dobs) = 1; + dpsize = prod(ens(dps)); + can = cell(1, dpsize); + for i=1:dpsize + if isempty(W) + can{i} = linear_gaussian_to_cpot(m(:,i), C(:,:,i), [], cdom, ns, cnodes, evidence); + else + can{i} = linear_gaussian_to_cpot(m(:,i), C(:,:,i), W(:,:,i), cdom, ns, cnodes, evidence); + end + end + pot = cgpot(ddom, cdom, ens, can); + + case 'scg', + [m, C, W] = gaussian_CPD_params_given_dps(CPD, domain, evidence); + cobs = myintersect(cdom, odom); + dobs = myintersect(ddom, odom); + ens = ns; % effective node size + ens(cobs) = 0; + ens(dobs) = 1; + dpsize = prod(ens(dps)); + cpsize = size(W, 2); % cts parents size + ss = size(m, 1); % self size + cheaddom = self; + ctaildom = cps(:)'; + pot_array = cell(1, dpsize); + for i=1:dpsize + pot_array{i} = scgcpot(ss, cpsize, 1, m(:,i), W(:,:,i), C(:,:,i)); + end + pot = scgpot(ddom, cheaddom, ctaildom, ens, pot_array); + + otherwise, + error(['unrecognized pot_type' pot_type]) +end + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/convert_to_table.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/convert_to_table.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,38 @@ +function T = convert_to_table(CPD, domain, evidence) +% CONVERT_TO_TABLE Convert a Gaussian CPD to a table +% T = convert_to_table(CPD, domain, evidence) + + +sz = CPD.sizes; +ns = zeros(1, max(domain)); +ns(domain) = sz; + +odom = domain(~isemptycell(evidence(domain))); +ps = domain(1:end-1); +cps = ps(CPD.cps); +dps = ps(CPD.dps); +self = domain(end); +cdom = [cps(:)' self]; +ddom = dps; +cnodes = cdom; + +[m, C, W] = gaussian_CPD_params_given_dps(CPD, domain, evidence); + + +ns(odom) = 1; +dpsize = prod(ns(dps)); +self = domain(end); +assert(myismember(self, odom)); +self_val = evidence{self}; +T = zeros(dpsize, 1); +if length(cps) > 0 + assert(~any(isemptycell(evidence(cps)))); + cps_vals = cat(1, evidence{cps}); + for i=1:dpsize + T(i) = gaussian_prob(self_val, m(:,i) + W(:,:,i)*cps_vals, C(:,:,i)); + end +else + for i=1:dpsize + T(i) = gaussian_prob(self_val, m(:,i), C(:,:,i)); + end +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/display.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/display.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,4 @@ +function display(CPD) + +disp('gaussian_CPD object'); +disp(struct(CPD)); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/gaussian_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/gaussian_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,161 @@ +function CPD = gaussian_CPD(bnet, self, varargin) +% GAUSSIAN_CPD Make a conditional linear Gaussian distrib. +% +% CPD = gaussian_CPD(bnet, node, ...) will create a CPD with random parameters, +% where node is the number of a node in this equivalence class. + +% To define this CPD precisely, call the continuous (cts) parents (if any) X, +% the discrete parents (if any) Q, and this node Y. Then the distribution on Y is: +% - no parents: Y ~ N(mu, Sigma) +% - cts parents : Y|X=x ~ N(mu + W x, Sigma) +% - discrete parents: Y|Q=i ~ N(mu(i), Sigma(i)) +% - cts and discrete parents: Y|X=x,Q=i ~ N(mu(i) + W(i) x, Sigma(i)) +% +% The list below gives optional arguments [default value in brackets]. +% (Let ns(i) be the size of node i, X = ns(X), Y = ns(Y) and Q = prod(ns(Q)).) +% Parameters will be reshaped to the right size if necessary. +% +% mean - mu(:,i) is the mean given Q=i [ randn(Y,Q) ] +% cov - Sigma(:,:,i) is the covariance given Q=i [ repmat(100*eye(Y,Y), [1 1 Q]) ] +% weights - W(:,:,i) is the regression matrix given Q=i [ randn(Y,X,Q) ] +% cov_type - if 'diag', Sigma(:,:,i) is diagonal [ 'full' ] +% tied_cov - if 1, we constrain Sigma(:,:,i) to be the same for all i [0] +% clamp_mean - if 1, we do not adjust mu(:,i) during learning [0] +% clamp_cov - if 1, we do not adjust Sigma(:,:,i) during learning [0] +% clamp_weights - if 1, we do not adjust W(:,:,i) during learning [0] +% cov_prior_weight - weight given to I prior for estimating Sigma [0.01] +% cov_prior_entropic - if 1, we also use an entropic prior for Sigma [0] +% +% e.g., CPD = gaussian_CPD(bnet, i, 'mean', [0; 0], 'clamp_mean', 1) + +if nargin==0 + % This occurs if we are trying to load an object from a file. + CPD = init_fields; + clamp = 0; + CPD = class(CPD, 'gaussian_CPD', generic_CPD(clamp)); + return; +elseif isa(bnet, 'gaussian_CPD') + % This might occur if we are copying an object. + CPD = bnet; + return; +end +CPD = init_fields; + +CPD = class(CPD, 'gaussian_CPD', generic_CPD(0)); + +args = varargin; +ns = bnet.node_sizes; +ps = parents(bnet.dag, self); +dps = myintersect(ps, bnet.dnodes); +cps = myintersect(ps, bnet.cnodes); +fam_sz = ns([ps self]); + +CPD.self = self; +CPD.sizes = fam_sz; + +% Figure out which (if any) of the parents are discrete, and which cts, and how big they are +% dps = discrete parents, cps = cts parents +CPD.cps = find_equiv_posns(cps, ps); % cts parent index +CPD.dps = find_equiv_posns(dps, ps); +ss = fam_sz(end); +psz = fam_sz(1:end-1); +dpsz = prod(psz(CPD.dps)); +cpsz = sum(psz(CPD.cps)); + +% set default params +CPD.mean = randn(ss, dpsz); +CPD.cov = 100*repmat(eye(ss), [1 1 dpsz]); +CPD.weights = randn(ss, cpsz, dpsz); +CPD.cov_type = 'full'; +CPD.tied_cov = 0; +CPD.clamped_mean = 0; +CPD.clamped_cov = 0; +CPD.clamped_weights = 0; +CPD.cov_prior_weight = 0.01; +CPD.cov_prior_entropic = 0; +nargs = length(args); +if nargs > 0 + CPD = set_fields(CPD, args{:}); +end + +% Make sure the matrices have 1 dimension per discrete parent. +% Bug fix due to Xuejing Sun 3/6/01 +CPD.mean = myreshape(CPD.mean, [ss ns(dps)]); +CPD.cov = myreshape(CPD.cov, [ss ss ns(dps)]); +CPD.weights = myreshape(CPD.weights, [ss cpsz ns(dps)]); + +% Precompute indices into block structured matrices +% to speed up CPD_to_lambda_msg and CPD_to_pi +cpsizes = CPD.sizes(CPD.cps); +CPD.cps_block_ndx = cell(1, length(cps)); +for i=1:length(cps) + CPD.cps_block_ndx{i} = block(i, cpsizes); +end + +%%%%%%%%%%% +% Learning stuff + +% expected sufficient statistics +CPD.Wsum = zeros(dpsz,1); +CPD.WYsum = zeros(ss, dpsz); +CPD.WXsum = zeros(cpsz, dpsz); +CPD.WYYsum = zeros(ss, ss, dpsz); +CPD.WXXsum = zeros(cpsz, cpsz, dpsz); +CPD.WXYsum = zeros(cpsz, ss, dpsz); + +% For BIC +CPD.nsamples = 0; +switch CPD.cov_type + case 'full', + % since symmetric + %ncov_params = ss*(ss-1)/2; + ncov_params = ss*(ss+1)/2; + case 'diag', + ncov_params = ss; + otherwise + error(['unrecognized cov_type ' cov_type]); +end +% params = weights + mean + cov +if CPD.tied_cov + CPD.nparams = ss*cpsz*dpsz + ss*dpsz + ncov_params; +else + CPD.nparams = ss*cpsz*dpsz + ss*dpsz + dpsz*ncov_params; +end + +% for speeding up maximize_params +CPD.useC = exist('rep_mult'); + +clamped = CPD.clamped_mean & CPD.clamped_cov & CPD.clamped_weights; +CPD = set_clamped(CPD, clamped); + +%%%%%%%%%%% + +function CPD = init_fields() +% This ensures we define the fields in the same order +% no matter whether we load an object from a file, +% or create it from scratch. (Matlab requires this.) + +CPD.self = []; +CPD.sizes = []; +CPD.cps = []; +CPD.dps = []; +CPD.mean = []; +CPD.cov = []; +CPD.weights = []; +CPD.clamped_mean = []; +CPD.clamped_cov = []; +CPD.clamped_weights = []; +CPD.cov_type = []; +CPD.tied_cov = []; +CPD.Wsum = []; +CPD.WYsum = []; +CPD.WXsum = []; +CPD.WYYsum = []; +CPD.WXXsum = []; +CPD.WXYsum = []; +CPD.nsamples = []; +CPD.nparams = []; +CPD.cov_prior_weight = []; +CPD.cov_prior_entropic = []; +CPD.useC = []; +CPD.cps_block_ndx = []; diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/gaussian_CPD_params_given_dps.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/gaussian_CPD_params_given_dps.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,28 @@ +function [m, C, W] = gaussian_CPD_params_given_dps(CPD, domain, evidence) +% GAUSSIAN_CPD_PARAMS_GIVEN_EV_ON_DPS Extract parameters given evidence on all discrete parents +% function [m, C, W] = gaussian_CPD_params_given_ev_on_dps(CPD, domain, evidence) + +ps = domain(1:end-1); +dps = ps(CPD.dps); +if isempty(dps) + m = CPD.mean; + C = CPD.cov; + W = CPD.weights; +else + odom = domain(~isemptycell(evidence(domain))); + dops = myintersect(dps, odom); + dpvals = cat(1, evidence{dops}); + if length(dops) == length(dps) + dpsizes = CPD.sizes(CPD.dps); + dpval = subv2ind(dpsizes, dpvals(:)'); + m = CPD.mean(:, dpval); + C = CPD.cov(:, :, dpval); + W = CPD.weights(:, :, dpval); + else + map = find_equiv_posns(dops, dps); + index = mk_multi_index(length(dps), map, dpvals); + m = CPD.mean(:, index{:}); + C = CPD.cov(:, :, index{:}); + W = CPD.weights(:, :, index{:}); + end +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/get_field.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/get_field.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,19 @@ +function val = get_params(CPD, name) +% GET_PARAMS Get the parameters (fields) for a gaussian_CPD object +% val = get_params(CPD, name) +% +% The following fields can be accessed +% +% mean - mu(:,i) is the mean given Q=i +% cov - Sigma(:,:,i) is the covariance given Q=i +% weights - W(:,:,i) is the regression matrix given Q=i +% +% e.g., mean = get_params(CPD, 'mean') + +switch name + case 'mean', val = CPD.mean; + case 'cov', val = CPD.cov; + case 'weights', val = CPD.weights; + otherwise, + error(['invalid argument name ' name]); +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/learn_params.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/learn_params.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,31 @@ +function CPD = learn_params(CPD, fam, data, ns, cnodes) +%function CPD = learn_params(CPD, fam, data, ns, cnodes) +% LEARN_PARAMS Compute the maximum likelihood estimate of the params of a gaussian CPD given complete data +% CPD = learn_params(CPD, fam, data, ns, cnodes) +% +% data(i,m) is the value of node i in case m (can be cell array). +% We assume this node has a maximize_params method. + +ncases = size(data, 2); +CPD = reset_ess(CPD); +% make a fully observed joint distribution over the family +fmarginal.domain = fam; +fmarginal.T = 1; +fmarginal.mu = []; +fmarginal.Sigma = []; +if ~iscell(data) + cases = num2cell(data); +else + cases = data; +end +hidden_bitv = zeros(1, max(fam)); +for m=1:ncases + % specify (as a bit vector) which elements in the family domain are hidden + hidden_bitv = zeros(1, max(fmarginal.domain)); + ev = cases(:,m); + hidden_bitv(find(isempty(ev)))=1; + CPD = update_ess(CPD, fmarginal, ev, ns, cnodes, hidden_bitv); +end +CPD = maximize_params(CPD); + + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/log_prob_node.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/log_prob_node.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,49 @@ +function L = log_prob_node(CPD, self_ev, pev) +% LOG_PROB_NODE Compute prod_m log P(x(i,m)| x(pi_i,m), theta_i) for node i (gaussian) +% L = log_prob_node(CPD, self_ev, pev) +% +% self_ev(m) is the evidence on this node in case m. +% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents). +% (These may also be cell arrays.) + +if iscell(self_ev), usecell = 1; else usecell = 0; end + +use_log = 1; +ncases = length(self_ev); +nparents = length(CPD.sizes)-1; +assert(ncases == size(pev, 2)); + +if ncases == 0 + L = 0; + return; +end + +L = 0; +for m=1:ncases + if isempty(CPD.dps) + i = 1; + else + if usecell + dpvals = cat(1, pev{CPD.dps, m}); + else + dpvals = pev(CPD.dps, m); + end + i = subv2ind(CPD.sizes(CPD.dps), dpvals(:)'); + end + if usecell + y = self_ev{m}; + else + y = self_ev(m); + end + if length(CPD.cps) == 0 + L = L + gaussian_prob(y, CPD.mean(:,i), CPD.cov(:,:,i), use_log); + else + if usecell + x = cat(1, pev{CPD.cps, m}); + else + x = pev(CPD.cps, m); + end + L = L + gaussian_prob(y, CPD.mean(:,i) + CPD.weights(:,:,i)*x, CPD.cov(:,:,i), use_log); + end +end + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/maximize_params.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/maximize_params.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,68 @@ +function CPD = maximize_params(CPD, temp) +% MAXIMIZE_PARAMS Set the params of a CPD to their ML values (Gaussian) +% CPD = maximize_params(CPD, temperature) +% +% Temperature is currently ignored. + +if ~adjustable_CPD(CPD), return; end + + +if CPD.clamped_mean + cl_mean = CPD.mean; +else + cl_mean = []; +end + +if CPD.clamped_cov + cl_cov = CPD.cov; +else + cl_cov = []; +end + +if CPD.clamped_weights + cl_weights = CPD.weights; +else + cl_weights = []; +end + +[ssz psz Q] = size(CPD.weights); + +[ss cpsz dpsz] = size(CPD.weights); % ss = self size = ssz +if cpsz > CPD.nsamples + fprintf('gaussian_CPD/maximize_params: warning: input dimension (%d) > nsamples (%d)\n', ... + cpsz, CPD.nsamples); +end + +prior = repmat(CPD.cov_prior_weight*eye(ssz,ssz), [1 1 Q]); + + +[CPD.mean, CPD.cov, CPD.weights] = ... + clg_Mstep(CPD.Wsum, CPD.WYsum, CPD.WYYsum, [], CPD.WXsum, CPD.WXXsum, CPD.WXYsum, ... + 'cov_type', CPD.cov_type, 'clamped_mean', cl_mean, ... + 'clamped_cov', cl_cov, 'clamped_weights', cl_weights, ... + 'tied_cov', CPD.tied_cov, ... + 'cov_prior', prior); + +if 0 +CPD.mean = reshape(CPD.mean, [ss dpsz]); +CPD.cov = reshape(CPD.cov, [ss ss dpsz]); +CPD.weights = reshape(CPD.weights, [ss cpsz dpsz]); +end + +% Bug fix 11 May 2003 KPM +% clg_Mstep collapses all discrete parents into one mega-node +% but convert_to_CPT needs access to each parent separately +sz = CPD.sizes; +ss = sz(end); + +% Bug fix KPM 20 May 2003: +cpsz = sum(sz(CPD.cps)); +%if isempty(CPD.cps) +% cpsz = 0; +%else +% cpsz = sz(CPD.cps); +%end +dpsz = sz(CPD.dps); +CPD.mean = myreshape(CPD.mean, [ss dpsz]); +CPD.cov = myreshape(CPD.cov, [ss ss dpsz]); +CPD.weights = myreshape(CPD.weights, [ss cpsz dpsz]); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/maximize_params_debug.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/maximize_params_debug.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,189 @@ +function CPD = maximize_params(CPD, temp) +% MAXIMIZE_PARAMS Set the params of a CPD to their ML values (Gaussian) +% CPD = maximize_params(CPD, temperature) +% +% Temperature is currently ignored. + +if ~adjustable_CPD(CPD), return; end + +CPD1 = struct(new_maximize_params(CPD)); +CPD2 = struct(old_maximize_params(CPD)); +assert(approxeq(CPD1.mean, CPD2.mean)) +assert(approxeq(CPD1.cov, CPD2.cov)) +assert(approxeq(CPD1.weights, CPD2.weights)) + +CPD = new_maximize_params(CPD); + +%%%%%%% +function CPD = new_maximize_params(CPD) + +if CPD.clamped_mean + cl_mean = CPD.mean; +else + cl_mean = []; +end + +if CPD.clamped_cov + cl_cov = CPD.cov; +else + cl_cov = []; +end + +if CPD.clamped_weights + cl_weights = CPD.weights; +else + cl_weights = []; +end + +[ssz psz Q] = size(CPD.weights); + +prior = repmat(CPD.cov_prior_weight*eye(ssz,ssz), [1 1 Q]); +[CPD.mean, CPD.cov, CPD.weights] = ... + Mstep_clg('w', CPD.Wsum, 'YY', CPD.WYYsum, 'Y', CPD.WYsum, 'YTY', [], ... + 'XX', CPD.WXXsum, 'XY', CPD.WXYsum, 'X', CPD.WXsum, ... + 'cov_type', CPD.cov_type, 'clamped_mean', cl_mean, ... + 'clamped_cov', cl_cov, 'clamped_weights', cl_weights, ... + 'tied_cov', CPD.tied_cov, ... + 'cov_prior', prior); + + +%%%%%%%%%%% + +function CPD = old_maximize_params(CPD) + + +if ~adjustable_CPD(CPD), return; end + +%assert(approxeq(CPD.nsamples, sum(CPD.Wsum))); +assert(~any(isnan(CPD.WXXsum))) +assert(~any(isnan(CPD.WXYsum))) +assert(~any(isnan(CPD.WYYsum))) + +[self_size cpsize dpsize] = size(CPD.weights); + +% Append 1s to the parents, and derive the corresponding cross products. +% This is used when estimate the means and weights simultaneosuly, +% and when estimatting Sigma. +% Let x2 = [x 1]' +XY = zeros(cpsize+1, self_size, dpsize); % XY(:,:,i) = sum_l w(l,i) x2(l) y(l)' +XX = zeros(cpsize+1, cpsize+1, dpsize); % XX(:,:,i) = sum_l w(l,i) x2(l) x2(l)' +YY = zeros(self_size, self_size, dpsize); % YY(:,:,i) = sum_l w(l,i) y(l) y(l)' +for i=1:dpsize + XY(:,:,i) = [CPD.WXYsum(:,:,i) % X*Y + CPD.WYsum(:,i)']; % 1*Y + % [x * [x' 1] = [xx' x + % 1] x' 1] + XX(:,:,i) = [CPD.WXXsum(:,:,i) CPD.WXsum(:,i); + CPD.WXsum(:,i)' CPD.Wsum(i)]; + YY(:,:,i) = CPD.WYYsum(:,:,i); +end + +w = CPD.Wsum(:); +% Set any zeros to one before dividing +% This is valid because w(i)=0 => WYsum(:,i)=0, etc +w = w + (w==0); + +if CPD.clamped_mean + % Estimating B2 and then setting the last column (the mean) to the clamped mean is *not* equivalent + % to estimating B and then adding the clamped_mean to the last column. + if ~CPD.clamped_weights + B = zeros(self_size, cpsize, dpsize); + for i=1:dpsize + if det(CPD.WXXsum(:,:,i))==0 + B(:,:,i) = 0; + else + % Eqn 9 in table 2 of TR + %B(:,:,i) = CPD.WXYsum(:,:,i)' * inv(CPD.WXXsum(:,:,i)); + B(:,:,i) = (CPD.WXXsum(:,:,i) \ CPD.WXYsum(:,:,i))'; + end + end + %CPD.weights = reshape(B, [self_size cpsize dpsize]); + CPD.weights = B; + end +elseif CPD.clamped_weights % KPM 1/25/02 + if ~CPD.clamped_mean % ML estimate is just sample mean of the residuals + for i=1:dpsize + CPD.mean(:,i) = (CPD.WYsum(:,i) - CPD.weights(:,:,i) * CPD.WXsum(:,i)) / w(i); + end + end +else % nothing is clamped, so estimate mean and weights simultaneously + B2 = zeros(self_size, cpsize+1, dpsize); + for i=1:dpsize + if det(XX(:,:,i))==0 % fix by U. Sondhauss 6/27/99 + B2(:,:,i)=0; + else + % Eqn 9 in table 2 of TR + %B2(:,:,i) = XY(:,:,i)' * inv(XX(:,:,i)); + B2(:,:,i) = (XX(:,:,i) \ XY(:,:,i))'; + end + CPD.mean(:,i) = B2(:,cpsize+1,i); + CPD.weights(:,:,i) = B2(:,1:cpsize,i); + end +end + +% Let B2 = [W mu] +if cpsize>0 + B2(:,1:cpsize,:) = reshape(CPD.weights, [self_size cpsize dpsize]); +end +B2(:,cpsize+1,:) = reshape(CPD.mean, [self_size dpsize]); + +% To avoid singular covariance matrices, +% we use the regularization method suggested in "A Quasi-Bayesian approach to estimating +% parameters for mixtures of normal distributions", Hamilton 91. +% If the ML estimate is Sigma = M/N, the MAP estimate is (M+gamma*I) / (N+gamma), +% where gamma >=0 is a smoothing parameter (equivalent sample size of I prior) + +gamma = CPD.cov_prior_weight; + +if ~CPD.clamped_cov + if CPD.cov_prior_entropic % eqn 12 of Brand AI/Stat 99 + Z = 1-temp; + % When temp > 1, Z is negative, so we are dividing by a smaller + % number, ie. increasing the variance. + else + Z = 0; + end + if CPD.tied_cov + S = zeros(self_size, self_size); + % Eqn 2 from table 2 in TR + for i=1:dpsize + S = S + (YY(:,:,i) - B2(:,:,i)*XY(:,:,i)); + end + %denom = CPD.nsamples + gamma + Z; + denom = CPD.nsamples + Z; + S = (S + gamma*eye(self_size)) / denom; + if strcmp(CPD.cov_type, 'diag') + S = diag(diag(S)); + end + CPD.cov = repmat(S, [1 1 dpsize]); + else + for i=1:dpsize + % Eqn 1 from table 2 in TR + S = YY(:,:,i) - B2(:,:,i)*XY(:,:,i); + %denom = w(i) + gamma + Z; + denom = w(i) + Z; + S = (S + gamma*eye(self_size)) / denom; + CPD.cov(:,:,i) = S; + end + if strcmp(CPD.cov_type, 'diag') + for i=1:dpsize + CPD.cov(:,:,i) = diag(diag(CPD.cov(:,:,i))); + end + end + end +end + + +check_covars = 0; +min_covar = 1e-5; +if check_covars % prevent collapsing to a point + for i=1:dpsize + if min(svd(CPD.cov(:,:,i))) < min_covar + disp(['resetting singular covariance for node ' num2str(CPD.self)]); + CPD.cov(:,:,i) = CPD.init_cov(:,:,i); + end + end +end + + + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/private/CPD_to_linear_gaussian.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/private/CPD_to_linear_gaussian.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,19 @@ +function [mu, Sigma, W] = CPD_to_linear_gaussian(CPD, domain, ns, cnodes, evidence) + +ps = domain(1:end-1); +dnodes = mysetdiff(1:length(ns), cnodes); +dps = myintersect(ps, dnodes); % discrete parents + +if isempty(dps) + Q = 1; +else + assert(~any(isemptycell(evidence(dps)))); + dpvals = cat(1, evidence{dps}); + Q = subv2ind(ns(dps), dpvals(:)'); +end + +mu = CPD.mean(:,Q); +Sigma = CPD.cov(:,:,Q); +W = CPD.weights(:,:,Q); + + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/private/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/private/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,2 @@ +/CPD_to_linear_gaussian.m/1.1.1.1/Wed May 29 15:59:52 2002// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/private/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/private/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@gaussian_CPD/private diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/private/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/private/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/reset_ess.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/reset_ess.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,11 @@ +function CPD = reset_ess(CPD) +% RESET_ESS Reset the Expected Sufficient Statistics for a Gaussian CPD. +% CPD = reset_ess(CPD) + +CPD.nsamples = 0; +CPD.Wsum = zeros(size(CPD.Wsum)); +CPD.WYsum = zeros(size(CPD.WYsum)); +CPD.WYYsum = zeros(size(CPD.WYYsum)); +CPD.WXsum = zeros(size(CPD.WXsum)); +CPD.WXXsum = zeros(size(CPD.WXXsum)); +CPD.WXYsum = zeros(size(CPD.WXYsum)); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/sample_node.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/sample_node.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,22 @@ +function y = sample_node(CPD, pev) +% SAMPLE_NODE Draw a random sample from P(Xi | x(pi_i), theta_i) (gaussian) +% y = sample_node(CPD, parent_evidence) +% +% pev{i} is the value of the i'th parent (if there are any parents) +% y is the sampled value (a scalar or vector) + +if length(CPD.dps)==0 + i = 1; +else + dpvals = cat(1, pev{CPD.dps}); + i = subv2ind(CPD.sizes(CPD.dps), dpvals(:)'); +end + +if length(CPD.cps) == 0 + y = gsamp(CPD.mean(:,i), CPD.cov(:,:,i), 1); +else + pev = pev(:); + x = cat(1, pev{CPD.cps}); + y = gsamp(CPD.mean(:,i) + CPD.weights(:,:,i)*x(:), CPD.cov(:,:,i), 1); +end +y = y(:); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/set_fields.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/set_fields.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,43 @@ +function CPD = set_fields(CPD, varargin) +% SET_PARAMS Set the parameters (fields) for a gaussian_CPD object +% CPD = set_params(CPD, name/value pairs) +% +% The following optional arguments can be specified in the form of name/value pairs: +% +% mean - mu(:,i) is the mean given Q=i +% cov - Sigma(:,:,i) is the covariance given Q=i +% weights - W(:,:,i) is the regression matrix given Q=i +% cov_type - if 'diag', Sigma(:,:,i) is diagonal +% tied_cov - if 1, we constrain Sigma(:,:,i) to be the same for all i +% clamp_mean - if 1, we do not adjust mu(:,i) during learning +% clamp_cov - if 1, we do not adjust Sigma(:,:,i) during learning +% clamp_weights - if 1, we do not adjust W(:,:,i) during learning +% clamp - if 1, we do not adjust any params +% cov_prior_weight - weight given to I prior for estimating Sigma +% cov_prior_entropic - if 1, we also use an entropic prior for Sigma [0] +% +% e.g., CPD = set_params(CPD, 'mean', [0;0]) + +args = varargin; +nargs = length(args); +for i=1:2:nargs + switch args{i}, + case 'mean', CPD.mean = args{i+1}; + case 'cov', CPD.cov = args{i+1}; + case 'weights', CPD.weights = args{i+1}; + case 'cov_type', CPD.cov_type = args{i+1}; + %case 'tied_cov', CPD.tied_cov = strcmp(args{i+1}, 'yes'); + case 'tied_cov', CPD.tied_cov = args{i+1}; + case 'clamp_mean', CPD.clamped_mean = args{i+1}; + case 'clamp_cov', CPD.clamped_cov = args{i+1}; + case 'clamp_weights', CPD.clamped_weights = args{i+1}; + case 'clamp', clamp = args{i+1}; + CPD.clamped_mean = clamp; + CPD.clamped_cov = clamp; + CPD.clamped_weights = clamp; + case 'cov_prior_weight', CPD.cov_prior_weight = args{i+1}; + case 'cov_prior_entropic', CPD.cov_prior_entropic = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gaussian_CPD/update_ess.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gaussian_CPD/update_ess.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,88 @@ +function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv) +% UPDATE_ESS Update the Expected Sufficient Statistics of a Gaussian node +% function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv) + +%if nargin < 6 +% hidden_bitv = zeros(1, max(fmarginal.domain)); +% hidden_bitv(find(isempty(evidence)))=1; +%end + +dom = fmarginal.domain; +self = dom(end); +ps = dom(1:end-1); +cps = myintersect(ps, cnodes); +dps = mysetdiff(ps, cps); + +CPD.nsamples = CPD.nsamples + 1; +[ss cpsz dpsz] = size(CPD.weights); % ss = self size +[ss dpsz] = size(CPD.mean); + +% Let X be the cts parent (if any), Y be the cts child (self). + +if ~hidden_bitv(self) & ~any(hidden_bitv(cps)) & all(hidden_bitv(dps)) + % Speedup for the common case that all cts nodes are observed, all discrete nodes are hidden + % Since X and Y are observed, SYY = 0, SXX = 0, SXY = 0 + % Since discrete parents are hidden, we do not need to add evidence to w. + w = fmarginal.T(:); + CPD.Wsum = CPD.Wsum + w; + y = evidence{self}; + Cyy = y*y'; + if ~CPD.useC + WY = repmat(w(:)',ss,1); % WY(y,i) = w(i) + WYY = repmat(reshape(WY, [ss 1 dpsz]), [1 ss 1]); % WYY(y,y',i) = w(i) + %CPD.WYsum = CPD.WYsum + WY .* repmat(y(:), 1, dpsz); + CPD.WYsum = CPD.WYsum + y(:) * w(:)'; + CPD.WYYsum = CPD.WYYsum + WYY .* repmat(reshape(Cyy, [ss ss 1]), [1 1 dpsz]); + else + W = w(:)'; + W2 = reshape(W, [1 1 dpsz]); + CPD.WYsum = CPD.WYsum + rep_mult(W, y(:), size(CPD.WYsum)); + CPD.WYYsum = CPD.WYYsum + rep_mult(W2, Cyy, size(CPD.WYYsum)); + end + if cpsz > 0 % X exists + x = cat(1, evidence{cps}); x = x(:); + Cxx = x*x'; + Cxy = x*y'; + WX = repmat(w(:)',cpsz,1); % WX(x,i) = w(i) + WXX = repmat(reshape(WX, [cpsz 1 dpsz]), [1 cpsz 1]); % WXX(x,x',i) = w(i) + WXY = repmat(reshape(WX, [cpsz 1 dpsz]), [1 ss 1]); % WXY(x,y,i) = w(i) + if ~CPD.useC + CPD.WXsum = CPD.WXsum + WX .* repmat(x(:), 1, dpsz); + CPD.WXXsum = CPD.WXXsum + WXX .* repmat(reshape(Cxx, [cpsz cpsz 1]), [1 1 dpsz]); + CPD.WXYsum = CPD.WXYsum + WXY .* repmat(reshape(Cxy, [cpsz ss 1]), [1 1 dpsz]); + else + CPD.WXsum = CPD.WXsum + rep_mult(W, x(:), size(CPD.WXsum)); + CPD.WXXsum = CPD.WXXsum + rep_mult(W2, Cxx, size(CPD.WXXsum)); + CPD.WXYsum = CPD.WXYsum + rep_mult(W2, Cxy, size(CPD.WXYsum)); + end + end + return; +end + +% general (non-vectorized) case +fullm = add_evidence_to_gmarginal(fmarginal, evidence, ns, cnodes); % slow! + +if dpsz == 1 % no discrete parents + w = 1; +else + w = fullm.T(:); +end + +CPD.Wsum = CPD.Wsum + w; +xi = 1:cpsz; +yi = (cpsz+1):(cpsz+ss); +for i=1:dpsz + muY = fullm.mu(yi, i); + SYY = fullm.Sigma(yi, yi, i); + CPD.WYsum(:,i) = CPD.WYsum(:,i) + w(i)*muY; + CPD.WYYsum(:,:,i) = CPD.WYYsum(:,:,i) + w(i)*(SYY + muY*muY'); % E[X Y] = Cov[X,Y] + E[X] E[Y] + if cpsz > 0 + muX = fullm.mu(xi, i); + SXX = fullm.Sigma(xi, xi, i); + SXY = fullm.Sigma(xi, yi, i); + CPD.WXsum(:,i) = CPD.WXsum(:,i) + w(i)*muX; + CPD.WXXsum(:,:,i) = CPD.WXXsum(:,:,i) + w(i)*(SXX + muX*muX'); + CPD.WXYsum(:,:,i) = CPD.WXYsum(:,:,i) + w(i)*(SXY + muX*muY'); + end +end + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,8 @@ +/README/1.1.1.1/Wed May 29 15:59:52 2002// +/adjustable_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002// +/display.m/1.1.1.1/Wed May 29 15:59:52 2002// +/generic_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002// +/learn_params.m/1.1.1.1/Thu Jun 10 01:53:20 2004// +/log_prior.m/1.1.1.1/Wed May 29 15:59:52 2002// +/set_clamped.m/1.1.1.1/Wed May 29 15:59:52 2002// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/CVS/Entries.Log --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/CVS/Entries.Log Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +A D/Old//// diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@generic_CPD diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/Old/BIC_score_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/Old/BIC_score_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,26 @@ +function score = BIC_score_CPD(CPD, fam, data, ns, cnodes) +% BIC_score_CPD Compute the BIC score of a generic CPD +% score = BIC_score_CPD(CPD, fam, data, ns, cnodes) +% +% We assume this node has a maximize_params method + +ncases = size(data, 2); +CPD = reset_ess(CPD); +% make a fully observed joint distribution over the family +fmarginal.domain = fam; +fmarginal.T = 1; +fmarginal.mu = []; +fmarginal.Sigma = []; +if ~iscell(data) + cases = num2cell(data); +else + cases = data; +end +for m=1:ncases + CPD = update_ess(CPD, fmarginal, cases(:,m), ns, cnodes); +end +CPD = maximize_params(CPD); +self = fam(end); +ps = fam(1:end-1); +L = log_prob_node(CPD, cases(self,:), cases(ps,:)); +score = L - 0.5*CPD.nparams*log(ncases); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/Old/CPD_to_dpots.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/Old/CPD_to_dpots.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,16 @@ +function pots = CPD_to_dpots(CPD, domain, ns, cnodes, evidence) +% CPD_TO_DPOTS Convert the CPD to several discrete potentials, for different instantiations (generic) +% pots = CPD_to_dpots(CPD, domain, ns, cnodes, evidence) +% +% domain(:,i) is the domain of the i'th instantiation of CPD. +% node_sizes(i) is the size of node i. +% cnodes = all the cts nodes +% evidence{i} is the evidence on the i'th node. +% +% This just calls CPD_to_dpot for each domain. + +nCPDs = size(domain,2); +pots = cell(1,nCPDs); +for i=1:nCPDs + pots{i} = CPD_to_dpot(CPD, domain(:,i), ns, cnodes, evidence); +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/Old/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/Old/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,3 @@ +/BIC_score_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002// +/CPD_to_dpots.m/1.1.1.1/Wed May 29 15:59:52 2002// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/Old/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/Old/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@generic_CPD/Old diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/Old/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/Old/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/README --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/README Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,2 @@ +A generic CPD implements general purpose functions like 'display', +that subtypes can inherit. diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/adjustable_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/adjustable_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,5 @@ +function p = adjustable_CPD(CPD) +% ADJUSTABLE_CPD Does this CPD have any adjustable params? (generic) +% p = adjustable_CPD(CPD) + +p = ~CPD.clamped; diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/display.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/display.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,3 @@ +function display(CPD) + +disp(struct(CPD)); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/generic_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/generic_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,8 @@ +function CPD = generic_CPD(clamped) +% GENERIC_CPD Virtual constructor for generic CPD +% CPD = discrete_CPD(clamped) + +if nargin < 1, clamped = 0; end + +CPD.clamped = clamped; +CPD = class(CPD, 'generic_CPD'); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/learn_params.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/learn_params.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,32 @@ +function CPD = learn_params(CPD, fam, data, ns, cnodes) +% LEARN_PARAMS Compute the maximum likelihood estimate of the params of a generic CPD given complete data +% CPD = learn_params(CPD, fam, data, ns, cnodes) +% +% data(i,m) is the value of node i in case m (can be cell array). +% We assume this node has a maximize_params method. + +%error('no longer supported') % KPM 1 Feb 03 + +if 1 +ncases = size(data, 2); +CPD = reset_ess(CPD); +% make a fully observed joint distribution over the family +fmarginal.domain = fam; +fmarginal.T = 1; +fmarginal.mu = []; +fmarginal.Sigma = []; +if ~iscell(data) + cases = num2cell(data); +else + cases = data; +end +hidden_bitv = zeros(1, max(fam)); +for m=1:ncases + % specify (as a bit vector) which elements in the family domain are hidden + hidden_bitv = zeros(1, max(fmarginal.domain)); + ev = cases(:,m); + hidden_bitv(find(isempty(evidence)))=1; + CPD = update_ess(CPD, fmarginal, ev, ns, cnodes, hidden_bitv); +end +CPD = maximize_params(CPD); +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/log_prior.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/log_prior.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,5 @@ +function L = log_prior(CPD) +% LOG_PRIOR Return log P(theta) for a generic CPD - we return 0 +% L = log_prior(CPD) + +L = 0; diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@generic_CPD/set_clamped.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@generic_CPD/set_clamped.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,3 @@ +function CPD = set_clamped(CPD, bit) + +CPD.clamped = bit; diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/CPD_to_lambda_msg.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/CPD_to_lambda_msg.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,62 @@ +function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence) +% CPD_TO_LAMBDA_MSG Compute lambda message (gmux) +% lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence) +% Pearl p183 eq 4.52 + +% Let Y be this node, X1..Xn be the cts parents and M the discrete switch node. +% e.g., for n=3, M=1 +% +% X1 X2 X3 M +% \ +% \ +% Y +% +% So the only case in which we send an informative message is if p=1=M. +% To the other cts parents, we send the "know nothing" message. + +switch msg_type + case 'd', + error('gaussian_CPD can''t create discrete msgs') + case 'g', + cps = ps(CPD.cps); + cpsizes = CPD.sizes(CPD.cps); + self_size = CPD.sizes(end); + i = find_equiv_posns(p, cps); % p is n's i'th cts parent + psz = cpsizes(i); + dps = ps(CPD.dps); + M = evidence{dps}; + if isempty(M) + error('gmux node must have observed discrete parent') + end + P = msg{n}.lambda.precision; + if all(P == 0) | (cps(M) ~= p) % if we know nothing, or are sending to a disconnected parent + lam_msg.precision = zeros(psz, psz); + lam_msg.info_state = zeros(psz, 1); + return; + end + % We are sending a message to the only effectively connected parent. + % There are no other incoming pi messages. + Bmu = CPD.mean(:,M); + BSigma = CPD.cov(:,:,M); + Bi = CPD.weights(:,:,M); + if (det(P) > 0) | isinf(P) + if isinf(P) % Y is observed + Sigma_lambda = zeros(self_size, self_size); % infinite precision => 0 variance + mu_lambda = msg{n}.lambda.mu; % observed_value; + else + Sigma_lambda = inv(P); + mu_lambda = Sigma_lambda * msg{n}.lambda.info_state; + end + C = inv(Sigma_lambda + BSigma); + lam_msg.precision = Bi' * C * Bi; + lam_msg.info_state = Bi' * C * (mu_lambda - Bmu); + else + % method that uses matrix inversion lemma + A = inv(P + inv(BSigma)); + C = P - P*A*P; + lam_msg.precision = Bi' * C * Bi; + D = eye(self_size) - P*A; + z = msg{n}.lambda.info_state; + lam_msg.info_state = Bi' * (D*z - D*P*Bmu); + end +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/CPD_to_pi.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/CPD_to_pi.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,18 @@ +function pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence) +% CPD_TO_PI Compute the pi vector (gaussian) +% function pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence) + +switch msg_type + case 'd', + error('gaussian_CPD can''t create discrete msgs') + case 'g', + dps = ps(CPD.dps); + k = evidence{dps}; + if isempty(k) + error('gmux node must have observed discrete parent') + end + m = msg{n}.pi_from_parent{k}; + B = CPD.weights(:,:,k); + pi.mu = CPD.mean(:,k) + B * m.mu; + pi.Sigma = CPD.cov(:,:,k) + B * m.Sigma * B'; +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,7 @@ +/CPD_to_lambda_msg.m/1.1.1.1/Wed May 29 15:59:52 2002// +/CPD_to_pi.m/1.1.1.1/Wed May 29 15:59:54 2002// +/convert_to_pot.m/1.1.1.1/Wed May 29 15:59:52 2002// +/display.m/1.1.1.1/Wed May 29 15:59:54 2002// +/gmux_CPD.m/1.1.1.1/Wed May 29 15:59:54 2002// +/sample_node.m/1.1.1.1/Wed May 29 15:59:54 2002// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/CVS/Entries.Log --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/CVS/Entries.Log Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +A D/Old//// diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@gmux_CPD diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/Old/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/Old/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,2 @@ +/gmux_CPD.m/1.1.1.1/Wed May 29 15:59:54 2002// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/Old/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/Old/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@gmux_CPD/Old diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/Old/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/Old/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/Old/gmux_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/Old/gmux_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,92 @@ +function CPD = gmux_CPD(bnet, self, varargin) +% GMUX_CPD Make a Gaussian multiplexer node +% +% CPD = gmux_CPD(bnet, node, ...) is used similarly to gaussian_CPD, +% except we assume there is exactly one discrete parent (call it M) +% which is used to select which cts parent to pass through to the output. +% i.e., we define P(Y=y|M=m, X1, ..., XK) = N(y | W*x(m) + mu, Sigma) +% where Y represents this node, and the Xi's are the cts parents. +% All the Xi must have the same size, and the num values for M must be K. +% +% Currently the params for this kind of CPD cannot be learned. +% +% Optional arguments [ default in brackets ] +% +% mean - mu [zeros(Y,1)] +% cov - Sigma [eye(Y,Y)] +% weights - W [ randn(Y,X) ] + +if nargin==0 + % This occurs if we are trying to load an object from a file. + CPD = init_fields; + clamp = 0; + CPD = class(CPD, 'gmux_CPD', generic_CPD(clamp)); + return; +elseif isa(bnet, 'gmux_CPD') + % This might occur if we are copying an object. + CPD = bnet; + return; +end +CPD = init_fields; + +CPD = class(CPD, 'gmux_CPD', generic_CPD(1)); + +ns = bnet.node_sizes; +ps = parents(bnet.dag, self); +dps = myintersect(ps, bnet.dnodes); +cps = myintersect(ps, bnet.cnodes); +fam_sz = ns([ps self]); + +CPD.self = self; +CPD.sizes = fam_sz; + +% Figure out which (if any) of the parents are discrete, and which cts, and how big they are +% dps = discrete parents, cps = cts parents +CPD.cps = find_equiv_posns(cps, ps); % cts parent index +CPD.dps = find_equiv_posns(dps, ps); +if length(CPD.dps) ~= 1 + error('gmux must have exactly 1 discrete parent') +end +ss = fam_sz(end); +cpsz = fam_sz(CPD.cps(1)); % in gaussian_CPD, cpsz = sum(fam_sz(CPD.cps)) +if ~all(fam_sz(CPD.cps) == cpsz) + error('all cts parents must have same size') +end +dpsz = fam_sz(CPD.dps); +if dpsz ~= length(cps) + error(['the arity of the mux node is ' num2str(dpsz) ... + ' but there are ' num2str(length(cps)) ' cts parents']); +end + +% set default params +CPD.mean = zeros(ss, 1); +CPD.cov = eye(ss); +CPD.weights = randn(ss, cpsz); + +args = varargin; +nargs = length(args); +for i=1:2:nargs + switch args{i}, + case 'mean', CPD.mean = args{i+1}; + case 'cov', CPD.cov = args{i+1}; + case 'weights', CPD.weights = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end +end + +%%%%%%%%%%% + +function CPD = init_fields() +% This ensures we define the fields in the same order +% no matter whether we load an object from a file, +% or create it from scratch. (Matlab requires this.) + +CPD.self = []; +CPD.sizes = []; +CPD.cps = []; +CPD.dps = []; +CPD.mean = []; +CPD.cov = []; +CPD.weights = []; + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/convert_to_pot.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/convert_to_pot.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,37 @@ +function pot = convert_to_pot(CPD, pot_type, domain, evidence) +% CONVERT_TO_POT Convert a gmux CPD to a Gaussian potential +% pot = convert_to_pot(CPD, pot_type, domain, evidence) + +switch pot_type + case {'d', 'u', 'cg', 'scg'}, + error(['can''t convert gmux to potential of type ' pot_type]) + + case {'c','g'}, + % We create a large weight matrix with zeros in all blocks corresponding + % to the non-chosen parents, since they are effectively disconnected. + % The chosen parent is determined by the value, m, of the discrete parent. + % Thus the potential is as large as the whole family. + ps = domain(1:end-1); + dps = ps(CPD.dps); % CPD.dps is an index, not a node number (because of param tying) + cps = ps(CPD.cps); + m = evidence{dps}; + if isempty(m) + error('gmux node must have observed discrete parent') + end + bs = CPD.sizes(CPD.cps); + b = block(m, bs); + sum_cpsz = sum(CPD.sizes(CPD.cps)); + selfsz = CPD.sizes(end); + W = zeros(selfsz, sum_cpsz); + W(:,b) = CPD.weights(:,:,m); + + ns = zeros(1, max(domain)); + ns(domain) = CPD.sizes; + self = domain(end); + cdom = [cps(:)' self]; + pot = linear_gaussian_to_cpot(CPD.mean(:,m), CPD.cov(:,:,m), W, domain, ns, cdom, evidence); + + otherwise, + error(['unrecognized pot_type' pot_type]) +end + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/display.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/display.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,4 @@ +function display(CPD) + +disp('gmux_CPD object'); +disp(struct(CPD)); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/gmux_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/gmux_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,95 @@ +function CPD = gmux_CPD(bnet, self, varargin) +% GMUX_CPD Make a Gaussian multiplexer node +% +% CPD = gmux_CPD(bnet, node, ...) is used similarly to gaussian_CPD, +% except we assume there is exactly one discrete parent (call it M) +% which is used to select which cts parent to pass through to the output. +% i.e., we define P(Y=y|M=m, X1, ..., XK) = N(y | W(m)*x(m) + mu(m), Sigma(m)) +% where Y represents this node, and the Xi's are the cts parents. +% All the Xi must have the same size, and the num values for M must be K. +% +% Currently the params for this kind of CPD cannot be learned. +% +% Optional arguments [ default in brackets ] +% +% mean - mu(:,i) is the mean given M=i [ zeros(Y,K) ] +% cov - Sigma(:,:,i) is the covariance given M=i [ repmat(1*eye(Y,Y), [1 1 K]) ] +% weights - W(:,:,i) is the regression matrix given M=i [ randn(Y,X,K) ] + +if nargin==0 + % This occurs if we are trying to load an object from a file. + CPD = init_fields; + clamp = 0; + CPD = class(CPD, 'gmux_CPD', generic_CPD(clamp)); + return; +elseif isa(bnet, 'gmux_CPD') + % This might occur if we are copying an object. + CPD = bnet; + return; +end +CPD = init_fields; + +CPD = class(CPD, 'gmux_CPD', generic_CPD(1)); + +ns = bnet.node_sizes; +ps = parents(bnet.dag, self); +dps = myintersect(ps, bnet.dnodes); +cps = myintersect(ps, bnet.cnodes); +fam_sz = ns([ps self]); + +CPD.self = self; +CPD.sizes = fam_sz; + +% Figure out which (if any) of the parents are discrete, and which cts, and how big they are +% dps = discrete parents, cps = cts parents +CPD.cps = find_equiv_posns(cps, ps); % cts parent index +CPD.dps = find_equiv_posns(dps, ps); +if length(CPD.dps) ~= 1 + error('gmux must have exactly 1 discrete parent') +end +ss = fam_sz(end); +cpsz = fam_sz(CPD.cps(1)); % in gaussian_CPD, cpsz = sum(fam_sz(CPD.cps)) +if ~all(fam_sz(CPD.cps) == cpsz) + error('all cts parents must have same size') +end +dpsz = fam_sz(CPD.dps); +if dpsz ~= length(cps) + error(['the arity of the mux node is ' num2str(dpsz) ... + ' but there are ' num2str(length(cps)) ' cts parents']); +end + +% set default params +%CPD.mean = zeros(ss, 1); +%CPD.cov = eye(ss); +%CPD.weights = randn(ss, cpsz); +CPD.mean = zeros(ss, dpsz); +CPD.cov = 1*repmat(eye(ss), [1 1 dpsz]); +CPD.weights = randn(ss, cpsz, dpsz); + +args = varargin; +nargs = length(args); +for i=1:2:nargs + switch args{i}, + case 'mean', CPD.mean = args{i+1}; + case 'cov', CPD.cov = args{i+1}; + case 'weights', CPD.weights = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end +end + +%%%%%%%%%%% + +function CPD = init_fields() +% This ensures we define the fields in the same order +% no matter whether we load an object from a file, +% or create it from scratch. (Matlab requires this.) + +CPD.self = []; +CPD.sizes = []; +CPD.cps = []; +CPD.dps = []; +CPD.mean = []; +CPD.cov = []; +CPD.weights = []; + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@gmux_CPD/sample_node.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@gmux_CPD/sample_node.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,10 @@ +function y = sample_node(CPD, pev) +% SAMPLE_NODE Draw a random sample from P(Xi | x(pi_i), theta_i) (gmux) +% y = sample_node(CPD, parent_evidence) +% +% parent_ev{i} is the value of the i'th parent + +dpval = pev{CPD.dps}; +x = pev{CPD.cps(dpval)}; +y = gsamp(CPD.mean(:,dpval) + CPD.weights(:,:,dpval)*x(:), CPD.cov(:,:,dpval), 1); +y = y(:); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmm2Q_CPD/CPD_to_CPT.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmm2Q_CPD/CPD_to_CPT.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,35 @@ +function CPT = CPD_to_CPT(CPD) +% Compute the big CPT for an HHMM Q node (including F parents) +% by combining internal transprob and startprob +% function CPT = CPD_to_CPT(CPD) + +Qsz = CPD.Qsz; + +if ~isempty(CPD.Fbelow_ndx) + if ~isempty(CPD.Fself_ndx) % general case + error('not implemented') + else % no F from self, hence no startprob (top level) + nps = length(CPD.dom_sz)-1; % num parents + CPT = 0*myones(CPD.dom_sz); + % when Fself=1, the CPT(i,j) = delta(i,j) for all k + for k=1:prod(CPD.Qpsizes) + Qps_vals = ind2subv(CPD.Qpsizes, k); + ndx = mk_multi_index(nps+1, [CPD.Fbelow_ndx CPD.Qps_ndx], [1 Qps_vals]); + CPT(ndx{:}) = eye(Qsz); % CPT(:,2,k,:) or CPT(:,k,2,:) etc + end + ndx = mk_multi_index(nps+1, CPD.Fbelow_ndx, 2); + CPT(ndx{:}) = CPD.transprob; % we assume transprob is in topo order + end +else % no F signal from below + if ~isempty(CPD.Fself_ndx) % bottom level + nps = length(CPD.dom_sz)-1; % num parents + CPT = 0*myones(CPD.dom_sz); + ndx = mk_multi_index(nps+1, CPD.Fself_ndx, 1); + CPT(ndx{:}) = CPD.transprob; + ndx = mk_multi_index(nps+1, CPD.Fself_ndx, 2); + CPT(ndx{:}) = CPD.startprob; + else % no F from self + error('An hhmmQ node without any F parents is just a tabular_CPD') + end +end + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmm2Q_CPD/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmm2Q_CPD/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,6 @@ +/CPD_to_CPT.m/1.1.1.1/Tue Sep 24 12:46:46 2002// +/hhmm2Q_CPD.m/1.1.1.1/Tue Sep 24 22:34:40 2002// +/maximize_params.m/1.1.1.1/Tue Sep 24 22:44:36 2002// +/reset_ess.m/1.1.1.1/Tue Sep 24 22:36:16 2002// +/update_ess.m/1.1.1.1/Tue Sep 24 22:43:30 2002// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmm2Q_CPD/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmm2Q_CPD/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@hhmm2Q_CPD diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmm2Q_CPD/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmm2Q_CPD/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmm2Q_CPD/hhmm2Q_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmm2Q_CPD/hhmm2Q_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,65 @@ +function CPD = hhmm2Q_CPD(bnet, self, varargin) +% HHMMQ_CPD Make the CPD for a Q node in a 2 level hierarchical HMM +% CPD = hhmmQ_CPD(bnet, self, ...) +% +% Fself(t-1) Qps +% \ | +% \ v +% Qold(t-1) -> Q(t) +% / +% / +% Fbelow(t-1) +% +% +% optional args [defaults] +% +% Fself - node number <= ss +% Fbelow - node number <= ss +% Qps - node numbers (all <= 2*ss) - uses 2TBN indexing +% transprob - CPT for when Fbelow=2 and Fself=1 +% startprob - CPT for when Fbelow=2 and Fself=2 +% If Fbelow=1, we cannot change state. + +ss = bnet.nnodes_per_slice; +ns = bnet.node_sizes(:); + +% set default arguments +Fself = []; +Fbelow = []; +Qps = []; +startprob = []; +transprob = []; + +for i=1:2:length(varargin) + switch varargin{i}, + case 'Fself', Fself = varargin{i+1}; + case 'Fbelow', Fbelow = varargin{i+1}; + case 'Qps', Qps = varargin{i+1}; + case 'transprob', transprob = varargin{i+1}; + case 'startprob', startprob = varargin{i+1}; + end +end + +ps = parents(bnet.dag, self); +old_self = self-ss; +ndsz = ns(:)'; +CPD.dom_sz = [ndsz(ps) ns(self)]; +CPD.Fself_ndx = find_equiv_posns(Fself, ps); +CPD.Fbelow_ndx = find_equiv_posns(Fbelow, ps); +Qps = mysetdiff(ps, [Fself Fbelow old_self]); +CPD.Qps_ndx = find_equiv_posns(Qps, ps); +CPD.old_self_ndx = find_equiv_posns(old_self, ps); + +Qps = ps(CPD.Qps_ndx); +CPD.Qsz = ns(self); +CPD.Qpsizes = ns(Qps); + +CPD.transprob = transprob; +CPD.startprob = startprob; +CPD.start_counts = []; +CPD.trans_counts = []; + +CPD = class(CPD, 'hhmm2Q_CPD', discrete_CPD(0, CPD.dom_sz)); + + + diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmm2Q_CPD/maximize_params.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmm2Q_CPD/maximize_params.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,10 @@ +function CPD = maximize_params(CPD, temp) +% MAXIMIZE_PARAMS Set the params of a hhmmQ node to their ML/MAP values. +% CPD = maximize_params(CPD, temperature) + +if sum(CPD.start_counts(:)) > 0 + CPD.startprob = mk_stochastic(CPD.start_counts); +end +if sum(CPD.trans_counts(:)) > 0 + CPD.transprob = mk_stochastic(CPD.trans_counts); +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmm2Q_CPD/reset_ess.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmm2Q_CPD/reset_ess.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,12 @@ +function CPD = reset_ess(CPD) +% RESET_ESS Reset the Expected Sufficient Statistics of a hhmm2 Q node. +% CPD = reset_ess(CPD) + +domsz = CPD.dom_sz; +domsz(CPD.Fself_ndx) = 1; +domsz(CPD.Fbelow_ndx) = 1; +Qdom_sz = domsz; +Qdom_sz(Qdom_sz==1)=[]; % get rid of dimensions of size 1 + +CPD.start_counts = zeros(Qdom_sz); +CPD.trans_counts = zeros(Qdom_sz); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmm2Q_CPD/update_ess.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmm2Q_CPD/update_ess.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,26 @@ +function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv) + +marg = add_ev_to_dmarginal(fmarginal, evidence, ns); + +nps = length(CPD.dom_sz)-1; % num parents + +if ~isempty(CPD.Fbelow_ndx) + if ~isempty(CPD.Fself_ndx) % general case + ndx = mk_multi_index(nps+1, [CPD.Fbelow_ndx CPD.Fself_ndx], [2 1]); + CPD.trans_counts = CPD.trans_counts + squeeze(marg.T(ndx{:})); + ndx = mk_multi_index(nps+1, [CPD.Fbelow_ndx CPD.Fself_ndx], [2 2]); + CPD.start_counts = CPD.start_counts + squeeze(marg.T(ndx{:})); + else % no F from self, hence no startprob (top level) + ndx = mk_multi_index(nps+1, CPD.Fbelow_ndx, 2); + CPD.trans_counts = CPD.trans_counts + squeeze(marg.T(ndx{:})); + end +else % no F signal from below + if ~isempty(CPD.Fself_ndx) % self F (bottom level) + ndx = mk_multi_index(nps+1, CPD.Fself_ndx, 1); + CPD.trans_counts = CPD.trans_counts + squeeze(marg.T(ndx{:})); + ndx = mk_multi_index(nps+1, CPD.Fself_ndx, 2); + CPD.start_counts = CPD.start_counts + squeeze(marg.T(ndx{:})); + else % no F from self or below + error('no F signal') + end +end diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmmF_CPD/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmmF_CPD/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,7 @@ +/hhmmF_CPD.m/1.1.1.1/Mon Jun 24 23:38:24 2002// +/log_prior.m/1.1.1.1/Wed May 29 15:59:54 2002// +/maximize_params.m/1.1.1.1/Wed May 29 15:59:54 2002// +/reset_ess.m/1.1.1.1/Wed May 29 15:59:54 2002// +/update_CPT.m/1.1.1.1/Mon Jun 24 22:45:04 2002// +/update_ess.m/1.1.1.1/Mon Jun 24 23:54:30 2002// +D/Old//// diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmmF_CPD/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmmF_CPD/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@hhmmF_CPD diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmmF_CPD/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmmF_CPD/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmmF_CPD/Old/CVS/Entries --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmmF_CPD/Old/CVS/Entries Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,7 @@ +/hhmmF_CPD.m/1.1.1.1/Mon Jun 24 22:35:06 2002// +/log_prior.m/1.1.1.1/Mon Jun 24 22:35:06 2002// +/maximize_params.m/1.1.1.1/Mon Jun 24 22:35:06 2002// +/reset_ess.m/1.1.1.1/Mon Jun 24 22:35:06 2002// +/update_CPT.m/1.1.1.1/Mon Jun 24 22:35:06 2002// +/update_ess.m/1.1.1.1/Mon Jun 24 22:35:06 2002// +D diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmmF_CPD/Old/CVS/Repository --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmmF_CPD/Old/CVS/Repository Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +FullBNT/BNT/CPDs/@hhmmF_CPD/Old diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmmF_CPD/Old/CVS/Root --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmmF_CPD/Old/CVS/Root Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmmF_CPD/Old/hhmmF_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmmF_CPD/Old/hhmmF_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,76 @@ +function CPD = hhmmF_CPD(bnet, self, Qnodes, d, D, varargin) +% HHMMF_CPD Make the CPD for an F node at depth D of a D-level hierarchical HMM +% CPD = hhmmF_CPD(bnet, self, Qnodes, d, D, ...) +% +% Q(d-1) +% \ +% \ +% F(d) +% / | +% / | +% Q(d) F(d+1) +% +% We assume nodes are ordered (numbered) as follows: +% Q(1), ... Q(d), F(d+1), F(d) +% +% F(d)=2 means level d has finished. The prob this happens depends on Q(d) +% and optionally on Q(d-1), Q(d=1), ..., Q(1). +% Also, level d can only finish if the level below has finished +% (hence the F(d+1) -> F(d) arc). +% +% If d=D, there is no F(d+1), so F(d) is just a regular tabular_CPD. +% If all models always finish in the same state (e.g., their last), +% we don't need to condition on the state of parent models (Q(d-1), ...) +% +% optional args [defaults] +% +% termprob - termprob(k,i,2) = prob finishing given Q(d)=i and Q(1:d-1)=k [ finish in last state ] +% +% hhmmF_CPD is a subclass of tabular_CPD so we inherit inference methods like CPD_to_pot, etc. +% +% We create an isolated tabular_CPD with no F parent to learn termprob +% so we can avail of e.g., entropic or Dirichlet priors. +% +% For details, see "Linear-time inference in hierarchical HMMs", Murphy and Paskin, NIPS'01. + + +ps = parents(bnet.dag, self); +Qps = myintersect(ps, Qnodes); +F = mysetdiff(ps, Qps); +CPD.Q = Qps(end); % Q(d) +assert(CPD.Q == Qnodes(d)); +CPD.Qps = Qps(1:end-1); % all Q parents except Q(d), i.e., calling context + +ns = bnet.node_sizes(:); +CPD.Qsizes = ns(Qnodes); +CPD.d = d; +CPD.D = D; + +Qsz = ns(CPD.Q); +Qpsz = prod(ns(CPD.Qps)); + +% set default arguments +p = 0.9; +%termprob(k,i,t) Might terminate if i=Qsz; will not terminate if i % We sum over the possibilities that F(d+1) = 1 or 2 + +obs_self = ~hidden_bitv(Q); +if obs_self + self_val = evidence{Q}; +end + +if isempty(Qps) % independent of parent context + counts = zeros(Qsz, 2); + %fmarginal.T(Q(d), F(d+1), F(d)) + if obs_self + marg = myreshape(fmarginal.T, [1 2 2]); + counts(self_val,:) = marg(1,2,:); + %counts(self_val,:) = marg(1,1,:) + marg(1,2,:); + else + marg = myreshape(fmarginal.T, [Qsz 2 2]); + counts = squeeze(marg(:,2,:)); + %counts = squeeze(marg(:,2,:)) + squeeze(marg(:,1,:)); + end +else + counts = zeros(Qpsz, Qsz, 2); + %fmarginal.T(Q(1:d-1), Q(d), F(d+1), F(d)) + obs_Qps = ~any(hidden_bitv(Qps)); % we assume that all or none of the Q parents are observed + if obs_Qps + Qps_val = subv2ind(Qpsz, cat(1, evidence{Qps})); + end + if obs_self & obs_Qps + marg = myreshape(fmarginal.T, [1 1 2 2]); + counts(Qps_val, self_val, :) = squeeze(marg(1,1,2,:)); + %counts(Qps_val, self_val, :) = squeeze(marg(1,1,2,:)) + squeeze(marg(1,1,1,:)); + elseif ~obs_self & obs_Qps + marg = myreshape(fmarginal.T, [1 Qsz 2 2]); + counts(Qps_val, :, :) = squeeze(marg(1,:,2,:)); + %counts(Qps_val, :, :) = squeeze(marg(1,:,2,:)) + squeeze(marg(1,:,1,:)); + elseif obs_self & ~obs_Qps + error('not yet implemented') + else + marg = myreshape(fmarginal.T, [Qpsz Qsz 2 2]); + counts(:, :, :) = squeeze(marg(:,:,2,:)); + %counts(:, :, :) = squeeze(marg(:,:,2,:)) + squeeze(marg(:,:,1,:)); + end +end + +CPD.sub_CPD_term = update_ess_simple(CPD.sub_CPD_term, counts); diff -r 12abff5474c8 -r b5b38998ef3b _FullBNT/BNT/CPDs/@hhmmF_CPD/hhmmF_CPD.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/_FullBNT/BNT/CPDs/@hhmmF_CPD/hhmmF_CPD.m Fri Apr 11 15:54:25 2014 +0100 @@ -0,0 +1,73 @@ +function CPD = hhmmF_CPD(bnet, self, Qself, Fbelow, varargin) +% HHMMF_CPD Make the CPD for an F node in a hierarchical HMM +% CPD = hhmmF_CPD(bnet, self, Qself, Fbelow, ...) +% +% Qps +% \ +% \ +% Fself +% / | +% / | +% Qself Fbelow +% +% We assume nodes are ordered (numbered) as follows: Qps, Q, Fbelow, F +% All nodes numbers should be from slice 1. +% +% If Fbelow if missing, this becomes a regular tabular_CPD. +% Qps may be omitted. +% +% optional args [defaults] +% +% Qps - node numbers. +% termprob - termprob(k,i,2) = prob finishing given Q(d)=i and Q(1:d-1)=k [ finish in last state wp 0.9] +% +% hhmmF_CPD is a subclass of tabular_CPD so we inherit inference methods like CPD_to_pot, etc. +% +% We create an isolated tabular_CPD with no F parent to learn termprob +% so we can avail of e.g., entropic or Dirichlet priors. +% +% For details, see "Linear-time inference in hierarchical HMMs", Murphy and Paskin, NIPS'01. + + + +Qps = []; +% get parents +for i=1:2:length(varargin) + switch varargin{i}, + case 'Qps', Qps = varargin{i+1}; + end +end + +ns = bnet.node_sizes(:); +Qsz = ns(Qself); +Qpsz = prod(ns(Qps)); +CPD.Qsz = Qsz; +CPD.Qpsz = Qpsz; + +ps = parents(bnet.dag, self); +CPD.Fbelow_ndx = find_equiv_posns(Fbelow, ps); +CPD.Qps_ndx = find_equiv_posns(Qps, ps); +CPD.Qself_ndx = find_equiv_posns(Qself, ps); + +% set default arguments +p = 0.9; +%termprob(k,i,t) Might terminate if i=Qsz; will not terminate if i