Julia Codes

Olivier de Mouzon, Thibault Laurent, Michel Lebreton and Dominique Lepelley

This document presents the Julia codes used to obtain the computational results included in the paper "Exploring the Effects on the Electoral College of National and Regional Popular Vote Interstate Compact: An Electoral Engineering Perspective".

Codes used in section 3

Function phi_IC

The function phi_IC computes the value of $\Phi_{IC}(K,L)$ (formula given in section 3.1) for given integers $K$ and $L$:

In [1]:
function phi_IC(L, K)
 Q = round(Int64, (K+1)/2)   
 H = K-L
 cste = lgammavect[H+1] - H*log(2)
 res = 0.0
    if L > Q
        res = 1/sqrt(pi*L)
    else
     for i=(Q-L):(Q-1)
      res = res + exp(cste - lgammavect[i+1] - lgammavect[H-i+1])
     end
     res = res/sqrt(pi*L)
    end
end
Out[1]:
phi_IC (generic function with 1 method)

Function psi_IC

The function psi_IC computes the value of $\Psi_{IC}(K,L)$ (formula given in section 3.2) for given integers $K$ and $L$:

In [2]:
function psi_IC(L, K)
 Q = round(Int64, (K+1)/2)
 H = K-L
 res =  exp(lgammavect[H] - lgammavect[Q] - 
    lgammavect[H-Q+1] - log(sqrt(pi)) - (H-1)*log(2))
end
Out[2]:
psi_IC (generic function with 1 method)

We will compute $\Delta_{IC}(K,L)$ thanks to phi_IC and psi_IC (see formula in section 3.3).

Computation of $\Phi_{IC}(51, L)$, $\Psi_{IC}(51, L)$ and $\Delta_{IC}(51, L)$

In [3]:
# choose a value of K
K = 51

# compute the values of loggamma needed 
lgammavect = lgamma(1:(K+1))

# initialization
Q = round(Int64, (K+1)/2)
inside_IC = zeros(Float64, K)
outside_IC = zeros(Float64, K)
total_IC = zeros(Float64, K)

# computation of phi, psi and delta from 1 to Q
for i=1:(Q-1)
 inside_IC[i] = phi_IC(i, K)
 outside_IC[i] = psi_IC(i, K)
 total_IC[i] = (K-i)*outside_IC[i] + i*inside_IC[i]
end

# computation of phi, psi and delta from Q to K
for i=Q:K
 inside_IC[i] = sqrt(1/(i*pi))
 total_IC[i] = i*inside_IC[i]
end

Plot given in Figure 1

In [7]:
using PyPlot
In [8]:
# The values of L 
L = round(Int64, linspace(1, K, K))

# plot the figure with total_IC 
figure(figsize = (8, 10))
subplot(2, 1, 1)
plot(L, total_IC, linewidth = 2.0, linestyle = "-", label = "total",
    color = "red")
scatter(L[indmin(total_IC)], minimum(total_IC), 35, color = "red")
plot([0; indmax(total_IC.>total_IC[1])], 
    [minimum(total_IC[total_IC.>total_IC[1]]); minimum(total_IC[total_IC.>total_IC[1]])],
    linewidth = 2.0, color = "red", linestyle = "--")
plot([indmax(total_IC.>total_IC[1]); indmax(total_IC.>total_IC[1])], [0; minimum(total_IC[total_IC.>total_IC[1]])],
    linewidth = 2.0, color = "red", linestyle = "--")
plot([indmin(total_IC) ; indmin(total_IC)], [0; minimum(total_IC)], color = "red",
    linewidth = 2.0, linestyle = "--")
# legend
title(L"\Delta_{IC}(51,L)")
ax = gca()
ax[:set_xlim]([0, K])
ax[:set_ylim]([0, maximum(total_IC)])
grid("on")
Mx = matplotlib[:ticker][:MultipleLocator](5) # Define interval of major ticks
ax[:xaxis][:set_major_locator](Mx) # Set interval of major ticks
mx = matplotlib[:ticker][:MultipleLocator](1) # Define interval of minor ticks
ax[:xaxis][:set_minor_locator](mx) # Set interval of minor ticks
legend(loc="upper right", fancybox="true")

