[k_solutions,k_probvalues] = Find_kMPEs(bnet,k,Card)
Find_kMPEs: Given a Bayesian network, find the k most probable
configurations of the network.
This is essentially the algorithm introduced by Nilsson in:
D. Nilsson. An efficient algorithm for finding the {M} most probable configurations in probabilistic expert systems},
Statistics and Computing, Vol 2., 1998, Pp. 159--173
The algorithm computes the junction tree of the BN and apply max-propagation and dynamic programming
for finding the k MPCs. Nilsson proposed two schedules for
finding the subsequent maxima. This implementation
corresponds to the first proposed schedule. Currently, it
works for binary variables
INPUTS
bnet: Bayesian network
k: Number of most probable configurations to find
Card: Cardinality of the variables
OUTPUTS
k_solutions: List of k most probable configurations
k_probvalues: Probalities given to each configuration during the process
(Not necessarily the exact probabilities given by the original JT)
Last version 8/26/2008. Roberto Santana (roberto.santana@ehu.es)0001 function[k_solutions,k_probvalues] = Find_kMPEs(bnet,k,Card) 0002 % [k_solutions,k_probvalues] = Find_kMPEs(bnet,k,Card) 0003 % Find_kMPEs: Given a Bayesian network, find the k most probable 0004 % configurations of the network. 0005 % This is essentially the algorithm introduced by Nilsson in: 0006 % D. Nilsson. An efficient algorithm for finding the {M} most probable configurations in probabilistic expert systems}, 0007 % Statistics and Computing, Vol 2., 1998, Pp. 159--173 0008 % The algorithm computes the junction tree of the BN and apply max-propagation and dynamic programming 0009 % for finding the k MPCs. Nilsson proposed two schedules for 0010 % finding the subsequent maxima. This implementation 0011 % corresponds to the first proposed schedule. Currently, it 0012 % works for binary variables 0013 % INPUTS 0014 % bnet: Bayesian network 0015 % k: Number of most probable configurations to find 0016 % Card: Cardinality of the variables 0017 % OUTPUTS 0018 % k_solutions: List of k most probable configurations 0019 % k_probvalues: Probalities given to each configuration during the process 0020 % (Not necessarily the exact probabilities given by the original JT) 0021 % Last version 8/26/2008. Roberto Santana (roberto.santana@ehu.es) 0022 0023 n = size(Card,2); % Number of variables 0024 [AccCard] = FindAccCard(n,Card); 0025 0026 JT = jtree_inf_engine(bnet, 'maximize', 1); % The original JT is found from the BN 0027 for i=1:n, 0028 auxvar{i} = []; % All variables are hidden 0029 end, 0030 0031 mpe_solution = find_mpe(JT,auxvar); % The most probable configuration is computed 0032 [JT, logprob] = enter_evidence(JT, mpe_solution); % The probability is computed 0033 prob_value = exp(logprob); 0034 0035 DataJTs(1,:) = [0,0,prob_value]; % The fields of DataJTs are: index,marca,prob_value of the JT 0036 AllJTs{1} = JT; 0037 AllConfs{1} = mpe_solution; 0038 nJT = 1; 0039 first = 1; 0040 no_visited = 1; 0041 while ( first <= k & ~isempty(no_visited)) % The procedure ends when the next JT is k or all JTs have been visited 0042 % i.e. many potential configurations have probability zero 0043 no_visited = find(DataJTs(:,1)==0); % All JTs that have not been visited 0044 if(~isempty(no_visited)) 0045 [val,Best_JT_index] = max(DataJTs(no_visited,3)); % The index of the JT with highest probability 0046 0047 k_probvalues(first) = val; 0048 k_solutions(first,:) = AllConfs{no_visited(Best_JT_index)}(:); 0049 Best_JT = AllJTs{no_visited(Best_JT_index)}; % The JT with highest probability 0050 Marca = DataJTs(no_visited(Best_JT_index),2); 0051 DataJTs(no_visited(Best_JT_index),1) = 1; % Now it has been visited 0052 0053 for i=1:n, 0054 auxvar{i} = []; % All variables are hidden 0055 end, 0056 0057 best_conf = k_solutions(first,:); %find_mpe(Best_JT,auxvar); 0058 [valindex] = NumconvertCard(cell2num(k_solutions(first,:))-1,n,AccCard)+1; 0059 0060 0061 for i=Marca+1:n 0062 0063 NewJT = Best_JT; 0064 0065 % To extend to the discrete case, the next steps should be 0066 % modified to enter not only positive envidence (p(x_i)=1) 0067 % but also negative evidence (p(x_i)=0) 0068 for j=1:i, 0069 if (j<i) 0070 auxvar{j} = cell2num(best_conf(j)); 0071 else 0072 auxvar{j} = (3 - cell2num(best_conf(j))); 0073 end 0074 end 0075 0076 mpe_solution = find_mpe(NewJT,auxvar); % The most probable configuration is computed 0077 [NewJT, logprob] = enter_evidence(NewJT, mpe_solution); % The probability is computed 0078 prob_value = exp(logprob); 0079 0080 if prob_value>0 % Only configurations with positive probabilities are considered 0081 nJT = nJT + 1; 0082 AllConfs{nJT} = mpe_solution; 0083 AllJTs{nJT} = NewJT; % The JT has been updated with the evidence 0084 DataJTs(nJT,:) = [0,i,prob_value]; 0085 end 0086 end, 0087 first = first + 1; 0088 end 0089 end 0090 0091 k_solutions = cell2num(k_solutions); 0092 return 0093