python类polyfit()的实例源码

def fit(self):
        # Quadratic fit to get an initial guess for the parameters. 
        # Thanks: https://github.com/materialsproject/pymatgen
        # -> pymatgen/io/abinitio/eos.py 
        a, b, c = np.polyfit(self.volume, self.energy, 2)
        v0 = -b/(2*a)
        e0 = a*v0**2 + b*v0 + c
        b0 = 2*a*v0
        b1 = 4  # b1 is usually a small number like 4
        if not self.volume.min() < v0 and v0 < self.volume.max():
            raise StandardError('The minimum volume of a fitted parabola is not in the input volumes')

        # need to use lst2dct and dct2lst here to keep the order of parameters 
        pp0_dct = dict(e0=e0, b0=b0, b1=b1, v0=v0)
        target = lambda pp, v: self.energy - self.func(v, self.func.lst2dct(pp))
        pp_opt, ierr = leastsq(target, 
                               self.func.dct2lst(pp0_dct), 
                               args=(self.volume,))
        self.params = self.func.lst2dct(pp_opt)

def draw_bs_pairs_linreg(x, y, size=1):
    """Perform pairs bootstrap for linear regression."""

    # Set up array of indices to sample from: inds
    inds = np.arange(len(x))

    # Initialize replicates: bs_slope_reps, bs_intercept_reps
    bs_slope_reps = np.empty(size)
    bs_intercept_reps = np.empty(size)

    # Generate replicates
    for i in range(size):
        bs_inds = np.random.choice(inds, size=len(inds))
        bs_x, bs_y = x[bs_inds], y[bs_inds]
        # noinspection PyTupleAssignmentBalance
        bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x, bs_y, 1)

    return bs_slope_reps, bs_intercept_reps

def SLOPE(self, param):
        class Context:
            def __init__(self, N):
                self.N = N
                self.q = deque([], self.N)
                self.x = [i for i in range(self.N)]

            def handleInput(self, value):
                if len(self.q) < self.N:
                    self.q.append(value)
                    return 0
                self.q.append(value)
                z1 = np.polyfit(self.x, self.q, 1)
                return z1[0]
        ctx = Context(param[1])
        result = param[0].apply(ctx.handleInput)
        return result

def FORCAST(self, param):
        class Context:
            def __init__(self, N):
                self.N = N
                self.q = deque([], self.N)
                self.x = [i for i in range(self.N)]

            def handleInput(self, value):
                if len(self.q) < self.N:
                    self.q.append(value)
                    return np.NaN
                z1 = np.polyfit(self.x, self.q, 1)
                fn = np.poly1d(z1)
                y = fn(self.N + 1)
                self.q.append(value)
                return y

        ctx = Context(param[1])
        result = param[0].apply(ctx.handleInput)
        return result

#????

def minScalErr(stec,el,z,thisBias):
    """
    this determines the slope of the vTEC vs. Elevation line, which
    should be minimized in the minimum scalloping technique for
    receiver bias removal
    inputs:
        stec - time indexed Series of slant TEC values
        el - corresponding elevation values, also Series
        z - mapping function values to convert to vTEC from entire file, may
            contain nans, Series
        thisBias - the bias to be tested and minimized
    """

    intel=np.asarray(el[stec.index],int) # bin the elevation values into int
    sTEC=np.asarray(stec,float)
    zmap = z[stec.index]
    c=np.array([(i,np.average((sTEC[intel==i]-thisBias)
                              /zmap[intel==i])) for i in np.unique(intel) if i>30])

    return np.polyfit(c[:,0],c[:,1],1)[0]

def draw_bs_pairs_linreg(x, y, size=1):
    """Perform pairs bootstrap for linear regression."""
    # Set up array of indices to sample from: inds
    inds = np.arange(len(x))

    # Initialize replicates: bs_slope_reps, bs_intercept_reps
    bs_slope_reps = np.empty(size)
    bs_intercept_reps = np.empty(size)

    # Generate replicates
    for i in range(size):
        bs_inds = np.random.choice(inds, size=len(inds))
        bs_x, bs_y = x[bs_inds], y[bs_inds]
        bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x, bs_y, 1)

    return bs_slope_reps, bs_intercept_reps