# plot the figure with inside_IC and outside_IC
subplot(2,1,2)
plot(L, inside_IC, linewidth = 2.0, linestyle = "-", label = "inside",
    color = "magenta")
plot(L, outside_IC, linewidth = 2.0, linestyle = "-", label = "outside",
    color = "cyan")
scatter(L[indmax(inside_IC)], maximum(inside_IC), 35, color = "magenta")
scatter(L[26], 0, 35, color = "cyan")
plot([indmax(inside_IC); indmax(inside_IC)], [-0.01; maximum(inside_IC)], 
    color = "magenta", linestyle = "--")
plot([0; indmax(inside_IC)], [maximum(inside_IC); maximum(inside_IC)], 
    color = "magenta", linestyle = "--")
plot([26; 26], [-0.01; 0], color = "cyan", linestyle = "--")
# legend
title(L"$\Phi_{IC}(51,L)$ and $\Psi_{IC}(51,L)$")
ax = gca()
ax[:set_xlim]([0, K])
ax[:set_ylim]([-0.01, maximum(inside_IC) + 0.01])
grid("on")
Mx = matplotlib[:ticker][:MultipleLocator](5) # Define interval of major ticks
ax[:xaxis][:set_major_locator](Mx) # Set interval of major ticks

mx = matplotlib[:ticker][:MultipleLocator](1) # Define interval of minor ticks
ax[:xaxis][:set_minor_locator](mx) # Set interval of minor ticks

legend(loc = "upper right", fancybox = "true")
ax = gca()
xlabel("Number of states inside the coalition, L")
Out[8]:
PyObject <matplotlib.text.Text object at 0x7fba94265a50>

Codes used in section 4

Computation of $c_L$

To compute $c_L$ (see section 4.1), we use the formula given in Le Breton et al. (2016) :

$$c_L=\displaystyle\sum_{m=0}^{\frac{L-1}{2}}(-1)^m \exp((L+1)\ln(L)-(L-1)\ln(2L)-\ln\Gamma(L-m+1) -\ln\Gamma(m+1) + (L-1)\ln(L-2m))$$

To compute $c_L$ for $L=1,...,1001$, we need to increase the precision machine with Julia until $mantissa=1024$.

In [9]:
Lmax = 1001

# here we give the number of the mantrix
setprecision(1024)

# if BigFloat
lgammavectBF = ones(BigFloat, Lmax+1)

@time for k in 1:(Lmax+1)
 lgammavectBF[k] = lgamma(BigFloat(k))
end

cL = zeros(BigFloat, Lmax)
cLFloat64 = zeros(Float64, Lmax)

for K=1:Lmax
  if(isodd(K))
    loopmax = round(Int64, (K-1)/2)
   else
    loopmax = round(Int64, K/2)
  end

  cste = (K+1)*log(BigFloat(K)) - (K-1)*log(BigFloat(2*K))
  powerm = ones(Int64, loopmax+1)
  for i=1:(loopmax)
   powerm[i+1] = powerm[i]*(-1)
  end

  for m=0:loopmax
    cL[K] = cL[K] + powerm[m+1]*exp(cste - lgammavectBF[K-m+1] - 
        lgammavectBF[m+1] + (K-1)*log(BigFloat(K-2*m)))
    cLFloat64[K] = convert(Float64, cL[K])
  end
end

  1.979948 seconds (622.74 k allocations: 14.512 MB, 0.34% gc time)

Function phi_IAC_star

The function phi_IAC_star computes the value of $\Phi_{IAC*}(K,L)$ (formula given in section 4.1) for given integers $K$ and $L$:

In [10]:
function phi_IAC_star(L, K)
 Q = round(Int64, (K+1)/2)
 H = K-L
 cste = lgammavect[H+1] - H*log(2)
 res = 0.0
  for i=(Q-L):(Q-1)
    res = res + exp(cste - lgammavect[i+1] - lgammavect[H-i+1])
  end
 res = cLFloat64[L]/L*res
end
Out[10]:
phi_IAC_star (generic function with 1 method)

