def weighted_linear_fit(x, y, w):
    """
    return linear fit by calcualating
    y_fit = a + b*x
    page 104 Data reduction and error analysis for the physicxal sciences Philip R. Bevington

    Parameters
    ----------
    x (1d numpy array)
    y (1d numpy array)
    w (1d numpy array)

    Returns
    -------
    a
    b
    sigma_a
    sigma_b

    Examples
    --------
    >>> a, b , sigma_a, sigma_v = weighted_linear_fit(x,y,w)
    """
    from numpy import isnan, nan, sum, sqrt
    Sx = sum(w ** 2 * x * 1.0)
    Sx2 = sum((w * x) ** 2.0)
    Sy = sum(w ** 2 * y * 1.0)
    Sy2 = sum(w ** 2 * y ** 2.0)
    Sxy = sum(w ** 2 * x * y * 1.0)
    Sw2 = sum(w ** 2)
    N = x.shape[0]
    if N >= 2:
        Delta = Sw2 * Sx2 - Sx ** 2
        a = 1.0 / Delta * (Sx2 * Sy - Sx * Sxy)
        b = 1.0 / Delta * (Sw2 * Sxy - Sx * Sy)
    else:
        a = None
        b = None
    if N > 2:
        sigma_a = sqrt(1 / Delta * Sx2)
        sigma_b = sqrt(1 / Delta * Sw2)
    else:
        sigma_a = sigma_b = None
    print('Delta = ', Delta)
    print('Sw2 = ', Sw2)
    print('Sx = ', Sx)
    print('Sx2 = ', Sx2)
    return a, b, sigma_a, sigma_b
