diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD')
18 files changed, 560 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Entries new file mode 100644 index 00000000..1e0984dc --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Entries @@ -0,0 +1,11 @@ +/convert_to_pot.m/1.1.1.1/Wed May 29 15:59:54 2002// +/convert_to_table.m/1.1.1.1/Tue Mar 30 17:19:22 2004// +/display.m/1.1.1.1/Wed May 29 15:59:54 2002// +/get_field.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// +/sample_node.m/1.1.1.1/Wed May 29 15:59:54 2002// +/set_fields.m/1.1.1.1/Wed May 29 15:59:54 2002// +/softmax_CPD.m/1.1.1.1/Tue Jan 7 16:25:14 2003// +/update_ess.m/1.1.1.1/Wed May 29 15:59:54 2002// +D diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Entries.Log new file mode 100644 index 00000000..b2cd71e0 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Entries.Log @@ -0,0 +1 @@ +A D/private//// diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Repository new file mode 100644 index 00000000..d5dac28b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/CPDs/@softmax_CPD diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/convert_to_pot.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/convert_to_pot.m new file mode 100644 index 00000000..518f4a50 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/convert_to_pot.m @@ -0,0 +1,58 @@ +function pot = convert_to_pot(CPD, pot_type, domain, evidence) +% CONVERT_TO_POT Convert a softmax CPD to a potential +% pots = 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; + +odom = domain(~isemptycell(evidence(domain))); +T = convert_to_table(CPD, domain, evidence); + +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 + ps = domain(1:end-1); + dps = ps(CPD.dpndx); + cps = ps(CPD.cpndx); + ddom = [dps CPD.self]; + cdom = cps; + pot = cgpot(ddom, cdom, 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 --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/convert_to_table.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/convert_to_table.m new file mode 100644 index 00000000..f703d79b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/convert_to_table.m @@ -0,0 +1,52 @@ +function T = convert_to_table(CPD, domain, evidence) +% CONVERT_TO_TABLE Convert a softmax CPD to a table, incorporating any evidence +% T = convert_to_table(CPD, domain, evidence) + +self = domain(end); +ps = domain(1:end-1); +cnodes = domain(CPD.cpndx); +cps = myintersect(ps, cnodes); +dps = domain(CPD.dpndx); +dps_as_cps = domain(CPD.dps_as_cps.ndx); +all_dps = union(dps,dps_as_cps); +odom = domain(~isemptycell(evidence(domain))); +if ~isempty(cps), assert(myismember(cps, odom)); end % all cts parents must be observed + +ns = zeros(1, max(domain)); +ns(domain) = CPD.sizes; +ens = ns; % effective node sizes +ens(odom) = 1; + +% dpsize >= glimsz because the glm parameters are tied across the dps_as_cps parents +dpsize = prod(ens(all_dps)); % size of ALL self'discrete parents +dpvals = cat(1, evidence{myintersect(all_dps, odom)}); +cpvals = cat(1, evidence{cps}); +if ~isempty(dps_as_cps), + separator = CPD.dps_as_cps.separator; + dp_as_cpmap = find_equiv_posns(dps_as_cps, all_dps); + dops_map = find_equiv_posns(myintersect(all_dps, odom), all_dps); + puredp_map = find_equiv_posns(dps, all_dps); + subs = ind2subv(ens(all_dps), 1:prod(ens(all_dps))); + if ~isempty(dops_map), subs(:,dops_map) = subs(:,dops_map)+repmat(dpvals(:)',[size(subs,1) 1])-1; end +end + +[w,b] = extract_params(CPD); +T = zeros(dpsize, ns(self)); +for i=1:dpsize, + active_glm = i; + dp_as_cpvals=zeros(1,sum(ns(dps_as_cps))); + if ~isempty(dps_as_cps), + active_glm = max([1,subv2ind(ns(dps), subs(i,puredp_map))]); + % Extract the params compatible with the observations (if any) on the 'pure' discrete parents (if any) + where_one = separator + subs(i,dp_as_cpmap); + % and get in the dp_as_cp parents... + dp_as_cpvals(where_one)=1; + end + T(i,:) = normalise(exp([dp_as_cpvals(:); cpvals(:)]'*w(:,:,active_glm) + b(:,active_glm)')); +end +if myismember(self, odom) + r = evidence{self}; + T = T(:,r); +end + +T = myreshape(T, ens(domain)); diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/display.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/display.m new file mode 100644 index 00000000..06a0f02c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/display.m @@ -0,0 +1,4 @@ +function display(CPD) + +disp('softmax_CPD object'); +disp(struct(CPD)); diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/get_field.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/get_field.m new file mode 100644 index 00000000..240f1fd7 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/get_field.m @@ -0,0 +1,18 @@ +function val = get_params(CPD, name) +% GET_PARAMS Get the parameters (fields) for a softmax_CPD object +% val = get_params(CPD, name) +% +% The following fields can be accessed +% +% weights - W(X,Y,Q) +% offset - b(Y,Q) +% +% e.g., W = get_params(CPD, 'weights') + +[W, b] = extract_params(CPD); +switch name + case 'weights', val = W; + case 'offset', val = b; + otherwise, + error(['invalid argument name ' name]); +end diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/maximize_params.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/maximize_params.m new file mode 100644 index 00000000..15c94dd5 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/maximize_params.m @@ -0,0 +1,41 @@ +function CPD = maximize_params(CPD, temp) +% MAXIMIZE_PARAMS Set the params of a CPD to their ML values (dsoftmax) using IRLS +% CPD = maximize_params(CPD, temperature) +% temperature parameter is ignored + +% Written by Pierpaolo Brutti + +if ~adjustable_CPD(CPD), return; end +options = foptions; + +if CPD.verbose + options(1) = 1; +else + options(1) = -1; +end +%options(1) = CPD.verbose; + +options(2) = CPD.wthresh; +options(3) = CPD.llthresh; +options(5) = CPD.approx_hess; +options(14) = CPD.max_iter; + +dpsize = size(CPD.self_vals,3); +for i=1:dpsize, + mask=find(CPD.eso_weights(:,:,i)>0); % for adapting the parameters we use only positive weighted example + if ~isempty(mask), + if ~isempty(CPD.dps_as_cps.ndx), + puredp_map = find_equiv_posns(CPD.dpndx, union(CPD.dpndx, CPD.dps_as_cps.ndx)); % find the glm structure + subs = ind2subv(CPD.sizes(union(CPD.dpndx, CPD.dps_as_cps.ndx)),i); % that corrisponds to the + active_glm = max([1,subv2ind(CPD.sizes(CPD.dpndx), subs(puredp_map))]); % i-th 'fictitious' example + + CPD.glim{active_glm} = netopt_weighted(CPD.glim{active_glm}, options, CPD.parent_vals(mask',:,i),... + CPD.self_vals(mask',:,i), CPD.eso_weights(mask',:,i), 'scg'); + else + alfa = 0.4; if CPD.solo, alfa = 1; end % learning step = 1 <=> self is all alone in the net + CPD.glim{i} = glmtrain_weighted(CPD.glim{i}, options, CPD.parent_vals(mask',:),... + CPD.self_vals(mask',:,i), CPD.eso_weights(mask',:,i), alfa); + end + end + mask=[]; +end diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Entries new file mode 100644 index 00000000..b6610f0d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Entries @@ -0,0 +1,2 @@ +/extract_params.m/1.1.1.1/Wed May 29 15:59:54 2002// +D diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Repository new file mode 100644 index 00000000..1667449e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/CPDs/@softmax_CPD/private diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/extract_params.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/extract_params.m new file mode 100644 index 00000000..486af06e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/extract_params.m @@ -0,0 +1,18 @@ +function [W, b] = extract_params(CPD) + +% W(X,Y,Q), b(Y,Q) where Y = ns(self), X = ns(cps), Q = prod(ns(dps)) + +glimsz = prod(CPD.sizes(CPD.dpndx)); +ss = CPD.sizes(end); +cpsz = sum(CPD.sizes(CPD.cpndx)); +dp_as_cpsz = sum(CPD.sizes(CPD.dps_as_cps.ndx)); +W = zeros(dp_as_cpsz + cpsz, ss, glimsz); +b = zeros(ss, glimsz); + +for i=1:glimsz + W(:,:,i) = CPD.glim{i}.w1; + b(:,i) = CPD.glim{i}.b1(:); +end + +W = myreshape(W, [dp_as_cpsz + cpsz ss CPD.sizes(CPD.dpndx)]); +b = myreshape(b, [ss CPD.sizes(CPD.dpndx)]); diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/reset_ess.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/reset_ess.m new file mode 100644 index 00000000..abf7d54e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/reset_ess.m @@ -0,0 +1,8 @@ +function CPD = reset_ess(CPD) +% RESET_ESS Reset the Expected Sufficient Statistics for a CPD (dsoftmax) +% CPD = reset_ess(CPD) + +CPD.parent_vals = []; +CPD.eso_weights=[]; +CPD.self_vals = []; +CPD.nsamples = 0; diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/sample_node.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/sample_node.m new file mode 100644 index 00000000..1c519049 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/sample_node.m @@ -0,0 +1,14 @@ +function y = sample_node(CPD, pvals) +% SAMPLE_NODE Draw a random sample from P(Xi | x(pi_i), theta_i) (discrete) +% y = sample_node(CPD, parent_evidence) +% +% parent_evidence{i} is the value of the i'th parent + +n = length(pvals)+1; +dom = 1:n; +%evidence = cell(1,n); +%evidence(1:n-1) = pvals(:)'; +evidence = pvals; +evidence{end+1} = []; +T = convert_to_table(CPD, dom, evidence); +y = sample_discrete(T); diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/set_fields.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/set_fields.m new file mode 100644 index 00000000..6c64b197 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/set_fields.m @@ -0,0 +1,45 @@ +function CPD = set_params(CPD, varargin) +% SET_PARAMS Set the parameters (fields) for a softmax_CPD object +% CPD = set_params(CPD, name/value pairs) +% +% The following optional arguments can be specified in the form of name/value pairs: +% (Let ns(i) be the size of node i, X = ns(X), Y = ns(Y), Q1=ns(dps(1)), Q2=ns(dps(2)), ... +% where dps are the discrete parents; if there are no discrete parents, we set Q1=1.) +% +% weights - (W(:,j,a,b,...) - W(:,j',a,b,...)) is ppn to dec. boundary +% between j,j' given Q1=a,Q2=b,... [ randn(X,Y,Q1,Q2,...) ] +% offset - (offset(j,a,b,...) - offset(j',a,b,...)) is the offset to dec. boundary +% between j,j' given Q1=a,Q2=b,... [ randn(Y,Q1,Q2,...) ] +% clamped - 'yes' means don't adjust params during learning ['no'] +% max_iter - the maximum number of steps to take [10] +% verbose - 'yes' means print the LL at each step of IRLS ['no'] +% wthresh - convergence threshold for weights [1e-2] +% llthresh - convergence threshold for log likelihood [1e-2] +% approx_hess - 'yes' means approximate the Hessian for speed ['no'] +% +% e.g., CPD = set_params(CPD,'offset', zeros(ns(i),1)); + +args = varargin; +nargs = length(args); +glimsz = prod(CPD.sizes(CPD.dpndx)); +for i=1:2:nargs + switch args{i}, + case 'discrete', str='nothing to do'; + case 'clamped', CPD = set_clamped(CPD, strcmp(args{i+1}, 'yes')); + case 'max_iter', CPD.max_iter = args{i+1}; + case 'verbose', CPD.verbose = strcmp(args{i+1}, 'yes'); + case 'max_iter', CPD.max_iter = args{i+1}; + case 'wthresh', CPD.wthresh = args{i+1}; + case 'llthresh', CPD.llthresh = args{i+1}; + case 'approx_hess', CPD.approx_hess = strcmp(args{i+1}, 'yes'); + case 'weights', for q=1:glimsz, CPD.glim{q}.w1 = args{i+1}(:,:,q); end; + case 'offset', + if glimsz == 1 + CPD.glim{1}.b1 = args{i+1}; + else + for q=1:glimsz, CPD.glim{q}.b1 = args{i+1}(:,q); end; + end + otherwise, + error(['invalid argument name ' args{i}]); + end +end diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/softmax_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/softmax_CPD.m new file mode 100644 index 00000000..3d2e5153 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/softmax_CPD.m @@ -0,0 +1,187 @@ +function CPD = softmax_CPD(bnet, self, varargin) +% SOFTMAX_CPD Make a softmax (multinomial logit) CPD +% +% To define this CPD precisely, let W be an (m x n) matrix with W(i,:) = {i-th row of B} +% => we can define the following vectorial function: +% +% softmax: R^n |--> R^m +% softmax(z,i-th)=exp(W(i,:)*z)/sum_k(exp(W(k,:)*z)) +% +% (this constructor augments z with a one at the beginning to introduce an offset term (=bias, intercept)) +% Now call the continuous (cts) and always observed (obs) parents X, +% the discrete parents (if any) Q, and this node Y then we use the discrete parent(s) just to index +% the parameter vectors (c.f., conditional Gaussian nodes); that is: +% prob(Y=i | X=x, Q=j) = softmax(x,i-th|j) +% where '|j' means that we are using the j-th (m x n) parameters matrix W(:,:,j). +% If there are no discrete parents, this is a regular softmax node. +% If Y is binary, this is a logistic (sigmoid) function. +% +% CPD = softmax_CPD(bnet, node_num, ...) will create a softmax CPD with random parameters, +% where node is the number of a node in this equivalence class. +% +% The following optional arguments can be specified in the form of name/value pairs: +% [default value in brackets] +% (Let ns(i) be the size of node i, X = ns(X), Y = ns(Y), Q1=ns(dps(1)), Q2=ns(dps(2)), ... +% where dps are the discrete parents; if there are no discrete parents, we set Q1=1.) +% +% discrete - the discrete parents that we want to treat like the cts ones [ [] ]. +% This can be used to define sigmoid belief network - see below the reference. +% For example suppose that Y has one cts parents X and two discrete ones: Q, C1 where: +% -> Q is binary (1/2) and used just to index the parameters of 'self' +% -> C1 is ternary (1/2/3) and treated as a cts node <=> its values appear into the linear +% part of the softmax function +% then: +% prob(Y|X=x, Q=q, C1=c1)= softmax(W(:,:,q)' * y) +% where y = [1 | delta(C1,1) delta(C1,2) delta(C1,3) | x(:)']' and delta(Y,a)=indicator(Y=a). +% weights - (w(:,j,a,b,...) - w(:,j',a,b,...)) is ppn to dec. boundary +% between j,j' given Q1=a,Q2=b,... [ randn(X,Y,Q1,Q2,...) ] +% offset - (b(j,a,b,...) - b(j',a,b,...)) is the offset to dec. boundary +% between j,j' given Q1=a,Q2=b,... [ randn(Y,Q1,Q2,...) ] +% +% e.g., CPD = softmax_CPD(bnet, i, 'offset', zeros(ns(i),1)); +% +% The following fields control the behavior of the M step, which uses +% a weighted version of the Iteratively Reweighted Least Squares (WIRLS) if dps_as_cps=[]; or +% a weighted SCG otherwise, as implemented in Netlab, and modified by Pierpaolo Brutti. +% +% clamped - 'yes' means don't adjust params during learning ['no'] +% max_iter - the maximum number of steps to take [10] +% verbose - 'yes' means print the LL at each step of IRLS ['no'] +% wthresh - convergence threshold for weights [1e-2] +% llthresh - convergence threshold for log likelihood [1e-2] +% approx_hess - 'yes' means approximate the Hessian for speed ['no'] +% +% For backwards compatibility with BNT2, you can also specify the parameters in the following order +% softmax_CPD(bnet, self, w, b, clamped, max_iter, verbose, wthresh, llthresh, approx_hess) +% +% REFERENCE +% For details on the sigmoid belief nets, see: +% - Neal (1992). Connectionist learning of belief networks, Artificial Intelligence, 56, 71-113. +% - Saul, Jakkola, Jordan (1996). Mean field theory for sigmoid belief networks, Journal of Artificial Intelligence Reseach (4), pagg. 61-76. +% +% For details on the M step, see: +% - K. Chen, L. Xu, H. Chi (1999). Improved learning algorithms for mixtures of experts in multiclass +% classification. Neural Networks 12, pp. 1229-1252. +% - M.I. Jordan, R.A. Jacobs (1994). Hierarchical Mixtures of Experts and the EM algorithm. +% Neural Computation 6, pp. 181-214. +% - S.R. Waterhouse, A.J. Robinson (1994). Classification Using Hierarchical Mixtures of Experts. In Proc. IEEE +% Workshop on Neural Network for Signal Processing IV, pp. 177-186 + +if nargin==0 + % This occurs if we are trying to load an object from a file. + CPD = init_fields; + CPD = class(CPD, 'softmax_CPD', discrete_CPD(0, [])); + return; +elseif isa(bnet, 'softmax_CPD') + % This might occur if we are copying an object. + CPD = bnet; + return; +end +CPD = init_fields; + +assert(myismember(self, bnet.dnodes)); +ns = bnet.node_sizes; +ps = parents(bnet.dag, self); +dps = myintersect(ps, bnet.dnodes); +cps = myintersect(ps, bnet.cnodes); + +clamped = 0; +CPD = class(CPD, 'softmax_CPD', discrete_CPD(clamped, ns([ps self]))); + +dps_as_cpssz = 0; +dps_as_cps = []; +% determine if any discrete parents are to be treated as cts +if nargin >= 3 && isstr(varargin{1}) % might have passed in 'discrete' + for i=1:2:length(varargin) + if strcmp(varargin{i}, 'discrete') + dps_as_cps = varargin{i+1}; + assert(myismember(dps_as_cps, dps)); + dps = mysetdiff(dps, dps_as_cps); % put out the dps treated as cts + CPD.dps_as_cps.ndx = find_equiv_posns(dps_as_cps, ps); + CPD.dps_as_cps.separator = [0 cumsum(ns(dps_as_cps(1:end-1)))]; % concatenated dps_as_cps dims separators + dps_as_cpssz = sum(ns(dps_as_cps)); + break; + end + end +end +assert(~isempty(union(cps, dps_as_cps))); % It have to be at least a cts or a dps_as_cps parents +self_size = ns(self); +cpsz = sum(ns(cps)); +glimsz = prod(ns(dps)); +CPD.dpndx = find_equiv_posns(dps, ps); % it contains only the indeces of the 'pure' dps +CPD.cpndx = find_equiv_posns(cps, ps); + +CPD.self = self; +CPD.solo = (length(ns)<=2); +CPD.sizes = bnet.node_sizes([ps self]); + +% set default params +CPD.max_iter = 10; +CPD.verbose = 0; +CPD.wthresh = 1e-2; +CPD.llthresh = 1e-2; +CPD.approx_hess = 0; +CPD.glim = cell(1,glimsz); +for i=1:glimsz + CPD.glim{i} = glm(dps_as_cpssz + cpsz, self_size, 'softmax'); +end + +if nargin >= 3 + args = varargin; + nargs = length(args); + if ~isstr(args{1}) + % softmax_CPD(bnet, self, w, b, clamped, max_iter, verbose, wthresh, llthresh, approx_hess) + if nargs >= 1 && ~isempty(args{1}), CPD = set_fields(CPD, 'weights', args{1}); end + if nargs >= 2 && ~isempty(args{2}), CPD = set_fields(CPD, 'offset', args{2}); end + if nargs >= 3 && ~isempty(args{3}), CPD = set_clamped(CPD, args{3}); end + if nargs >= 4 && ~isempty(args{4}), CPD.max_iter = args{4}; end + if nargs >= 5 && ~isempty(args{5}), CPD.verbose = args{5}; end + if nargs >= 6 && ~isempty(args{6}), CPD.wthresh = args{6}; end + if nargs >= 7 && ~isempty(args{7}), CPD.llthresh = args{7}; end + if nargs >= 8 && ~isempty(args{8}), CPD.approx_hess = args{8}; end + else + CPD = set_fields(CPD, args{:}); + end +end + +% sufficient statistics +% Since dsoftmax is not in the exponential family, we must store all the raw data. +CPD.parent_vals = []; % X(l,:) = value of cts parents in l'th example +CPD.self_vals = []; % Y(l,:) = value of self in l'th example + +CPD.eso_weights=[]; % weights used by the WIRLS algorithm + +% For BIC +CPD.nsamples = 0; +if ~adjustable_CPD(CPD), + CPD.nparams=0; +else + [W, b] = extract_params(CPD); + CPD.nparams= prod(size(W)) + prod(size(b)); +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.glim = {}; +CPD.self = []; +CPD.solo = []; +CPD.max_iter = []; +CPD.verbose = []; +CPD.wthresh = []; +CPD.llthresh = []; +CPD.approx_hess = []; +CPD.sizes = []; +CPD.parent_vals = []; +CPD.eso_weights=[]; +CPD.self_vals = []; +CPD.nsamples = []; +CPD.nparams = []; +CPD.dpndx = []; +CPD.cpndx = []; +CPD.dps_as_cps.ndx = []; +CPD.dps_as_cps.separator = []; diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/update_ess.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/update_ess.m new file mode 100644 index 00000000..143c567c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/update_ess.m @@ -0,0 +1,97 @@ +function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv) +% UPDATE_ESS Update the Expected Sufficient Statistics of a softmax node +% function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv) +% +% fmarginal = overall posterior distribution of self and its parents +% fmarginal(i1,i2...,ik,s)=prob(Pa1=i1,...,Pak=ik, self=s| X) +% +% => 1) prob(self|Pa1,...,Pak)=fmarginal/prob(Pa1,...,Pak) with prob(Pa1,...,Pak)=sum{s,fmarginal} +% [self estimation -> CPD.self_vals] +% 2) prob(Pa1,...,Pak) [WIRLS weights -> CPD.eso_weights] +% +% Hidden_bitv is ignored + +% Written by Pierpaolo Brutti + +if ~adjustable_CPD(CPD), return; end + +domain = fmarginal.domain; +self = domain(end); +ps = domain(1:end-1); +cnodes = domain(CPD.cpndx); +cps = myintersect(domain, cnodes); +dps = mysetdiff(ps, cps); +dn_use = dps; +if isempty(evidence{self}) dn_use = [dn_use self]; end % if self is hidden we must consider its dimension +dps_as_cps = domain(CPD.dps_as_cps.ndx); +odom = domain(~isemptycell(evidence(domain))); + +ns = zeros(1, max(domain)); +ns(domain) = CPD.sizes; % CPD.sizes = bnet.node_sizes([ps self]); +ens = ns; % effective node sizes +ens(odom) = 1; +dpsize = prod(ns(dps)); + +% Extract the params compatible with the observations (if any) on the discrete parents (if any) +dops = myintersect(dps, odom); +dpvals = cat(1, evidence{dops}); + +subs = ind2subv(ens(dn_use), 1:prod(ens(dn_use))); +dpmap = find_equiv_posns(dops, dn_use); +if ~isempty(dpmap), subs(:,dpmap) = subs(:,dpmap)+repmat(dpvals(:)',[size(subs,1) 1])-1; end +supportedQs = subv2ind(ns(dn_use), subs); subs=subs(1:prod(ens(dps)),1:length(dps)); +Qarity = prod(ns(dn_use)); +if isempty(dn_use), Qarity = 1; end + +fullm.T = zeros(Qarity, 1); +fullm.T(supportedQs) = fmarginal.T(:); +rs_dim = CPD.sizes; rs_dim(CPD.cpndx) = 1; % +if ~isempty(evidence{self}), rs_dim(end)=1; end % reshaping the marginal +fullm.T = reshape(fullm.T, rs_dim); % + +% --------------------------------------------------------------------------------UPDATE-- + +CPD.nsamples = CPD.nsamples + 1; + +% 1) observations vector -> CPD.parents_vals --------------------------------------------- +cpvals = cat(1, evidence{cps}); + +if ~isempty(dps_as_cps), % ...get in the dp_as_cp parents... + separator = CPD.dps_as_cps.separator; + dp_as_cpmap = find_equiv_posns(dps_as_cps, dps); + for i=1:dpsize, + dp_as_cpvals=zeros(1,sum(ns(dps_as_cps))); + possible_vals = ind2subv(ns(dps),i); + ll=find(ismember(subs(:,dp_as_cpmap), possible_vals(dp_as_cpmap), 'rows')==1); + if ~isempty(ll), + where_one = separator + possible_vals(dp_as_cpmap); + dp_as_cpvals(where_one)=1; + end + CPD.parent_vals(CPD.nsamples,:,i) = [dp_as_cpvals(:); cpvals(:)]'; + end +else + CPD.parent_vals(CPD.nsamples,:) = cpvals(:)'; +end + +% 2) weights vector -> CPD.eso_weights ---------------------------------------------------- +if isempty(evidence{self}), % self is hidden + pesi=reshape(sum(fullm.T, length(rs_dim)),[dpsize,1]); +else + pesi=reshape(fullm.T,[dpsize,1]); +end +assert(approxeq(sum(pesi),1)); % check + +% 3) estimate (if R is hidden) or recover (if R is obs) self'value------------------------- +if isempty(evidence{self}) % P(self|Pa1,...,Pak)=fmarginal/prob(Pa1,...,Pak) + r=reshape(mk_stochastic(fullm.T), [dpsize ns(self)]); % matrix size: prod{j,ns(Paj)} x ns(self) +else + r = zeros(dpsize,ns(self)); + for i=1:dpsize, if pesi(i)~=0, r(i,evidence{self}) = 1; end; end +end +for i=1:dpsize, if pesi(i)~=0, assert(approxeq(sum(r(i,:)),1)); end; end % check + +% 4) save the previous values -------------------------------------------------------------- +for i=1:dpsize + CPD.eso_weights(CPD.nsamples,:,i)=pesi(i); + CPD.self_vals(CPD.nsamples,:,i) = r(i,:); +end |