senate2014_est.py

# Senate estimate 30 June 2014
# Copyright Sam Wang, Princeton University, [email protected]
# GNU license - free use but retain this header.
#
# Python translation 5 July 2014 by Ryan Buckley
#
# Dem/Ind safe: 45
# GOP safe: 46
# 9 races remain: AK, AR, CO, GA, IA, KY, LA, MI, NC
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats

def senate2014_est(margins, dem_safe, GOP_safe, un, bias = 0):
    """ Take polling margins, number of safe dem seats, number of safe GOP seats, un, and bias. Return all_dists, whose i-th element (indexed at 0) is the probability of democrats winning i contested seats """
    num_states = len(margins)
    assert( num_states == 100 - dem_safe - GOP_safe )

    # get the medians of the margins for each state
    margin_lens = np.zeros( len(margins) )
    medians = np.zeros( len(margins) )
    margin_mads = np.zeros( len(margins) )
    for i, state_margin in enumerate(margins):
        margin_lens[i] = len(state_margin)
        medians[i] = np.median(state_margin)
        margin_mads[i] = mad(state_margin)
    # calculate the standard error of the mean for each state
    sem = margin_mads / 0.6745 * 1.2 / np.sqrt(margin_lens)

    all_dists = np.zeros( num_states + 1 )
    all_probs = np.zeros( num_states )
    bias_values = np.arange(bias-2, bias+2.1, 0.2)

    for diff in bias_values:
        # SUBTLE CHANGE: max() in MATLAB becomes maximum() in NumPy
        z = ( medians + diff ) / np.maximum(np.sqrt( sem * sem + un*un),
                                            np.sqrt(1.0/500/margin_lens) )
        senate_probs = scipy.stats.t.cdf(z, margin_lens)

        # get dist of sums of independent random variables by convolution
        senate_dist = np.array([1-senate_probs[0], senate_probs[0]])
        for prob in senate_probs[1:num_states+1]:
            senate_dist = np.convolve(senate_dist, np.array([1-prob, prob]))

        all_probs = all_probs + senate_probs
        all_dists = all_dists + senate_dist

    all_probs =   all_probs / len(bias_values)
    all_dists = ( all_dists / len(bias_values) ) *100

    mode_dem = dem_safe + np.argmax(all_dists)
    mode_gop = 100 - mode_dem
    prob_dem = 0.88*all_dists[50-dem_safe] + np.sum(all_dists[51-dem_safe:])
    mean_seats = np.sum(all_dists*np.arange(num_states+1) ) / 100 + dem_safe
    print all_probs
    print all_dists
    print mean_seats
    print prob_dem

    arr_to_write = np.array([mode_dem, mode_gop, prob_dem])
    np.savetxt('Sen_estimates_2014.csv', arr_to_write, delimiter="," )
    return all_dists


def mad(arr):
    """ Calculate the mean absolute deviation of the input array. """
    mean = np.mean(arr)
    n = len(arr)
    return np.sum( np.abs(arr - mean) )*1.0/n

def draw_barplot(dem_safe, num_states, all_dists):
    plt.bar( range(dem_safe, dem_safe+num_states+1), all_dists,
                    color='red', align='center')
    plt.bar( range(50, dem_safe + num_states + 1 ),
                all_dists[ 50-dem_safe:num_states+1],
                color='blue', align='center')
    plt.xticks(  range(dem_safe, dem_safe+num_states+1 ) )
    plt.grid()
    plt.xlabel('Democratic/Independent Senate Seats')
    plt.ylabel('Probability (%)')
    plt.show()

def get_margins():
    """ Return a list of poll margins for contested states """
    # polls newest to oldest from HuffPost. last three polls, going
    # back no more than 6 weeks
    # If over 3 polls are available and the oldest polls are near one another in time, take them all
    margins_ak=np.array([5, -5]) # http://elections.huffingtonpost.com/pollster/2014-alaska-senate-sullivan-vs-begich
    margins_ar=np.array([-7, -8, -4]) # http://elections.huffingtonpost.com/pollster/2014-arkansas-senate-cotton-vs-pryor
    margins_co=np.array([1, 4, -2]) # http://elections.huffingtonpost.com/pollster/2014-colorado-senate-gardner-vs-udall
    margins_ga=np.array([2, 3, -5]) # http://elections.huffingtonpost.com/pollster/2014-georgia-senate-perdue-vs-nunn
    margins_ia=np.array([4, -6, -1, -5]) # http://elections.huffingtonpost.com/pollster/2014-iowa-senate-ernst-vs-braley
    margins_ky=np.array([2, 3, -7]) # http://elections.huffingtonpost.com/pollster/2014-kentucky-senate-mcconnell-vs-grimes
    margins_la=np.array([0, -6]) # http://elections.huffingtonpost.com/pollster/2014-louisiana-senate-cassidy-vs-landrieu
    margins_mi=np.array([3, 5, 9]) # http://elections.huffingtonpost.com/pollster/2014-michigan-senate-land-vs-peters
    margins_nc=np.array([1, 5, 7]) # http://elections.huffingtonpost.com/pollster/2014-north-carolina-senate-tillis-vs-hagan

    # Pass the state margins as a *list* of 1D NumPy arrays
    # using a 2D NumPy array might be nicer, but we have to deal with
    # different states having different numbers of polls
    return [margins_ak, margins_ar, margins_co, margins_ga, margins_ia,
               margins_ky, margins_la, margins_mi, margins_nc]

def main():
    margins = get_margins()
    num_states = len(margins)
    dem_safe = 45
    GOP_safe = 46

    all_dists = senate2014_est(margins, dem_safe, GOP_safe, un=1.5, bias =0)
    draw_barplot( dem_safe, num_states, all_dists )

    bias_range = np.arange(-5, 5.1, 0.2)
    bias_calc(bias_range, margins, dem_safe, GOP_safe, un=1.5)

def bias_calc(bias_range, margins, dem_safe, GOP_safe, un):
    num_states = len(margins)
    assert( num_states == 100 - dem_safe - GOP_safe )

    senate_means = np.zeros( len(bias_range) )
    for i, bias in enumerate(bias_range):
        all_dists = senate2014_est(margins, dem_safe,GOP_safe,un, bias=bias)
        senate_means[i] = ( np.sum( all_dists * np.arange(num_states+1) )
                                    / 100 + dem_safe )
    plt.plot(bias_range, senate_means,'.-k')
    plt.grid()
    plt.xlabel('Shift in national D-R vote margin (%)')
    plt.ylabel('Average Dem+Ind Senate seats')
    plt.axis( [-3, 3, 47.9, 53] )
    plt.show()

if __name__ == "__main__":
    main()