def draw_bs_pairs_linreg(x, y, size=1):
    """Perform pairs bootstrap for linear regression."""

    # Set up array of indices to sample from: inds
    inds = np.arange(len(x))

    # Initialize replicates: bs_slope reps, bs_intercept_reps
    bs_slope_reps = np.empty(size)
    bs_intercept_reps = np.empty(size)

    # Generate replicates
    for i in range(size):
        bs_inds = np.random.choice(inds, size=len(inds))
        bs_x, bs_y = x[bs_inds], y[bs_inds]
        bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x, bs_y, 1)

    return bs_slope_reps, bs_intercept_reps

def draw_bs_pairs(x, y, func, size=1):
    """Perform pairs bootstrap for linear regression."""

    # Set up array of indices to sample from: inds
    inds = np.arange(len(x))
    # Initialize replicates
    bs_replicates = np.empty(size)
    # bs_intercept_reps = ____

    # Generate replicates
    for i in range(size):
        bs_inds = np.random.choice(inds, len(inds))
        bs_x, bs_y = x[bs_inds], y[bs_inds]
        bs_replicates[i] = func(bs_x, bs_y)
        # bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x, bs_y, 1)

    return bs_replicates

def test_polyfit_build(self):
        # Ticket #628
        ref = [-1.06123820e-06, 5.70886914e-04, -1.13822012e-01,
               9.95368241e+00, -3.14526520e+02]
        x = [90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
             104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115,
             116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 129,
             130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
             146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157,
             158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
             170, 171, 172, 173, 174, 175, 176]
        y = [9.0, 3.0, 7.0, 4.0, 4.0, 8.0, 6.0, 11.0, 9.0, 8.0, 11.0, 5.0,
             6.0, 5.0, 9.0, 8.0, 6.0, 10.0, 6.0, 10.0, 7.0, 6.0, 6.0, 6.0,
             13.0, 4.0, 9.0, 11.0, 4.0, 5.0, 8.0, 5.0, 7.0, 7.0, 6.0, 12.0,
             7.0, 7.0, 9.0, 4.0, 12.0, 6.0, 6.0, 4.0, 3.0, 9.0, 8.0, 8.0,
             6.0, 7.0, 9.0, 10.0, 6.0, 8.0, 4.0, 7.0, 7.0, 10.0, 8.0, 8.0,
             6.0, 3.0, 8.0, 4.0, 5.0, 7.0, 8.0, 6.0, 6.0, 4.0, 12.0, 9.0,
             8.0, 8.0, 8.0, 6.0, 7.0, 4.0, 4.0, 5.0, 7.0]
        tested = np.polyfit(x, y, 4)
        assert_array_almost_equal(ref, tested)

def calculateKDP(phiDP, ran):
    # smooth phiDP field and take derivative
    numRan = ran.shape[0]
    kdp = np.ma.masked_all(phiDP.shape)
    smoothPhiDP = smPhiDP(phiDP, ran)

    # take derivative of kdp field
    winLen = 11
    rprof = ran[0:winLen*2-1]/1000.

    for i in range(numRan-winLen*3):
        numvalid = smoothPhiDP[:,i:i+winLen*2-1].count(axis=1)
        max_numv = np.max(numvalid)
        if max_numv==(winLen*2-1):
            kdp[numvalid==(winLen*2-1),i+winLen] = 0.5*np.polyfit(rprof,
                smoothPhiDP[numvalid==(winLen*2-1),i:i+winLen*2-1].transpose(), 1)[0]
    return kdp

# function for creating color map
#----------------------------------

def hplot(name, h, err, xlab="", ylab="", hstr="h", rate=None, fit=True, fig=True):
    from matplotlib.pyplot import figure, loglog, xlabel, ylabel, legend, show
    if fig is True:
        figure()
    loglog(h, err, 's-', label=name)
    if fit:
        p = numpy.polyfit(numpy.log(h), numpy.log(err), 1)
        C = numpy.exp(p[1])
        alpha = p[0]
        alg = [C*n**alpha for n in h]
        loglog(h, alg, '--', label="%.3f*%s^(%.3f)" %(C, hstr, alpha))
    if rate and h[0] != 0:
        alg = [err[0]/(h[0]**rate)*n**rate for n in h]
        loglog(h, alg, '--', label="%s^(%.3f)" %(hstr, rate))
    xlabel(xlab)
    ylabel(ylab)
    legend(loc='best')

# --- create interpolation for different h and orders ---

