wolffd@0: function CPD = gaussian_CPD(bnet, self, varargin)
wolffd@0: % GAUSSIAN_CPD Make a conditional linear Gaussian distrib.
wolffd@0: %
wolffd@0: % CPD = gaussian_CPD(bnet, node, ...) will create a CPD with random parameters,
wolffd@0: % where node is the number of a node in this equivalence class.
wolffd@0: 
wolffd@0: % To define this CPD precisely, call the continuous (cts) parents (if any) X,
wolffd@0: % the discrete parents (if any) Q, and this node Y. Then the distribution on Y is:
wolffd@0: % - no parents: Y ~ N(mu, Sigma)
wolffd@0: % - cts parents : Y|X=x ~ N(mu + W x, Sigma)
wolffd@0: % - discrete parents: Y|Q=i ~ N(mu(i), Sigma(i))
wolffd@0: % - cts and discrete parents: Y|X=x,Q=i ~ N(mu(i) + W(i) x, Sigma(i))
wolffd@0: %
wolffd@0: % The list below gives optional arguments [default value in brackets].
wolffd@0: % (Let ns(i) be the size of node i, X = ns(X), Y = ns(Y) and Q = prod(ns(Q)).)
wolffd@0: % Parameters will be reshaped to the right size if necessary.
wolffd@0: %
wolffd@0: % mean       - mu(:,i) is the mean given Q=i [ randn(Y,Q) ]
wolffd@0: % cov        - Sigma(:,:,i) is the covariance given Q=i [ repmat(100*eye(Y,Y), [1 1 Q]) ]
wolffd@0: % weights    - W(:,:,i) is the regression matrix given Q=i [ randn(Y,X,Q) ]
wolffd@0: % cov_type   - if 'diag', Sigma(:,:,i) is diagonal [ 'full' ]
wolffd@0: % tied_cov   - if 1, we constrain Sigma(:,:,i) to be the same for all i [0]
wolffd@0: % clamp_mean - if 1, we do not adjust mu(:,i) during learning [0]
wolffd@0: % clamp_cov  - if 1, we do not adjust Sigma(:,:,i) during learning [0]
wolffd@0: % clamp_weights - if 1, we do not adjust W(:,:,i) during learning [0]
wolffd@0: % cov_prior_weight - weight given to I prior for estimating Sigma [0.01]
wolffd@0: % cov_prior_entropic - if 1, we also use an entropic prior for Sigma [0]
wolffd@0: %
wolffd@0: % e.g., CPD = gaussian_CPD(bnet, i, 'mean', [0; 0], 'clamp_mean', 1)
wolffd@0: 
wolffd@0: if nargin==0
wolffd@0:   % This occurs if we are trying to load an object from a file.
wolffd@0:   CPD = init_fields;
wolffd@0:   clamp = 0;
wolffd@0:   CPD = class(CPD, 'gaussian_CPD', generic_CPD(clamp));
wolffd@0:   return;
wolffd@0: elseif isa(bnet, 'gaussian_CPD')
wolffd@0:   % This might occur if we are copying an object.
wolffd@0:   CPD = bnet;
wolffd@0:   return;
wolffd@0: end
wolffd@0: CPD = init_fields;
wolffd@0:  
wolffd@0: CPD = class(CPD, 'gaussian_CPD', generic_CPD(0));
wolffd@0: 
wolffd@0: args = varargin;
wolffd@0: ns = bnet.node_sizes;
wolffd@0: ps = parents(bnet.dag, self);
wolffd@0: dps = myintersect(ps, bnet.dnodes);
wolffd@0: cps = myintersect(ps, bnet.cnodes);
wolffd@0: fam_sz = ns([ps self]);
wolffd@0: 
wolffd@0: CPD.self = self;
wolffd@0: CPD.sizes = fam_sz;
wolffd@0: 
wolffd@0: % Figure out which (if any) of the parents are discrete, and which cts, and how big they are
wolffd@0: % dps = discrete parents, cps = cts parents
wolffd@0: CPD.cps = find_equiv_posns(cps, ps); % cts parent index
wolffd@0: CPD.dps = find_equiv_posns(dps, ps);
wolffd@0: ss = fam_sz(end);
wolffd@0: psz = fam_sz(1:end-1);
wolffd@0: dpsz = prod(psz(CPD.dps));
wolffd@0: cpsz = sum(psz(CPD.cps));
wolffd@0: 
wolffd@0: % set default params
wolffd@0: CPD.mean = randn(ss, dpsz);
wolffd@0: CPD.cov = 100*repmat(eye(ss), [1 1 dpsz]);    
wolffd@0: CPD.weights = randn(ss, cpsz, dpsz);
wolffd@0: CPD.cov_type = 'full';
wolffd@0: CPD.tied_cov = 0;
wolffd@0: CPD.clamped_mean = 0;
wolffd@0: CPD.clamped_cov = 0;
wolffd@0: CPD.clamped_weights = 0;
wolffd@0: CPD.cov_prior_weight = 0.01;
wolffd@0: CPD.cov_prior_entropic = 0;
wolffd@0: nargs = length(args);
wolffd@0: if nargs > 0
wolffd@0:   CPD = set_fields(CPD, args{:});
wolffd@0: end
wolffd@0: 
wolffd@0: % Make sure the matrices have 1 dimension per discrete parent.
wolffd@0: % Bug fix due to Xuejing Sun 3/6/01
wolffd@0: CPD.mean = myreshape(CPD.mean, [ss ns(dps)]);
wolffd@0: CPD.cov = myreshape(CPD.cov, [ss ss ns(dps)]);
wolffd@0: CPD.weights = myreshape(CPD.weights, [ss cpsz ns(dps)]);
wolffd@0: 
wolffd@0: % Precompute indices into block structured  matrices
wolffd@0: % to speed up CPD_to_lambda_msg and CPD_to_pi
wolffd@0: cpsizes = CPD.sizes(CPD.cps);
wolffd@0: CPD.cps_block_ndx = cell(1, length(cps));
wolffd@0: for i=1:length(cps)
wolffd@0:   CPD.cps_block_ndx{i} = block(i, cpsizes);
wolffd@0: end
wolffd@0: 
wolffd@0: %%%%%%%%%%% 
wolffd@0: % Learning stuff
wolffd@0: 
wolffd@0: % expected sufficient statistics 
wolffd@0: CPD.Wsum = zeros(dpsz,1);
wolffd@0: CPD.WYsum = zeros(ss, dpsz);
wolffd@0: CPD.WXsum = zeros(cpsz, dpsz);
wolffd@0: CPD.WYYsum = zeros(ss, ss, dpsz);
wolffd@0: CPD.WXXsum = zeros(cpsz, cpsz, dpsz);
wolffd@0: CPD.WXYsum = zeros(cpsz, ss, dpsz);
wolffd@0: 
wolffd@0: % For BIC
wolffd@0: CPD.nsamples = 0;
wolffd@0: switch CPD.cov_type
wolffd@0:  case 'full',
wolffd@0:   % since symmetric 
wolffd@0:     %ncov_params = ss*(ss-1)/2; 
wolffd@0:     ncov_params = ss*(ss+1)/2; 
wolffd@0:   case 'diag',
wolffd@0:     ncov_params = ss;
wolffd@0:   otherwise
wolffd@0:     error(['unrecognized cov_type ' cov_type]);
wolffd@0: end
wolffd@0: % params = weights + mean + cov
wolffd@0: if CPD.tied_cov
wolffd@0:   CPD.nparams = ss*cpsz*dpsz + ss*dpsz + ncov_params;
wolffd@0: else
wolffd@0:   CPD.nparams = ss*cpsz*dpsz + ss*dpsz + dpsz*ncov_params;
wolffd@0: end
wolffd@0: 
wolffd@0: % for speeding up maximize_params
wolffd@0: CPD.useC = exist('rep_mult');
wolffd@0: 
wolffd@0: clamped = CPD.clamped_mean & CPD.clamped_cov & CPD.clamped_weights;
wolffd@0: CPD = set_clamped(CPD, clamped);
wolffd@0: 
wolffd@0: %%%%%%%%%%%
wolffd@0: 
wolffd@0: function CPD = init_fields()
wolffd@0: % This ensures we define the fields in the same order 
wolffd@0: % no matter whether we load an object from a file,
wolffd@0: % or create it from scratch. (Matlab requires this.)
wolffd@0: 
wolffd@0: CPD.self = [];
wolffd@0: CPD.sizes = [];
wolffd@0: CPD.cps = [];
wolffd@0: CPD.dps = [];
wolffd@0: CPD.mean = [];
wolffd@0: CPD.cov = [];
wolffd@0: CPD.weights = [];
wolffd@0: CPD.clamped_mean = [];
wolffd@0: CPD.clamped_cov = [];
wolffd@0: CPD.clamped_weights = [];
wolffd@0: CPD.cov_type = [];
wolffd@0: CPD.tied_cov = [];
wolffd@0: CPD.Wsum = [];
wolffd@0: CPD.WYsum = [];
wolffd@0: CPD.WXsum = [];
wolffd@0: CPD.WYYsum = [];
wolffd@0: CPD.WXXsum = [];
wolffd@0: CPD.WXYsum = [];
wolffd@0: CPD.nsamples = [];
wolffd@0: CPD.nparams = [];            
wolffd@0: CPD.cov_prior_weight = [];
wolffd@0: CPD.cov_prior_entropic = [];
wolffd@0: CPD.useC = [];
wolffd@0: CPD.cps_block_ndx = [];