1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
|
function CPD = mlp_CPD(bnet, self, nhidden, w1, b1, w2, b2, clamped, max_iter, verbose, wthresh, llthresh)
% MLP_CPD Make a CPD from a Multi Layer Perceptron (i.e., feedforward neural network)
%
% We use a different MLP for each discrete parent combination (if there are any discrete parents).
% We currently assume this node (the child) is discrete.
%
% CPD = mlp_CPD(bnet, self, nhidden)
% will create a CPD with random parameters, where self is the number of this node and nhidden the number of the hidden nodes.
% The params are drawn from N(0, s*I), where s = 1/sqrt(n+1), n = length(X).
%
% CPD = mlp_CPD(bnet, self, nhidden, w1, b1, w2, b2) allows you to specify the params, where
% w1 = first-layer weight matrix
% b1 = first-layer bias vector
% w2 = second-layer weight matrix
% b2 = second-layer bias vector
% These are assumed to be the same for each discrete parent combination.
% If any of these are [], random values will be created.
%
% CPD = mlp_CPD(bnet, self, nhidden, w1, b1, w2, b2, clamped) allows you to prevent the params from being
% updated during learning (if clamped = 1). Default: clamped = 0.
%
% CPD = mlp_CPD(bnet, self, nhidden, w1, b1, w2, b2, clamped, max_iter, verbose, wthresh, llthresh)
% alllows you to specify params that control the M step:
% max_iter - the maximum number of steps to take (default: 10)
% verbose - controls whether to print (default: 0 means silent).
% wthresh - a measure of the precision required for the value of
% the weights W at the solution. Default: 1e-2.
% llthresh - a measure of the precision required of the objective
% function (log-likelihood) at the solution. Both this and the previous condition must
% be satisfied for termination. Default: 1e-2.
%
% For learning, we use a weighted version of scaled conjugated gradient in the M step.
if nargin==0
% This occurs if we are trying to load an object from a file.
CPD = init_fields;
CPD = class(CPD, 'mlp_CPD', discrete_CPD(0,[]));
return;
elseif isa(bnet, 'mlp_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);
dnodes = mysetdiff(1:length(bnet.dag), bnet.cnodes);
dps = myintersect(ps, dnodes);
cps = myintersect(ps, bnet.cnodes);
dpsz = prod(ns(dps));
cpsz = sum(ns(cps));
self_size = ns(self);
% discrete/cts parent index - which ones of my parents are discrete/cts?
CPD.dpndx = find_equiv_posns(dps, ps);
CPD.cpndx = find_equiv_posns(cps, ps);
CPD.mlp = cell(1,dpsz);
for i=1:dpsz
CPD.mlp{i} = mlp(cpsz, nhidden, self_size, 'softmax');
if nargin >=4 & ~isempty(w1)
CPD.mlp{i}.w1 = w1;
end
if nargin >=5 & ~isempty(b1)
CPD.mlp{i}.b1 = b1;
end
if nargin >=6 & ~isempty(w2)
CPD.mlp{i}.w2 = w2;
end
if nargin >=7 & ~isempty(b2)
CPD.mlp{i}.b2 = b2;
end
W1app(:,:,i)=CPD.mlp{i}.w1;
W2app(:,:,i)=CPD.mlp{i}.w2;
b1app(i,:)=CPD.mlp{i}.b1;
b2app(i,:)=CPD.mlp{i}.b2;
end
if nargin < 8, clamped = 0; end
if nargin < 9, max_iter = 10; end
if nargin < 10, verbose = 0; end
if nargin < 11, wthresh = 1e-2; end
if nargin < 12, llthresh = 1e-2; end
CPD.self = self;
CPD.max_iter = max_iter;
CPD.verbose = verbose;
CPD.wthresh = wthresh;
CPD.llthresh = llthresh;
% sufficient statistics
% Since MLP is not in the exponential family, we must store all the raw data.
%
CPD.W1=W1app; % Extract all the parameters of the node for handling discrete obs parents
CPD.W2=W2app; %
nparaW=[size(W1app) size(W2app)]; %
CPD.b1=b1app; %
CPD.b2=b2app; %
nparab=[size(b1app) size(b2app)]; %
CPD.sizes=bnet.node_sizes(:); % used in CPD_to_table to pump up the node sizes
CPD.parent_vals = []; % X(l,:) = value of cts parents in l'th example
CPD.eso_weights=[]; % weights used by the SCG algorithm
CPD.self_vals = []; % Y(l,:) = value of self in l'th example
% For BIC
CPD.nsamples = 0;
CPD.nparams=prod(nparaW)+prod(nparab);
CPD = class(CPD, 'mlp_CPD', discrete_CPD(clamped, ns([ps self])));
%%%%%%%%%%%
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.mlp = {};
CPD.self = [];
CPD.max_iter = [];
CPD.verbose = [];
CPD.wthresh = [];
CPD.llthresh = [];
CPD.approx_hess = [];
CPD.W1 = [];
CPD.W2 = [];
CPD.b1 = [];
CPD.b2 = [];
CPD.sizes = [];
CPD.parent_vals = [];
CPD.eso_weights=[];
CPD.self_vals = [];
CPD.nsamples = [];
CPD.nparams = [];
|