def default1():
    numRadialPts = 144;
    rhoMin = 0.1
    rhoMax = 0.9
    rho = np.linspace(rhoMin, rhoMax, numRadialPts);
    qbar = 0.854 + 2.184 * rho**2;

    minorRadius = 0.28;  # a
    majorRadius = 0.71;  # R0
    r = rho * minorRadius;
    safetyFactor = qbar / np.sqrt(1 - (r/majorRadius)**2);

    # fit
    polynomialDegree = 3;
    p = np.polyfit(rho, safetyFactor, polynomialDegree)
    qCoeffs = p[::-1]
    return qCoeffs

def PCSCH3(self,row):
        actuals=[];read=[]
        resistance = float(self.tbl.item(row,0).text())
        for a in np.linspace(.2e-3,1.5e-3,20):
            actuals.append( self.I.DAC.setCurrent(a) )
            time.sleep(0.001)
            read.append (self.I.get_average_voltage('CH3',samples=5))
        read = np.array(read)/resistance 
        actuals = np.array(actuals)
        sread = read*1e3
        sactuals = actuals*1e3
        self.DacCurves['PCS'].setData(sactuals,sread-sactuals)
        self.tbl.item(row,1).setText(string.join(['%.3f'%a for a in sread-sactuals],' '))
        if np.any(abs(read-actuals)>10e-6):self.setSuccess(self.tbl.item(row,1),0)
        else: self.setSuccess(self.tbl.item(row,1),1)

        fitvals = np.polyfit(self.I.DAC.CHANS['PCS'].VToCode(read),self.I.DAC.CHANS['PCS'].VToCode(actuals),1)
        print (fitvals)
        if list(fitvals):
            self.PCS_SLOPE = fitvals[0] #slope
            self.PCS_OFFSET = fitvals[1]    #offset

def PCSCH3(self,row):
        actuals=[];read=[]
        resistance = float(self.tbl.item(row,0).text())
        for a in np.linspace(.2e-3,1.5e-3,20):
            actuals.append( self.I.DAC.setCurrent(a) )
            time.sleep(0.001)
            read.append (self.I.get_voltage('CH3',samples=5))
        read = np.array(read)/resistance 
        actuals = np.array(actuals)
        sread = read*1e3
        sactuals = actuals*1e3
        self.DacCurves['PCS'].setData(sactuals,sread-sactuals)
        self.tbl.item(row,1).setText(string.join(['%.3f'%a for a in sread-sactuals],' '))
        if np.any(abs(read-actuals)>10e-6):self.setSuccess(self.tbl.item(row,1),0)
        else: self.setSuccess(self.tbl.item(row,1),1)

        fitvals = np.polyfit(self.I.DAC.CHANS['PCS'].VToCode(read),self.I.DAC.CHANS['PCS'].VToCode(actuals),1)
        print (fitvals)
        if list(fitvals):
            self.PCS_SLOPE = fitvals[0] #slope
            self.PCS_OFFSET = fitvals[1]    #offset

def get_linear_regression(xdata, ydata):
    """ Calculate the coefficients of the linear equation corresponding
    to the linear regression of a series of points. """
    try:
        import numpy
        return tuple(numpy.polyfit(xdata, ydata, 1))
    except ImportError:
        # numpy not available
        # try something approximate and simple
        datasize = len(xdata)
        sum_x = float(sum(xdata))
        sum_y = float(sum(ydata))
        sum_xx = float(sum(map(lambda x: x * x, xdata)))
        sum_products = float(sum([xdata[i] * ydata[i]
                                  for i in range(datasize)]))
        a = (sum_products - sum_x * sum_y / datasize) / (
            sum_xx - (sum_x * sum_x) / datasize)
        b = (sum_y - a * sum_x) / datasize
        return a, b

def polyfit_baseline(bands, intensities, poly_order=5, num_stdv=3.,
                     max_iter=200):
    '''Iteratively fits a polynomial, discarding far away points as peaks.
    Similar in spirit to ALS and related methods.
    Automated method for subtraction of fluorescence from biological Raman spectra
    Lieber & Mahadevan-Jansen 2003
    '''
    fit_pts = intensities.copy()
    # precalculate [x^p, x^p-1, ..., x^1, x^0]
    poly_terms = bands[:, None] ** np.arange(poly_order, -1, -1)
    for _ in range(max_iter):
        coefs = np.polyfit(bands, fit_pts.T, poly_order)
        baseline = poly_terms.dot(coefs).T
        diff = fit_pts - baseline
        thresh = diff.std(axis=-1) * num_stdv
        mask = diff > np.array(thresh, copy=False)[..., None]
        unfitted = np.count_nonzero(mask)
        if unfitted == 0:
            break
        fit_pts[mask] = baseline[mask]  # these points are peaks, discard
    else:
        print("Warning: polyfit_baseline didn't converge in %d iters" % max_iter)
    return baseline

