This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Numerics Volume, Chapter 1

The Banzhaf index of a voting system is defined in the following way [✶], [✶]: Given n parties (Subscript[p, k]), each having a (nonnegative) voting weight (Subscript[v, k]) (not necessary integer), (2^n)-1 possible coalitions can be formed through all subsets of the n parties. A coalition is winning, if its combined voting weight is larger than a given quota q (typically 50% or 66.67% of the sums of all weights (Subscript[ν, k])). Let ω be the total number of winning coalitions that can be formed, and let (Subscript[ω, k]) be the number of coaltions that contains the party (Subscript[p, k]). Then the Banzhaf index (Subscript[η, k]) of the party (Subscript[p, k]) is defined as (Subscript[η, k])=2(Subscript[ω, k])-ω. The normalized Banzhaf index (Subscript[(OverTilde[η]), k]) is (Subscript[(OverTilde[η]), k])=(Subscript[(OverTilde[η]), k])/(2^(n-1)). Implement a function that calculates the normalized Banzhaf indices for given party weights and quotas.

We start with implementing a functions that counts the winning coalitions in general and for each party. Because we want later apply our function to a list of 25 parties (meaning we will have to analyze 33554432 coalitions), we use the function Compile, for speed reasons. To avoid excessive memory use, we do not create a list of all possible coalitions and keep it, but rather generate and analyze the possible coalitions sequentially. Each possible coalitions corresponds naturally to the base 2 representation of an integer between 1 and (2^n)-1 (if the kth digits is a one, the party (Subscript[p, k]) is part of the coalition, otherwise not). For each quota, we compare the combined coalition power with the quotas and increment the counters for the winning coalitions. Because Compile cannot return ragged arrays, we define the subsidiary function winningCoalitionDataC that forms the nested lists of total winning coalitions and winning coalitions that the individual parties are members of.

winningCoalitionDataC[partyPowers_List, quotas_List] :=
Module[{resC = winningCoalitionDataMainC[partyPowers, quotas]},
 (* separate total number of winning coalitions from party counters *)
 Transpose[{Take[resC, Length[quotas]], 
            Partition[Drop[resC, Length[quotas]], Length[partyPowers]]}]]

winningCoalitionDataMainC =
Compile[{{partyPowers, _Real, 1}, {quotas, _Real, 1}},
Module[{o = Length[partyPowers], n = Length[quotas], L, l, weight, 
        winningCoalitionCounters, winningParticipatingCoalitionCounters}, 
        l = Array[1&, {o}]; L = Array[1&, {o}];
        (* winning coalition counters *)
        winningCoalitionCounters = Array[0&, {n}];
        winningParticipatingCoalitionCounters = Table[Array[0&, {o}], {n}];
 (* loop of all possible (2^n - 1) coalitions *)
 Do[(* coalition weight *)
    l = IntegerDigits[k, 2, o]; weight = l.partyPowers;
    (* compare colation weight with all quotas *)
    Do[If[weight > quotas[[j]], 
          (* increase coalition counters *)
          winningCoalitionCounters[[j]] =  winningCoalitionCounters[[j]] + 1;
          winningParticipatingCoalitionCounters[[j]] = 
                 winningParticipatingCoalitionCounters[[j]] + l],
       {j, n}],
    {k, 2^o}];
  (* form one list because Compile must return a tensor *)
  Join[winningCoalitionCounters, 
       Flatten[winningParticipatingCoalitionCounters]]],
                                          {{_IntegerDigits, _Integer, 1}}];

The function winningCoalitionDataMainC compiled flawlessly.

Union[Head /@ Flatten[winningCoalitionDataMainC[[4]]]]

Here is a simple example for five parties. Seven winning coalitions can be formed to reach a 2/3 majority.

winningCoalitionDataC[{5, 4, 3, 2, 1}, {10}]

Plus @@ {5, 4, 3, 2, 1}

Here is the list of EU country weights.

weightsEU = (* appriximate square root of population 
               measured in million people *)
            {9.09 (* Germany *), 7.72 (* France *), 7.69 (* U.K. *), 
             7.55 (* Italy *), 6.38 (* Spain *), 6.18 (* Poland *), 
             4.02 (* Netherlands *), 3.32 (* Greece *), 3.23 (* Portugal *), 
             3.22 (* Belgium *), 3.19 (* Czech Republic *), 3.19 (* Hungary *), 
             2.99 (* Sweden *), 2.84 (* Austria *), 2.32 (* Denmark *), 
             2.32 (* Slovakia *), 2.28 (* Finland *), 1.99 (* Ireland *), 
             1.86 (* Lithuania *), 1.53 (* Latvia *), 1.41 (* Slovenia *), 
             1.17 (* Estonia *), 0.89 (* Cyprus *), 0.67 (* Luxembourg *), 
             0.63 (* Malta *)};