$\Psi_{IAC*}$ is the same that $\Psi_{IC}$ case up to the multiplicative factor $\frac{1}{2r+1}$.

$\Delta_{IAC*}(K,L)$ is computed thanks $\Phi_{IAC*}$ and $\Psi_{IAC*}$.

Computation of $\Phi_{IAC*}(51, L)$, $\Psi_{IAC*}(51, L)$ and $\Delta_{IAC*}(51, L)$

In [11]:
# choose a value of K
K = 51; 

# compute the vector of loggamma
lgammavect = lgamma(1:(K+1))

#initialization
Q = round(Int64, (K+1)/2)
inside_IACstar = zeros(Float64, K)
outside_IACstar = zeros(Float64, K)
total_IACstar = zeros(Float64, K)

# computation from 1 to Q
for i=1:(Q-1)
 inside_IACstar[i] = phi_IAC_star(i, K);
 outside_IACstar[i] = sqrt(pi)*psi_IC(i, K);
 total_IACstar[i] = (K-i)*outside_IACstar[i] + i*inside_IACstar[i]
end
# computation from Q to K
for i=Q:K
 inside_IACstar[i] = cL[i]/i
 total_IACstar[i] = i*inside_IACstar[i]
end

Plot given in Figure 4

In [12]:
# values of L 
L = round(Int64, linspace(1, K, K))

# plot of total_IACstar
figure(figsize=(8,10))
subplot(2,1,1)
plot(L, total_IACstar, linewidth = 2.0, linestyle = "-", 
label = "total", color = "red")
scatter(L[indmin(total_IACstar)], minimum(total_IACstar), color="red")
plot([indmin(total_IACstar) ; indmin(total_IACstar)], [0; minimum(total_IACstar)], 
 color = "red", linestyle = "--", linewidth = 2.0)
# legend
title(L"\Delta_{IAC*}(51,L)")
ax = gca()
ax[:set_xlim]([0, K])
ax[:set_ylim]([0, maximum(total_IACstar)])
grid("on")
Mx = matplotlib[:ticker][:MultipleLocator](5) # Define interval of major ticks
ax[:xaxis][:set_major_locator](Mx) # Set interval of major ticks
mx = matplotlib[:ticker][:MultipleLocator](1) # Define interval of minor ticks
ax[:xaxis][:set_minor_locator](mx) # Set interval of minor ticks
legend(loc = "upper right", fancybox = "true")

# plot of inside_IACstar and outside_IACstar
subplot(2,1,2)
plot(L, inside_IACstar, linewidth = 2.0, linestyle = "-", 
    label = "inside", color = "magenta")
plot(L, outside_IACstar, linewidth = 2.0, linestyle = "-", 
    label = "outside", color = "cyan")
scatter(L[indmax(inside_IACstar)], maximum(inside_IACstar), 35, color = "magenta")
scatter(L[26], 0, 35, color = "cyan")
plot([indmax(inside_IACstar); indmax(inside_IACstar)], [-0.02; maximum(inside_IACstar)], 
    color = "magenta", linestyle = "--")
plot([0; indmax(inside_IACstar)], [maximum(inside_IACstar); maximum(inside_IACstar)], 
    color = "magenta", linestyle = "--")
plot([26; 26], [-0.02; 0], color = "cyan", linestyle = "--")
title(L"$\Phi_{IAC*}(51,L)$ and $\Psi_{IAC*}(51,L)$")
# legend
ax = gca()
ax[:set_xlim]([0, K])
ax[:set_ylim]([-0.02, maximum(inside_IACstar)+0.01])
legend(loc = "upper right", fancybox = "true")
ax = gca()
xlabel(L"D")
xlabel("Number of states inside the coalition, L")
grid("on")
Mx = matplotlib[:ticker][:MultipleLocator](5) # Define interval of major ticks
ax[:xaxis][:set_major_locator](Mx) # Set interval of major ticks
mx = matplotlib[:ticker][:MultipleLocator](1) # Define interval of minor ticks
ax[:xaxis][:set_minor_locator](mx); # Set interval of minor ticks