def band_center(spectrum, low_endmember=None, high_endmember=None, degree=3):
    x = spectrum.index
    y = spectrum

    if not low_endmember:
        low_endmember = x[0]
    if not high_endmember:
        high_endmember = x[-1]

    ny = y[low_endmember:high_endmember]

    fit = np.polyfit(ny.index, ny, degree)

    center_fit = Series(np.polyval(fit, ny.index), ny.index)
    center = band_minima(center_fit)

    return center, center_fit

def test_polyfit_build(self):
        # Ticket #628
        ref = [-1.06123820e-06, 5.70886914e-04, -1.13822012e-01,
               9.95368241e+00, -3.14526520e+02]
        x = [90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
             104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115,
             116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 129,
             130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
             146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157,
             158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
             170, 171, 172, 173, 174, 175, 176]
        y = [9.0, 3.0, 7.0, 4.0, 4.0, 8.0, 6.0, 11.0, 9.0, 8.0, 11.0, 5.0,
             6.0, 5.0, 9.0, 8.0, 6.0, 10.0, 6.0, 10.0, 7.0, 6.0, 6.0, 6.0,
             13.0, 4.0, 9.0, 11.0, 4.0, 5.0, 8.0, 5.0, 7.0, 7.0, 6.0, 12.0,
             7.0, 7.0, 9.0, 4.0, 12.0, 6.0, 6.0, 4.0, 3.0, 9.0, 8.0, 8.0,
             6.0, 7.0, 9.0, 10.0, 6.0, 8.0, 4.0, 7.0, 7.0, 10.0, 8.0, 8.0,
             6.0, 3.0, 8.0, 4.0, 5.0, 7.0, 8.0, 6.0, 6.0, 4.0, 12.0, 9.0,
             8.0, 8.0, 8.0, 6.0, 7.0, 4.0, 4.0, 5.0, 7.0]
        tested = np.polyfit(x, y, 4)
        assert_array_almost_equal(ref, tested)

def trendLine(self, axis_choose=None):
        stable_sec= int(self.record_sec_le.text())
        stable_count = int(stable_sec * (1/0.007))

        if axis_choose:
            axis = axis_choose
        else:
            axis = str(self.axis_combobox.currentText())

        x = self.raw_data['time'][:stable_count]
        y = self.raw_data[axis][:stable_count]
        coefficients = np.polyfit(x,y,1)
        p = np.poly1d(coefficients)
        coefficient_of_dermination = r2_score(y, p(x))

        self.trendLine_content1_label.setText("Trendline: " + str(p))
        self.trendLine_content2_label.setText("R: " + str(coefficient_of_dermination))
        return coefficients

def FitPolynomial(fig, x, y, n):

#solution, res = np.polynomial.polynomial.polyfit(x,y,order,full=True)

    solution, C_p = np.polyfit(x, y, n, cov=True)  # C_z is estimated covariance matrix

    # Do the interpolation for plotting:
    xrange = fig.get_xlim()
    t = np.linspace(xrange[0],xrange[1],100)
    # Matrix with rows 1, t, t**2, ...:
    TT = np.vstack([t**(n-i) for i in range(n+1)]).T
    yi = np.dot(TT, solution)  # matrix multiplication calculates the polynomial values
    C_yi = np.dot(TT, np.dot(C_p, TT.T)) # C_y = TT*C_z*TT.T
    sig_yi = np.sqrt(np.diag(C_yi))  # Standard deviations are sqrt of diagonal


    fig.plot(t,yi, 'k-')
    # fig.plot(t,yi+sig_yi, 'k--')
    # fig.plot(t,yi-sig_yi, 'k--')

    return solution

# Take the log of the y value

def compute_xvvr(self):
        """ Return xvv(r) matrix """
        r = np.array([i*self.dr for i in range(self.ngrid)])
        k = self.get_k()
        xvvr = [["" for i in range(self.nsites)] for j in range(self.nsites)]
        for i in range(self.nsites):
            for j in range(self.nsites):
                xvvk_ij = self.xvv_data[:,i,j]
                xvvr_ij = pubfft.sinfti(xvvk_ij*k, self.dr, -1)/r
