Error of value of polynomial with covariance matrix for coefficents

Question

From polyfit() function I have coefficents and their covariance matrix. Now for some argument I want a value of that polynomial with error. Apparently polyval() function doesn't accept covariance matrix (in matlab it does). From definition I can write my own function:

def sigma(x,pcov):
    s=0
    for i in range(len(pcov)):
        for j in range(len(pcov[0])):
            s+=x**(len(pcov)-i-1)*pcov[i][j]*x**(len(pcov[0])-j-1)
    return sqrt(s)

Where I have to do some funny stuff with indices because of reverse order of coefficents in polynomial. It doesn't seem pythonic though. Am I missing something? Or maybe support for this kind of operation is in some bigger library like SciPy?

Answer 1

Without loops, this is computed by multiplying the covariance matrix C on both sides by the vector x_powers made of the relevant powers of x: for example, [x**5, x**4, ..., x**0]. The setup of polyfit is included for completeness.

xdata = np.arange(10)
ydata = np.abs(xdata-2)    # some data to fit
degree = 5
p, C = np.polyfit(xdata, ydata, deg=degree, cov=True)
x = 5.4                    # some value of x to estimate error at
x_powers = x**np.arange(degree, -1, -1)   
x_error = np.sqrt(x_powers.dot(C).dot(x_powers))

Here x_error is the same as what your function returns.