BeliefMatrices.nb

Matrix calculus for belief functions

Minh Ha-Duong, 2002, 2003

Version :  25/3/2003

This package implements   Matrix calculus for belief functions, as described by Philippe Smets (2001)
http://iridia.ulb.ac.be/~psmets
Licensed under GPL version 2 or later.

TODO

Add {:Keywords:}
Replace b, pl, bel, q by english full words.
BeliefMatrices::usage
i defaults to Length[Ω]

Public interface

BeginPackage["BeliefMatrices`"]

\(Belief::usage = \*""Belief[m,i] lists the   (2^i)-1   values the belief function defined by the basic belief assignment  m  when the frame of reference has  i  elements. Belief of subset  k  is Belief[m,i][[k]], except for k=0 (the empty subset) in which case belief is 0.""
\)
\(Plausibility::usage = \*""Plausibility[m,i] lists the   (2^i)-1  values of the plausibility function defined by the basic belief assignment  m  when the frame of reference has  i  elements. Plausibility of subset  k  is Plausibility[m,i][[k]], except for k=0 (the empty subset) in which case plausibility is 0.""
\)
\(Commonality::usage = \*""Commonality[m,i] lists the  (2^i)-1  values of the commonality function defined by the  basic belief assignment  m  when the frame of reference has  i  elements. Commonality of subset  k  is Commonality[m,i][[k]], except for k=0 (the empty subset) in which case commonality is 1.""
\)
\(Implicability::usage = \*""Implicability[m,i]  lists the   (2^i)-1  values of the commonality function defined by the  basic belief assignment  m  when the frame of reference has  i  elements. Implicability of subset  k  is  Implicability[m,i][[k]], except for k=0 (the empty subset) in which case implicability is  m[0].""
\)
PignisticProbability::usage = "PignisticProbability[m,i]  lists the  \
i  values of the pignistic probability function defined by the  basic belief \
assignment  m  when the frame of reference has  i  elements. Pignistic \
probability of outcome  k  is PignisticProbability[m,i][[k\
]]."

bba::usage = "Symbol  bba  represents the the basic belief \
assignment (also noted  m  in the litterature)";
b::usage = "Symbol  b  represents the implicability \
function";
q::usage = "Symbol  q  represents the commonality function";
bel::usage = "Symbol  bel  represents the belief function";
\(\(pl::usage = "Symbol  pl  represents the plausibility function";\)\)

TransformMatrix::usage = "TransformMatrix[x,y,i] is the 2^i \
dimensional square matrix that obtains  x  when multiplied with  y. Vectors  \
x  and  y  are ordered according to usage in the belief function theory: In x\
[[k]] , the argument  k  must be read in \
binary as a bitfield representing a subset of {1, ..., i}.\n
Acceptable x and y are:\n
 bba     the basic belief assignment (also noted  m  in the litterature),\n
 b       the implicability function\n
 q       the commonality function\n
 bel     the belief function (with the convention bel[0] = 1)\n
 pl      the plausibility function (with the convention pl[0]=1).";

betPMatrix::usage = "betPMatrix[bba,i] is the i lines, 2^i columns \
matrix that obtains pignistic betting probabilities when multiplied with the \
basic belief assigment vector.";

SpecializationMatrix::usage = "SpecializationMatrix[m,i] is the  \
2^i  dimension square Depsterian specialization matrix associated with the \
basic belief assignment  m. They form a commutative and associative \
familly.";

    
GeneralizationMatrix::usage = "GeneralizationMatrix[m,i] is the  \
2^i  dimension square Depsterian generalization matrix associated with the \
basic belief assignment  m. They form a commutative and associative \
familly.";

    
ConditionatingMatrix::usage = "ConditionatingMatrix[k,i] is the  \
2^i  dimension square Dempsterian specialization matrix conditioning with \
respect to subset  k .";

NonInteractiveConjunction::usage = "NonInteractiveConjunction[m1, \
m2 , i] is the non-interactive conjunction of basic belief assignments m1 and \
m2, where  i  is the cardinality of the reference frame (the result uses  \
m1[k]  , for 0 ≤ k ≤ 2^i - 1 ).";

NonInteractiveDisjunction::usage = "NonInteractiveDisjunction[m1, \
m2 , i] is the non-interactive disjunction of basic belief assignments m1 and \
m2, where  i  is the cardinality of the reference frame (the result uses  \
m1[k]  , for 0 ≤ k ≤ 2^i - 1 ).";