#                n_pots_for_interp = 6
#                r_for_interp = r[1:n_pots_for_interp+1]
#                xvvr_for_interp = xvvr_ij[:n_pots_for_interp]
#                poly_coefs = np.polyfit(r_for_interp, xvvr_for_interp, 3)
#                poly_f = np.poly1d(poly_coefs)
#                xvvr[i][j] = [poly_f(0)]
                xvvr[i][j] = xvvr_ij
        return r, np.swapaxes(xvvr, 0, 2)

def compute_zr(self):
        """ Return z(r) matrix """
        r = np.array([i*self.dr for i in range(self.ngrid)])
        k, zk = self.compute_zk()
        print 'computed zk',zk.shape
        zr = [["" for i in range(self.nsites)] for j in range(self.nsites)]
        for i in range(self.nsites):
            for j in range(self.nsites):
                zk_ij = zk[1:,i,j]
                zr_ij = pubfft.sinfti(zk_ij*k[1:], self.dr, -1)/r[1:]
                #zr_ij = np.abs(fftpack.fft(zk_ij))
                n_pots_for_interp = 6
                r_for_interp = r[1:n_pots_for_interp+1]
                zr_for_interp = zr_ij[:n_pots_for_interp]
                poly_coefs = np.polyfit(r_for_interp, zr_for_interp, 3)
                poly_f = np.poly1d(poly_coefs)
                zr[i][j] = [poly_f(0)]
                zr[i][j].extend(zr_ij)
        return r, np.swapaxes(zr, 0, 2)

def smooth(x, y, weights):
    '''
    in case the NLF cannot be described by 
    a square root function
    commit bounded polynomial interpolation
    '''
    # Spline hard to smooth properly, therefore solfed with
    # bounded polynomal interpolation
    # ext=3: no extrapolation, but boundary value
#     return UnivariateSpline(x, y, w=weights,
#                             s=len(y)*weights.max()*100, ext=3)

#     return np.poly1d(np.polyfit(x,y,w=weights,deg=2))
    p = np.polyfit(x, y, w=weights, deg=2)
    if np.any(np.isnan(p)):
        # couldn't even do polynomial fit
        # as last option: assume constant noise
        my = np.average(y, weights=weights)
        return lambda x: my
    return lambda xint: np.poly1d(p)(np.clip(xint, x[0], x[-1]))

def __init__(self,xpixel,ypixel,library,zoomLevel):
        self.xpixel=xpixel
        self.ypixel=ypixel
        #load csv data. Rows are zoom levels, column are latitude
        #and values are meters per pixel
        with open(library,'rb') as csvfile:
            metersperpixel = list(csv.reader(csvfile,delimiter=","))
        latitudes=[]
        #convert to floats
        for element in metersperpixel[0][1:]:
            latitudes.append(float(element))
        for row in range(1,len(metersperpixel)):
            res_values=[]
            if int(metersperpixel[row][0])==zoomLevel:
            #convert to floats
                for element in metersperpixel[row][1:]:
                    res_values.append(float(element))
                self.fitvalues=\
                    numpy.polyfit(latitudes, #fit to latitutde values
                    res_values, 3)  #3rd degree polynomial fit
                print "Fit done."   
                break

def callback_click_ok(self,data):
        points = np.array(self.interal_data)    #[(1, 1), (2, 2), (3, 3), (7, 3), (9, 3)]
        # get x and y vectors
        x = points[:,0]
        y = points[:,1]

        # calculate polynomial
        terms=int(self.sp.value())
        self.ret = np.polyfit(x, y, terms)
        f = np.poly1d(self.ret)

        tot=""
        val=0
        for i in range(0,len(self.ret)):
            p=len(self.ret)-1-i
            tot=tot+str(self.ret[i])+"*pow(w,"+str(p)+")"+"+"
        tot=tot[:-1]
        self.ret_math=tot

        self.close()

def build_model_poly(detection_pairs, beacon_sdoa, nominal_sample_rate, deg=2):
    if len(detection_pairs) < deg + 1:
        # not enough beacon transmissions
        return None

    soa0 = np.array([d[0].soa for d in detection_pairs])
    soa1 = np.array([d[1].soa for d in detection_pairs])
    soa1at0 = soa1 + np.array(beacon_sdoa)
    coef = np.polyfit(soa1at0, soa0, deg)
    fit = np.poly1d(coef)
    # residuals = soa0 - fit(soa1at0)
    # print(np.mean(residuals))

    def evaluate(det0, det1):
        return (det0.soa - fit(det1.soa)) / nominal_sample_rate

    return evaluate