To obtain a 51% majority vote, approximately 15 million winning coalitions can be formed. Analyzing all possible 33 million coalitions takes only minutes.

winningCoalitionDataC[weightsEU, {0.51 (Plus @@ weightsEU)}] // Timing

Given the counts of the winning coalitions, it is straightforward to define the Banzhaf indices.

BanzhafIndices[partyPowers_, quotas_] := 
  ((2 #2 - #1)/2^(Length[partyPowers] - 1))& @@@ 
             winningCoalitionDataC[partyPowers, quotas]

BanzhafIndicesNormalized[partyPowers_, quotas_] := 
  #/(Plus @@ #)& /@ BanzhafIndices[partyPowers, quotas]

For a 51% majority vote, the Banzhaf indices deviate substantially from the country weights. The larger countries have too much weight, and the smaller to little.

BanzhafIndicesNormalized[weightsEU, {0.51 (Plus @@ weightsEU)}] // N

To obtain a quota such that the country weights agree approximately with their Banzhaf indices, we will minimize the sum of the squares of the differences between the country weights and their Banzhaf indices. The function votingResiduals defines the square root of these sum of squares.

votingResiduals[partyPowers_, quotas_] :=
Module[{bhN = BanzhafIndicesNormalized[partyPowers, quotas], 
        sqL = #/(Plus @@ #)&[partyPowers]},
       Sqrt[Plus @@ ((# - sqL)^2)]& /@ bhN]

The function votingResidualsTable generates a list of the just-defined sums of squares for pp quotas from the interval [quota0,quota1].

votingResidualsTable[partyPowers_, {quota0_, quota1_}, pp_] :=
Module[{relQuotas = Table[q, {q, quota0, quota1, (quota1 - quota0)/pp}], 
        weightSum = Plus @@ partyPowers},
 Transpose[{relQuotas, votingResiduals[partyPowers, relQuotas weightSum]}]]

For comparison, we 10 times choose randomly 18 party weights that are uniformly, exponentially, and normally distributed. Approximate equality between the given voting power and the Banzhaf indices happens for quotas around 60% … 65%.

SeedRandom[1]

Module[{partyPowers, parties = 18, elections = 10, quotas = 50}, 
Show[GraphicsArray[
 Function[𝓅, Graphics[
  Table[partyPowers = Table[𝓅[Random[]], {parties}];
        (* random color for each party-weight list *)
        {Hue[Random[]], PointSize[0.006], Point /@ 
               votingResidualsTable[partyPowers, {0.5, 0.75}, quotas]}, 
        {elections}], 
       PlotRange -> All, Frame -> True, Axes -> False, Frame -> True]] /@
       (* uniform, exponential, and normal distribution *)
       {Identity, -Log[1 - #]&, Abs[2 + InverseErf[2 # - 1]]&}]]];

For the current EU weights, we obtain approximate equality between the given voting power and the Banzhaf indices for a quota of 62%. Zooming into the minima (right graphic) shows that 62% is an excellent approximation for the optimal quota.

(res = votingResidualsTable[weightsEU, {0.5, 0.75}, 50]); // Timing
resD = votingResidualsTable[weightsEU, {0.615, 0.625}, 25];
Show[GraphicsArray[
      ListPlot[#, PlotRange -> All, Frame -> True, Axes -> False,
               DisplayFunction -> Identity]& /@ {res, resD}]]

Because the voting weights all have two decimal digits, not a single point-minimum for the sum of squares can be found. The next graphic shows a close-up around the minimum. The optimal quota is the interval [0.6201, 0.6202] and the square root of the differences becomes about 0.000309. Also, at this resolution, the sum of squares is no longer a monotonic function to the left and right of the optimal quota.

resD = votingResidualsTable[weightsEU, {0.619, 0.621}, 100];
Show[GraphicsArray[
      ListPlot[#, PlotRange -> All, Frame -> True, Axes -> False,
               DisplayFunction -> Identity]& /@ {resD}]]