\(\(\(AlphaJunction::usage = "AlphaJunction[m1 , m2 , i , α , type] is the α-junction of basic belief assignments m1 and m2, where  i  is the cardinality of the reference frame (the result uses  m1[k]  , for 0 ≤ k ≤ 2^i - 1 ), α is the junction parameter and type is either And for conjunction, or Or for disjunction.";\)\(
\)
    \)\)
AlphaJunctionMatrix::usage = "AlphaJunctionMatrix[m , i , α \
, type] is the α-junction matrix of basic belief assignment m, where  \
i  is the cardinality of the reference frame (the result uses  m1[k]  , for 0 \
≤ k ≤ 2^i - 1 ), α is the junction parameter and \
type is either And for conjunction, or Or for disjunction.";

\(\(Unspecificity::usage = "Unspecificity[m , i] is the relative unspecificity of basic belief assignment m, where  i  is the cardinality of the reference frame (the result uses  m[k]  , for 0 ≤ k ≤ 2^i - 1 ). Unspecificity is a generalization of cardinality.";\)\)

CautiousConjunction::usage = "CautiousConjunction[m1, m2 , i] is \
the cautious junction of basic belief assignments m1 and m2, where  i  is the \
cardinality of the reference frame (the result uses  m1[k]  , for 0 \
≤ k ≤ 2^i - 1 ).";

Begin["`Private`"]

The matrices

\(\(\(TransformMatrix[b, bba, 0] = {1};\)\(
\)
    \)\)
TransformMatrix[b, bba, i_]  /; 
        i > 0  := 
TransformMatrix[b, bba, 
          i] = 
Join[
Flatten /@ 
            Transpose[{TransformMatrix[b, bba, i - 1], 
                Table[0, {(2^(i - 1))}, {(2^(i - 1))}]}], 
Flatten /@ 
            Transpose[{TransformMatrix[b, bba, i - 1], 
                TransformMatrix[b, bba, 
                  i - 1]}]
]

    
\(\(\(TransformMatrix[j, i_] := \(TransformMatrix[j, i] = 
        Reverse[IdentityMatrix[2\^i]]\)\)\(
\)
    \)\)
TransformMatrix[q, bba, i_] := TransformMatrix[q, bba, i] = 
        TransformMatrix[j, i] .  TransformMatrix[b, bba, i] .  
          TransformMatrix[j, i]

    
\(\(\(TransformMatrix[b, q, i_] := \(TransformMatrix[b, q, i] = 
        TransformMatrix[b, bba, i] . \ TransformMatrix[bba, q, i]\)\)\(
\)
    \)\)
\(\(\(TransformMatrix[bel, bba, i_] := \(TransformMatrix[bel, bba, i] = \ 
        ReplacePart[TransformMatrix[b, bba, i] /. {1, l__} -> {0, l}, \ 
          1, {1, 1}]\)\)\(
\)
    \)\)
TransformMatrix[pl, bba, i_] := TransformMatrix[pl, bba, i] = 
        ReplacePart[( 
              Table[1, {(2^i)}, {(2^i)}] - 
                TransformMatrix[j, i]  . 
                  TransformMatrix[bel, bba, i])  /. {1, l__} -> {0, 
                l} , 1, {1, 1}]

    
TransformMatrix[bba, b, i_] := TransformMatrix[bba, b, i] = 
        Inverse[TransformMatrix[b, bba, i]]
TransformMatrix[bba, q, i_] := TransformMatrix[bba, q, i] = 
        Inverse[TransformMatrix[q, bba, i]]
\(\(\(TransformMatrix[q, b, i_] := \(TransformMatrix[q, b, i] = 
        Inverse[TransformMatrix[b, q, i]]\)\)\(
\)
    \)\)
TransformMatrix[bba, pl, i_] := TransformMatrix[bba, pl, i] = 
        Inverse[TransformMatrix[pl, bba, i]]
TransformMatrix[bba, bel, i_] := TransformMatrix[bba, bel, i] = 
        Inverse[TransformMatrix[bel, bba, i]]

Compare with Smets (2001) Tables 4, 5 and 6 page 7:
Matrices for bel and pl are modified so that they are no more singular

