# mrocklin/multipolyfit

Factored out model function. Made SymPy function

1 parent 972adb4 commit a536be4ca14caec8dae25af4da14abb96048d9ed committed Jun 12, 2012
Showing with 58 additions and 17 deletions.
1. +28 −14 multipolyfit.py
2. +30 −3 test.py
@@ -57,25 +57,39 @@ def multipolyfit(xs, y, deg, full=False, model_out=False, powers_out=False):
for i in range(num_covariates+1)]
# All combinations of degrees
- powers = itertools.combinations_with_replacement(generators, deg)
+ powers = map(sum, itertools.combinations_with_replacement(generators, deg))
# Raise data to specified degree pattern, stack in order
- A = hstack(asarray([as_tall((xs**sum(p)).prod(1)) for p in powers]))
+ A = hstack(asarray([as_tall((xs**p).prod(1)) for p in powers]))
beta = linalg.lstsq(A, y)[0]
if model_out:
- # Create a function that takes in many x values
- # and returns an approximate y value
- def model(*args):
- if len(args)!=(num_covariates):
- raise ValueError("Expected %d inputs"%num_covariates)
- xs = asarray((1,) + args)
- powers = itertools.combinations_with_replacement(generators, deg)
- return sum([coeff * (xs**sum(p)).prod()
- for p, coeff in zip(powers, beta)])
- return model
+ return mk_model(beta, powers)
+
if powers_out:
- return (beta,
- list(itertools.combinations_with_replacement(generators, deg)))
+ return beta, powers
return beta
+
+def mk_model(beta, powers):
+ """ Create a callable python function out of beta/powers from multipolyfit
+
+ This function is callable from within multipolyfit using the model_out flag
+ """
+ # Create a function that takes in many x values
+ # and returns an approximate y value
+ def model(*args):
+ num_covariates = len(powers[0]) - 1
+ if len(args)!=(num_covariates):
+ raise ValueError("Expected %d inputs"%num_covariates)
+ xs = asarray((1,) + args)
+ return sum([coeff * (xs**p).prod()
+ for p, coeff in zip(powers, beta)])

#### royguo Jul 19, 2013

Is it right to use inner product for `(xs**p)` here ? I can't find any references about this. 😄

+ return model
+
+def mk_sympy_function(beta, powers):
+ from sympy import symbols, Add, Mul, S
+ num_covariates = len(powers[0]) - 1
+ xs = (S.One,) + symbols('x0:%d'%num_covariates)
+ return Add(*[coeff * Mul(*[x**deg for x, deg in zip(xs, power)])
+ for power, coeff in zip(powers, beta)])
33 test.py
 @@ -1,10 +1,37 @@ -from multipolyfit import multipolyfit +from multipolyfit import multipolyfit, mk_sympy_function import numpy as np def test_multipolyfit(): - xs = np.array([[1,2,3,4,5],[1,1,1,1,1]]).T + xs = np.array([[1,2,3,4,5]]).T y = np.array([2,3,4,5,6]) + betas = multipolyfit(xs, y, 2) + + # 1 + x + 0x**2 + assert np.linalg.norm(betas - np.asarray((1,1,0))) < .0001 + + + + +def test_model(): + xs = np.array([[1,2,3,4,5], [1,-1,1,-1,1]]).T + y = np.array([2,2,4,4,6]) model = multipolyfit(xs, y, 2, model_out = True) + yhat = np.asarray(map(model, xs[:,0], xs[:,1])) - assert np.linalg.norm(y - np.asarray(yhat)) < .0001 + assert np.linalg.norm(y - yhat) < .0001 + +def test_sympy(): + try: + from sympy import symbols + xs = np.array([[1,2,3,4,5], [1,-1,1,-1,1]]).T + y = np.array([2,3,4,5,6]) + betas, powers = multipolyfit(xs, y, 2, powers_out=True) + + expr = mk_sympy_function(betas, powers) + x0, x1 = symbols('x0,x1') + symbol_sets = {tuple(sorted(arg.free_symbols)) for arg in expr.args} + assert symbol_sets == {tuple(), (x0,), (x1,), (x0, x1)} + except ImportError: + pass +