def ExpLine(*args, **kwargs):
        rate = kwargs.pop('rate', None)
        log_data = kwargs.pop('log_data', True)
        data = kwargs.pop('data', None)
        const = kwargs.pop('const', 1)
        if rate is None:
            assert data is not None and len(data) > 0, "rate or data must be given"
            x = np.array([d[0] for d in data])
            y = np.array([d[1] for d in data])
            if len(x) > 0 and len(y) > 0:
                if log_data:
                    rate = np.polyfit(np.log(x), np.log(y), 1)[0]
                else:
                    rate = np.polyfit(x, y, 1)[0]
        if "label" in kwargs:
            kwargs["label"] = kwargs.pop("label").format(rate=rate)
        return FunctionLine2D(*args, fn=lambda x, r=rate:
                              const*np.array(x)**r, data=data,
                              log_data=log_data, **kwargs)

def redo_fit(self):
        lx0, lx1 = self.power_plot_lr.getRegion()

        x0, x1 = 10**lx0, 10**lx1


        X = self.dat['power_meter_power']

        n = len(X)

        ii0 = np.argmin(np.abs(X[:n//2+1]-x0))
        ii1 = np.argmin(np.abs(X[:n//2+1]-x1))

        print(ii0,ii1)

        m, b = np.polyfit(np.log10(X[ii0:ii1]), np.log10(self.power_plot_y[ii0:ii1]), deg=1)
        print("fit", m,b) 

        fit_data = 10**(np.poly1d((m,b))(np.log10(X)))
        print("fit_data", fit_data)    
        self.power_fit_plotcurve.setData(X, fit_data)

        self.fit_text.setHtml("<h1>I<sup>{:1.2f}</sup></h1>".format(m))
        self.fit_text.setPos(0.5*(lx0+lx1), np.log10(fit_data[(ii0+ii1)//2]))

def test_polyfit_build(self):
        # Ticket #628
        ref = [-1.06123820e-06, 5.70886914e-04, -1.13822012e-01,
               9.95368241e+00, -3.14526520e+02]
        x = [90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
             104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115,
             116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 129,
             130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
             146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157,
             158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
             170, 171, 172, 173, 174, 175, 176]
        y = [9.0, 3.0, 7.0, 4.0, 4.0, 8.0, 6.0, 11.0, 9.0, 8.0, 11.0, 5.0,
             6.0, 5.0, 9.0, 8.0, 6.0, 10.0, 6.0, 10.0, 7.0, 6.0, 6.0, 6.0,
             13.0, 4.0, 9.0, 11.0, 4.0, 5.0, 8.0, 5.0, 7.0, 7.0, 6.0, 12.0,
             7.0, 7.0, 9.0, 4.0, 12.0, 6.0, 6.0, 4.0, 3.0, 9.0, 8.0, 8.0,
             6.0, 7.0, 9.0, 10.0, 6.0, 8.0, 4.0, 7.0, 7.0, 10.0, 8.0, 8.0,
             6.0, 3.0, 8.0, 4.0, 5.0, 7.0, 8.0, 6.0, 6.0, 4.0, 12.0, 9.0,
             8.0, 8.0, 8.0, 6.0, 7.0, 4.0, 4.0, 5.0, 7.0]
        tested = np.polyfit(x, y, 4)
        assert_array_almost_equal(ref, tested)

def scatter_crimes_population():
    """creates a scatter plot using the values of Population and Crimes Per 100000.
    iterates through the database and reads all the data in the given headers and creates plots for each data point
    """
    x = df["Number of Crimes"].values
    y = df["Population"].values
    assert len(x) == len(y)
    df["Crimes Per 100000"] = np.array([(x[i] / y[i]) * 100000.0 for i in range(len(x))], dtype="float32")

    ax = df.plot.scatter(y="Population", x="Crimes Per 100000")

    for index, row in df.iterrows():
        ax.annotate(row["Community Name"], (row["Crimes Per 100000"], row["Population"]),
                    size=7,
                    color='darkslategrey')

    x = df["Crimes Per 100000"].values
    y = df["Population"].values

    m, b = np.polyfit(x, y, 1)
    line = plt.plot(x, m * x + b, 'b--')
    plt.setp(line, color='orange', alpha=0.5, linewidth=2.0)
    plt.show()