Obtaining Legendre polynomial form once Legendre coefficients are determined

To simplify Ahmed's example\n\nIn [1]: from numpy.polynomial import Polynomial, Legendre\n\nIn [2]: p = Polynomial([0.5, 0.3, 0.1])\n\nIn [3]: x = np.random.rand(10) * 10\n\nIn [4]: y = p(x)\n\nIn [5]: pfit = Legendre.fit(x, y, 2)\n\nIn [6]: plot(*pfit.linspace())\nOut[6]: [<matplotlib.lines.Line2D at 0x7f815364f310>]\n\nIn [7]: plot(x, y, 'o')\nOut[7]: [<matplotlib.lines.Line2D at 0x7f81535d8bd0>]\n\n\nThe Legendre functions are scaled and offset, as the data should be confined to the interval [-1, 1] to get any advantage over the usual power basis. If you want the coefficients for plain old Legendre functions\n\nIn [8]: pfit.convert()\nOut[8]: Legendre([ 0.53333333, 0.3 , 0.06666667], [-1., 1.], [-1., 1.])\n\n\nBut that isn't recommended.",

Once you have a function, you can just generate a numpy array for the timepoints:\n\n>>> import numpy as np\n>>> timepoints = [1,3,7,15,16,17,19]\n>>> myarray = np.array(timepoints)\n>>> def mypolynomial(bins, pfinal): #pfinal is just the estimate of the final array (i'll do quadratic)\n... a,b,c = pfinal # obviously, for a*x^2 + b*x + c\n... return (a*bins**2) + b*bins + c\n>>> mypolynomial(myarray, (1,1,0))\narray([ 2, 12, 56, 240, 272, 306, 380])\n\n\nIt automatically evaluates it for each timepoint is in the numpy array.\n\nNow all you have to do is rewrite mypolynomial to go from a simple quadratic example to a proper one for a Legendre polynomial. Treat the function as if it were evaluating a float to return the value, and when called on the numpy array it will automatically evaluate it for each value.\n\nEDIT:\nLet's say I wanted to generalize this to all standard polynomials:\n\n>>> import numpy as np\n>>> timepoints = [1,3,7,15,16,17,19]\n>>> myarray = np.array(timepoints)\n>>> def mypolynomial(bins, pfinal): #pfinal is just the estimate of the final array (i'll do quadratic)\n>>> hist = np.zeros((1, len(myarray))) # define blank return\n... for i in range(len(pfinal)):\n... # fixed a typo here, was pfinal[-i] which would give -0 rather than -1, since negative indexing starts at -1, not -0\n... const = pfinal[-i-1] # negative index to go from 0 exponent to highest exponent\n... hist += const*(bins**i)\n... return hist\n>>> mypolynomial(myarray, (1,1,0))\narray([ 2, 12, 56, 240, 272, 306, 380])\n\n\nEDIT2: Typo fix\n\nEDIT3:\n\n@Ahmed is perfect right when he states Homer's rule is good for numerical stability. The implementation here would be as follows:\n\n>>> def horner(coeffs, x):\n... acc = 0\n... for c in coeffs:\n... acc = acc * x + c\n... return acc\n>>> horner((1,1,0), myarray)\narray([ 2, 12, 56, 240, 272, 306, 380])\n\n\nSlightly modified to keep the same argument order as before, from the code here:\nhttp://rosettacode.org/wiki/Horner%27s_rule_for_polynomial_evaluation#Python",