|
wolffd@0
|
1 function prob = quickscore(fpos, fneg, inhibit, prior, leak)
|
|
wolffd@0
|
2 % QUICKSCORE Heckerman's algorithm for BN2O networks.
|
|
wolffd@0
|
3 % prob = quickscore(fpos, fneg, inhibit, prior, leak)
|
|
wolffd@0
|
4 %
|
|
wolffd@0
|
5 % Consider a BN2O (Binary Node 2-layer Noisy-or) network such as QMR with
|
|
wolffd@0
|
6 % dieases on the top and findings on the bottom. (We assume all findings are observed,
|
|
wolffd@0
|
7 % since hidden leaves can be marginalized away.)
|
|
wolffd@0
|
8 % This algorithm takes O(2^|fpos|) time to compute the marginal on all the diseases.
|
|
wolffd@0
|
9 %
|
|
wolffd@0
|
10 % Inputs:
|
|
wolffd@0
|
11 % fpos = the positive findings (a vector of numbers in {1, ..., Nfindings})
|
|
wolffd@0
|
12 % fneg = the negative findings (a vector of numbers in {1, ..., Nfindings})
|
|
wolffd@0
|
13 % inhibit(i,j) = inhibition prob. for finding i, disease j, or 1.0 if j is not a parent.
|
|
wolffd@0
|
14 % prior(j) = prior prob. disease j is ON. We assume prior(off) = 1-prior(on).
|
|
wolffd@0
|
15 % leak(i) = inhibition prob. for the leak node for finding i
|
|
wolffd@0
|
16 %
|
|
wolffd@0
|
17 % Output:
|
|
wolffd@0
|
18 % prob(d) = Pr(disease d = on | ev)
|
|
wolffd@0
|
19 %
|
|
wolffd@0
|
20 % For details, see
|
|
wolffd@0
|
21 % - Heckerman, "A tractable inference algorithm for diagnosing multiple diseases", UAI89.
|
|
wolffd@0
|
22 % - Rish and Dechter, "On the impact of causal independence", UCI tech report, 1998.
|
|
wolffd@0
|
23 %
|
|
wolffd@0
|
24 % Note that this algorithm is numerically unstable, since it adds a large number of positive and
|
|
wolffd@0
|
25 % negative terms and hopes that some of them exactly cancel.
|
|
wolffd@0
|
26 %
|
|
wolffd@0
|
27 % For matlab experts, use 'mex' to compile C_quickscore, which has identical behavior to this function.
|
|
wolffd@0
|
28
|
|
wolffd@0
|
29 [nfindings ndiseases] = size(inhibit);
|
|
wolffd@0
|
30
|
|
wolffd@0
|
31 % make the first disease be always on, for the leak term
|
|
wolffd@0
|
32 Pon = [1 prior(:)'];
|
|
wolffd@0
|
33 Poff = 1-Pon;
|
|
wolffd@0
|
34 Uon = [leak(:) inhibit]; % U(f,d) = Pr(f=0|d=1)
|
|
wolffd@0
|
35 Uoff = [leak(:) ones(nfindings, ndiseases)]; % Uoff(f,d) = Pr(f=0|d=0)
|
|
wolffd@0
|
36 ndiseases = ndiseases + 1;
|
|
wolffd@0
|
37
|
|
wolffd@0
|
38 npos = length(fpos);
|
|
wolffd@0
|
39 post = zeros(ndiseases, 2);
|
|
wolffd@0
|
40 % post(d,1) = alpha Pr(d=off), post(d,2) = alpha Pr(d=m)
|
|
wolffd@0
|
41
|
|
wolffd@0
|
42 FP = length(fpos);
|
|
wolffd@0
|
43 %allbits = logical(dec2bitv(0:(2^FP - 1), FP));
|
|
wolffd@0
|
44 allbits = logical(ind2subv(2*ones(1,FP), 1:(2^FP))-1);
|
|
wolffd@0
|
45
|
|
wolffd@0
|
46 for si=1:2^FP
|
|
wolffd@0
|
47 bits = allbits(si,:);
|
|
wolffd@0
|
48 fprime = fpos(bits);
|
|
wolffd@0
|
49 fmask = zeros(1, nfindings);
|
|
wolffd@0
|
50 fmask(fneg)=1;
|
|
wolffd@0
|
51 fmask(fprime)=1;
|
|
wolffd@0
|
52 fmask = logical(fmask);
|
|
wolffd@0
|
53 p = 1;
|
|
wolffd@0
|
54 pterm = zeros(1, ndiseases);
|
|
wolffd@0
|
55 ptermOff = zeros(1, ndiseases);
|
|
wolffd@0
|
56 ptermOn = zeros(1, ndiseases);
|
|
wolffd@0
|
57 for d=1:ndiseases
|
|
wolffd@0
|
58 ptermOff(d) = prod(Uoff(fmask,d));
|
|
wolffd@0
|
59 ptermOn(d) = prod(Uon(fmask,d));
|
|
wolffd@0
|
60 pterm(d) = Poff(d)*ptermOff(d) + Pon(d)*ptermOn(d);
|
|
wolffd@0
|
61 end
|
|
wolffd@0
|
62 p = prod(pterm);
|
|
wolffd@0
|
63 sign = (-1)^(length(fprime));
|
|
wolffd@0
|
64 for d=1:ndiseases
|
|
wolffd@0
|
65 myp = p / pterm(d);
|
|
wolffd@0
|
66 post(d,1) = post(d,1) + sign*(myp * ptermOff(d));
|
|
wolffd@0
|
67 post(d,2) = post(d,2) + sign*(myp * ptermOn(d));
|
|
wolffd@0
|
68 end
|
|
wolffd@0
|
69 end
|
|
wolffd@0
|
70
|
|
wolffd@0
|
71 post(:,1) = post(:,1) .* Poff(:);
|
|
wolffd@0
|
72 post(:,2) = post(:,2) .* Pon(:);
|
|
wolffd@0
|
73 post = mk_stochastic(post);
|
|
wolffd@0
|
74 prob = post(2:end,2)'; % skip the leak term
|
|
wolffd@0
|
75
|
|
wolffd@0
|
76
|
