from mathutils import PolynomialFunction

def quadfit(xdata, ydata):
    """
    Quadratic fit
    """    
    p = (a0, a1, a2) = fit_polynomial(ydata, xdata, 2)
    f = PolynomialFunction(p)
    if a1 != 0:
        x_ext = -a1 / (2 * a0)
        y_ext = f.value(x_ext)
    else:
        x_ext = None
        y_ext = None
    yhat = [f.value(val) for val in xdata]
    ybar = sum(ydata)/len(ydata)
    ssreg = sum([(val - ybar)**2 for val in yhat])
    sstot = sum([(val - ybar)**2 for val in ydata])
    R2 = ssreg / sstot
    x1 = min(xdata)
    x2 = max(xdata)
    x_fit = frange(x1, x2, (x2-x1)/100, True)
    y_fit = [f.value(val) for val in x_fit]
    return (p, x_ext, y_ext, R2, x_fit, y_fit)

xx=(1.0,2.0,3.0,4.0)
yy=(1.1,3.9,8.5,15.8)
quadfit(xx, yy)