{TransformMatrix[b, bba, 3] // MatrixForm,  
      TransformMatrix[bba, b, 3] // MatrixForm}
{TransformMatrix[q, bba, 3] // MatrixForm,  
      TransformMatrix[bba, q, 3] // MatrixForm}
{TransformMatrix[b, q, 3] // MatrixForm,  
      TransformMatrix[q, b, 3] // MatrixForm}
{TransformMatrix[bel, bba, 3] // MatrixForm, 
      TransformMatrix[bba, bel, 3] // MatrixForm}
{TransformMatrix[pl, bba, 3] // MatrixForm, 
      TransformMatrix[bba, pl, 3] // MatrixForm}

Pignistic probabilities

supS[i_] := 
      Transpose[
        Table[Reverse[PadLeft[IntegerDigits[k, 2], i]], {k, 0, 
            (2^i) - 1}]]

cardAInv[i_] := cardAInv[i] = 
        Join[{0}, 
          Table[(1/DigitCount)[k, 2, 1], {k, 1, 
              (2^i) - 1}]]

betPMatrix[bba, i_] := 
      supS[i] . DiagonalMatrix[cardAInv[i]]

Check.
(note: need to do also the probability -> possibility transformation. See Smets (2000) theorem 7.1...)

supS[3] // MatrixForm
cardAInv[3]
betPMatrix[bba, 3] // MatrixForm

TODO also the probability -> possibility transformation. See Smets (2000) theorem 7.1...)

Direct access to belief functions

Clear[vec]
vec[bba, m_, i_] := 
      Table[m[k], {k, 0, (2^i) - 1}]
vec[x_, m_, i_] := 
      TransformMatrix[x, bba, i] . vec[bba, m, i]

Belief[m_, i_] := vec[bel, m, i] // Rest
Plausibility[m_, i_] := vec[pl, m, i] // Rest
Commonality[m_, i_] := vec[q, m, i] // Rest
Implicability[m_, i_] := vec[b, m, i] // Rest
PignisticProbability[m_,  i_] := 
      betPMatrix[bba, i] . vec[bba, m, i]

Generalization, specialization and conditioning

\(\(\(SpecializationMatrix[m_, i_] := 
      TransformMatrix[bba, q, i] . DiagonalMatrix[vec[q, m, i]] . 
        TransformMatrix[q, bba, i]\)\(
\)
    \)\)
GeneralizationMatrix[m_, i_] := 
      TransformMatrix[bba, b, i] . DiagonalMatrix[vec[b, m, i]] . 
        TransformMatrix[b, bba, i]

characteristicBba[event_][x_] := 
      KroneckerDelta[event, x]

    
ConditionatingMatrix[k_, i_] := 
      SpecializationMatrix[characteristicBba[k], i]

Examples from Smets (2001) tables 10 and 11:

vec[bba, mSmets, 3]  = {0.10, 0.05, 0.20, 0.15, 0.27, 0.08, 0.13, 
          0.02};
SpecializationMatrix[mSmets, 3] // Chop // 
      MatrixForm
GeneralizationMatrix[mSmets, 3] // Chop // MatrixForm

Example, conditioning wrt to {a,b} when Ω={a,b,c} see Smets (2001) p. 15.

ConditionatingMatrix[3, 3] // MatrixForm

Non-interactive basic combination

NonInteractiveConjunction[m1_,  m2_, i_] := 
      SpecializationMatrix[m1, i] . 
        vec[bba, m2, i]
NonInteractiveDisjunction[m1_,  m2_, i_] := 
      GeneralizationMatrix[m1, i] . vec[bba, m2, i]

Alpha-junctions

\(\(\(AlphaJunction[m1_, \ m2_, i_, α_, \ type_\ ] := 
      AlphaJunctionMatrix[m1, i, α, \ type] . vec[bba, m2, i]\)\(
\)
    \)\)
AlphaJunctionMatrix[m_, i_, α_,  type_] := 
      Sum[m[k] KMatrix[k, i, α, type],  {k, 0, 
          (2^i) - 1}]

bitUnion[a_, b_] := BitOr[a, b]
bitIntersection[a_, b_] := BitAnd[a, b]
bitSubstraction[a_, b_] := b - a

      
bitInQ[x_, a_] := BitAnd[x, a] ≠ 0
bitCardinal[x_] := 
      Plus @@ RealDigits[x, 
            2][[1]]      
singletonQ[x_] := bitCardinal[x] == 1

KMatrix[compl[x_], i_, α_,  And] /;  
        singletonQ[x] := 
Table[
        Which[
(b == 
                bitUnion[a, x])  &&  (!  bitInQ[x, a]), 
          1, 
(b == a)  &&  (!  
                bitInQ[x, b]), α, 
(a == 
                bitUnion[b, x])  &&  (!  bitInQ[x, b]), 
          1 - α, 
True,  0], 
{a, 
          0, (2^i) - 1}, {b, 0, (2^i) - 1}]

    
KMatrix[j_, i_, α_, And] := 
 
      IdentityMatrix[(2^i)]  . 
        Dot @@ (KMatrix[compl[#], i, α, And] & /@ 
              Select[Table[
                  (2^k), {k, 0, i - 1}], ! bitInQ[#, j] &])

KMatrix[x_, i_, α_, Or] /;  
        singletonQ[x] := 
Table[
        Which[
(a == 
                bitUnion[b, x])  &&  (!  bitInQ[x, b]), 
          1, 
(b == a)  &&  (bitInQ[x, 
                b]), α, 
(b == 
                bitUnion[a, x])  &&  (!  bitInQ[x, a]), 
          1 - α, 
True,  0], 
{a, 
          0, (2^i) - 1}, {b, 0, (2^i) - 1}]

    
KMatrix[j_, i_, α_, Or] := 
 
      IdentityMatrix[(2^i)]  . 
        Dot @@ (KMatrix[#, i, α, Or] & /@ 
              Select[Table[(2^k), {k, 0, i - 1}], 
                bitInQ[#, j] &])

Example 12.1 from Smets p. 21

\(Table[
      MatrixForm[
        KMatrix[k, 2, α, 
            And] /. {1 - α -> α\&_}], {k, 0, 
        3}]\)

Check that α=1 corresponds to specialization

MatrixForm[AlphaJunctionMatrix[m, 2, 1, And]] == 
      MatrixForm[SpecializationMatrix[m, 2]]

Example 12.2 from Smets

\(Table[
      MatrixForm[
        KMatrix[k, 2, α, 
            Or] /. {1 - α -> α\&_}], {k, 0, 3}]\)

Check that  α=1 corresponds to generalization

MatrixForm[AlphaJunctionMatrix[m, 3, 1, Or]] == 
      MatrixForm[GeneralizationMatrix[m, 3]]

Uncertainty criteria

Unspecificity[m_, i_] := 
      vec[bba, m, i] . Table[ bitCardinal[k], {k, 0, (2^i) - 1}]

This unspecificity generalizes cardinality. As discussed in Rocha (1997) Relative uncertainty and evidence sets, it could be divided by Lengh[Ω] to give an index of relative unspecificity. Absolut unspecificity uses Log[Length[subset[k]]] instead.

Unspecificity[m, 3]
Unspecificity[mSmets, 3] == 1.3

Cautious combination using maximum unspecificity

CautiousConjunction[m1_,  m2_, i_] := 
      vec[bba, m, 
          i] /. cautiousConjunctionProgram[m1, m2, 
            i][[2]]

    
cautiousConjunctionProgram[m1_,  m2_, 
          i_] := 
optimum = 
          ConstrainedMax[
Unspecificity[m, 
              i], 
constraintList[m, m1, m2, mm1, mm2, 
              i], 
Join[vec[bba, m, i], vec[bba, mm1, i], 
              vec[bba, mm2, 
                i]]
];

    
constraintList[m_, m1_, m2_, mm1_, mm2_, i_] := 
      Flatten[{Table[
            mm1[k] ≥ 0,  {k, 0, 
              (2^i) - 1}], 
Table[
            mm2[k] ≥ 0,  {k, 0, 
              (2^i) - 1}], 
Sum[
              mm1[k],  {k, 0, (2^i) - 1}]  ==  
            1 , 
Sum[
              mm2[k],  {k, 0, (2^i) - 1}]  ==  1 , 
          Transpose [{vec[bba, m, i], 
                nonInteractiveConjunction[m1,  mm1, i]}] /. 
            List[a_, b_] -> 
              Equal[a, b], 
Transpose [{vec[bba, m, i], 
                nonInteractiveConjunction[m2,  mm2, i]}] /. 
            List[a_, b_] -> Equal[a, b]
}]

Check that cautious combination is idempotent

Do[mSmets[i] = vec[bba, mSmets, 3][[
          i + 1]], {i, 0, 7}]
vec[bba, mSmets, 3]

    
CautiousConjunction[mSmets, mSmets, 3]

Hyper Cautious combination when m is a possibility function.

hyperCautiousConjunction[m1_,  m2_, 
        i_] := 
TransformMatrix[bba, q, 
          i] . (Transpose[{vec[q, m1, i], vec[q, m2, i]}] /. 
            List[a_, b_] -> Min[a, b])

Example: Smets (2000) QEPT seen as Hyper Cautious TBM tables 1-4
(there are two errors in the paper)

mO[_] = 0;
mO[2] = 0.1;
mO[3] = 0.6;
mO[7] = 0.3;

    
mT[_] = 0;
mT[4] = 0.2;
mT[5] = 0.2;
mT[7] = 0.6;

    
hyperCautiousConjunction[mO, mT, 3]

    
TransformMatrix[q, bba, 3] . %

End[]

EndPackage[]