matthiasm@8: function score = score_dags(data, ns, dags, varargin)
matthiasm@8: % SCORE_DAGS Compute the score of one or more DAGs
matthiasm@8: % score = score_dags(data, ns, dags, varargin)
matthiasm@8: %
matthiasm@8: % data{i,m} = value of node i in case m (can be a cell array).
matthiasm@8: % node_sizes(i) is the number of size of node i.
matthiasm@8: % dags{g} is the g'th dag
matthiasm@8: % score(g) is the score of the i'th dag
matthiasm@8: %
matthiasm@8: % The following optional arguments can be specified in the form of name/value pairs:
matthiasm@8: % [default value in brackets]
matthiasm@8: %
matthiasm@8: % scoring_fn - 'bayesian' or 'bic' [ 'bayesian' ]
matthiasm@8: %              Currently, only networks with all tabular nodes support Bayesian scoring.
matthiasm@8: % type       - type{i} is the type of CPD to use for node i, where the type is a string
matthiasm@8: %              of the form 'tabular', 'noisy_or', 'gaussian', etc. [ all cells contain 'tabular' ]
matthiasm@8: % params     - params{i} contains optional arguments passed to the CPD constructor for node i,
matthiasm@8: %              or [] if none.  [ all cells contain {'prior', 1}, meaning use uniform Dirichlet priors ]
matthiasm@8: % discrete   - the list of discrete nodes [ 1:N ]
matthiasm@8: % clamped    - clamped(i,m) = 1 if node i is clamped in case m [ zeros(N, ncases) ]
matthiasm@8: %
matthiasm@8: % e.g., score = score_dags(data, ns, mk_all_dags(n), 'scoring_fn', 'bic', 'params', []);
matthiasm@8: %
matthiasm@8: % If the DAGs have a lot of families in common, we can cache the sufficient statistics,
matthiasm@8: % making this potentially more efficient than scoring the DAGs one at a time.
matthiasm@8: % (Caching is not currently implemented, however.)
matthiasm@8: 
matthiasm@8: [n ncases] = size(data);
matthiasm@8: 
matthiasm@8: % set default params
matthiasm@8: type = cell(1,n);
matthiasm@8: params = cell(1,n);
matthiasm@8: for i=1:n
matthiasm@8:   type{i} = 'tabular';
matthiasm@8:   params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 };
matthiasm@8: end
matthiasm@8: scoring_fn = 'bayesian';
matthiasm@8: discrete = 1:n;
matthiasm@8: clamped = zeros(n, ncases);
matthiasm@8: 
matthiasm@8: args = varargin;
matthiasm@8: nargs = length(args);
matthiasm@8: for i=1:2:nargs
matthiasm@8:   switch args{i},
matthiasm@8:    case 'scoring_fn', scoring_fn = args{i+1};
matthiasm@8:    case 'type',       type = args{i+1}; 
matthiasm@8:    case 'discrete',   discrete = args{i+1}; 
matthiasm@8:    case 'clamped',    clamped = args{i+1}; 
matthiasm@8:    case 'params',     if isempty(args{i+1}), params = cell(1,n); else params = args{i+1};  end
matthiasm@8:   end
matthiasm@8: end
matthiasm@8: 
matthiasm@8: NG = length(dags);
matthiasm@8: score = zeros(1, NG);
matthiasm@8: for g=1:NG
matthiasm@8:   dag = dags{g};
matthiasm@8:   for j=1:n
matthiasm@8:     u = find(clamped(j,:)==0);    
matthiasm@8:     ps = parents(dag, j);
matthiasm@8:     score(g) = score(g) + score_family(j, ps, type{j}, scoring_fn, ns, discrete, data(:,u), params{j});
matthiasm@8:   end
matthiasm@8: end