[PYTHON] 고정 점으로 다항식을 맞추는 법

고정 점으로 다항식을 맞추는 법

파이썬에서 numpy (최소 제곱을 사용)를 사용하여 피팅을 해왔습니다.

데이터를 고정 소수점을 통해 강제로 가져 오는 방법이 있는지 궁금합니다. 그렇지 않다면 파이썬 (또는 다른 언어로 링크 할 수있는 다른 라이브러리가 있습니까? 예 : c)?

NOTE 여기서 언급했듯이 하나의 고정 소수점을 원점으로 이동시키고 상수 항을 0으로 강제하는 것이 가능하다는 것을 알고 있지만보다 일반적으로 2 개 이상의 고정 점에 대해 궁금해하고 있습니다.

http://www.physicsforums.com/showthread.php?t=523360

해결법

1. ==============================

   1.고정 된 점을 사용하여 수학적으로 올바른 방법을 사용하려면 라그랑 지 승수를 사용해야합니다. 기본적으로, 최소화하려는 목표 함수를 수정합니다. 일반적으로 잔차의 제곱의 합계이며 모든 고정 소수점에 대한 추가 매개 변수가 추가됩니다. scipy의 최소화 도구 중 하나에 수정 된 목적 함수를 제공하는 데 성공하지 못했습니다. 그러나 다항식 적합성의 경우 펜과 종이로 세부 사항을 파악하고 문제를 선형 방정식 시스템의 해로 변환 할 수 있습니다.

   고정 된 점을 사용하여 수학적으로 올바른 방법을 사용하려면 라그랑 지 승수를 사용해야합니다. 기본적으로, 최소화하려는 목표 함수를 수정합니다. 일반적으로 잔차의 제곱의 합계이며 모든 고정 소수점에 대한 추가 매개 변수가 추가됩니다. scipy의 최소화 도구 중 하나에 수정 된 목적 함수를 제공하는 데 성공하지 못했습니다. 그러나 다항식 적합성의 경우 펜과 종이로 세부 사항을 파악하고 문제를 선형 방정식 시스템의 해로 변환 할 수 있습니다.

   def polyfit_with_fixed_points(n, x, y, xf, yf) :
       mat = np.empty((n + 1 + len(xf),) * 2)
       vec = np.empty((n + 1 + len(xf),))
       x_n = x**np.arange(2 * n + 1)[:, None]
       yx_n = np.sum(x_n[:n + 1] * y, axis=1)
       x_n = np.sum(x_n, axis=1)
       idx = np.arange(n + 1) + np.arange(n + 1)[:, None]
       mat[:n + 1, :n + 1] = np.take(x_n, idx)
       xf_n = xf**np.arange(n + 1)[:, None]
       mat[:n + 1, n + 1:] = xf_n / 2
       mat[n + 1:, :n + 1] = xf_n.T
       mat[n + 1:, n + 1:] = 0
       vec[:n + 1] = yx_n
       vec[n + 1:] = yf
       params = np.linalg.solve(mat, vec)
       return params[:n + 1]
   

   그것이 작동하는지 테스트하려면 다음을 시도하십시오. 여기서 n은 점의 수, d는 다항식의 차수, f는 고정 된 점의 수입니다.

   n, d, f = 50, 8, 3
   x = np.random.rand(n)
   xf = np.random.rand(f)
   poly = np.polynomial.Polynomial(np.random.rand(d + 1))
   y = poly(x) + np.random.rand(n) - 0.5
   yf = np.random.uniform(np.min(y), np.max(y), size=(f,))
   params = polyfit_with_fixed_points(d, x , y, xf, yf)
   poly = np.polynomial.Polynomial(params)
   xx = np.linspace(0, 1, 1000)
   plt.plot(x, y, 'bo')
   plt.plot(xf, yf, 'ro')
   plt.plot(xx, poly(xx), '-')
   plt.show()
   

   물론 피팅 된 다항식은 점들을 정확히 통과합니다 :

   >>> yf
   array([ 1.03101335,  2.94879161,  2.87288739])
   >>> poly(xf)
   array([ 1.03101335,  2.94879161,  2.87288739])
   

2. ==============================

   2.curve_fit ()를 사용하면 sigma 인수를 사용하여 모든 점에 가중치를 부여 할 수 있습니다. 다음 예제는 첫 번째, 중간, 마지막 점이 매우 작은 시그마를 제공하므로 피팅 결과가이 세 점에 매우 가깝습니다.

   curve_fit ()를 사용하면 sigma 인수를 사용하여 모든 점에 가중치를 부여 할 수 있습니다. 다음 예제는 첫 번째, 중간, 마지막 점이 매우 작은 시그마를 제공하므로 피팅 결과가이 세 점에 매우 가깝습니다.

   N = 20
   x = np.linspace(0, 2, N)
   np.random.seed(1)
   noise = np.random.randn(N)*0.2
   sigma =np.ones(N)
   sigma[[0, N//2, -1]] = 0.01
   pr = (-2, 3, 0, 1)
   y = 1+3.0*x**2-2*x**3+0.3*x**4 + noise
   
   def f(x, *p):
       return np.poly1d(p)(x)
   
   p1, _ = optimize.curve_fit(f, x, y, (0, 0, 0, 0, 0), sigma=sigma)
   p2, _ = optimize.curve_fit(f, x, y, (0, 0, 0, 0, 0))
   
   x2 = np.linspace(0, 2, 100)
   y2 = np.poly1d(p)(x2)
   plot(x, y, "o")
   plot(x2, f(x2, *p1), "r", label=u"fix three points")
   plot(x2, f(x2, *p2), "b", label=u"no fix")
   legend(loc="best")
   

3. ==============================

   3.하나의 간단하고 직접적인 방법은 다음과 같이 구속 조건에 큰 M이 가중되는 제한된 최소 제곱을 이용하는 것입니다.

   하나의 간단하고 직접적인 방법은 다음과 같이 구속 조건에 큰 M이 가중되는 제한된 최소 제곱을 이용하는 것입니다.

   from numpy import dot
   from numpy.linalg import solve
   from numpy.polynomial.polynomial import Polynomial as P, polyvander as V
   
   def clsq(A, b, C, d, M= 1e5):
       """A simple constrained least squared solution of Ax= b, s.t. Cx= d,
       based on the idea of weighting constraints with a largish number M."""
       return solve(dot(A.T, A)+ M* dot(C.T, C), dot(A.T, b)+ M* dot(C.T, d))
   
   def cpf(x, y, x_c, y_c, n, M= 1e5):
       """Constrained polynomial fit based on clsq solution."""
       return P(clsq(V(x, n), y, V(x_c, n), y_c, M))
   

   분명히 이것은 모든 포괄적 인은 총알 해결책은 아니지만 분명히 ([0, 4, 24, 124, 624, 3124]의 M에 대한) 간단한 예제를 통해 합리적으로 잘 작동하는 것처럼 보입니다.

   In []: x= linspace(-6, 6, 23)
   In []: y= sin(x)+ 4e-1* rand(len(x))- 2e-1
   In []: x_f, y_f= linspace(-(3./ 2)* pi, (3./ 2)* pi, 4), array([1, -1, 1, -1])
   In []: n, x_s= 5, linspace(-6, 6, 123)    
   
   In []: plot(x, y, 'bo', x_f, y_f, 'bs', x_s, sin(x_s), 'b--')
   Out[]: <snip>
   
   In []: for M in 5** (arange(6))- 1:
      ....:     plot(x_s, cpf(x, y, x_f, y_f, n, M)(x_s))
      ....: 
   Out[]: <snip>
   
   In []: ylim([-1.5, 1.5])
   Out[]: <snip>
   In []: show()
   

   다음과 같이 출력을 생성합니다.

   편집 : 추가 된 '정확한'솔루션 :

   from numpy import dot
   from numpy.linalg import solve
   from numpy.polynomial.polynomial import Polynomial as P, polyvander as V
   from scipy.linalg import qr 
   
   def solve_ns(A, b): return solve(dot(A.T, A), dot(A.T, b))
   
   def clsq(A, b, C, d):
       """An 'exact' constrained least squared solution of Ax= b, s.t. Cx= d"""
       p= C.shape[0]
       Q, R= qr(C.T)
       xr, AQ= solve(R[:p].T, d), dot(A, Q)
       xaq= solve_ns(AQ[:, p:], b- dot(AQ[:, :p], xr))
       return dot(Q[:, :p], xr)+ dot(Q[:, p:], xaq)
   
   def cpf(x, y, x_c, y_c, n):
       """Constrained polynomial fit based on clsq solution."""
       return P(clsq(V(x, n), y, V(x_c, n), y_c))
   

   적합성 테스트 :

   In []: x= linspace(-6, 6, 23)
   In []: y= sin(x)+ 4e-1* rand(len(x))- 2e-1
   In []: x_f, y_f= linspace(-(3./ 2)* pi, (3./ 2)* pi, 4), array([1, -1, 1, -1])
   In []: n, x_s= 5, linspace(-6, 6, 123)
   In []: p= cpf(x, y, x_f, y_f, n)
   In []: p(x_f)
   Out[]: array([ 1., -1.,  1., -1.])
   

from https://stackoverflow.com/questions/15191088/how-to-do-a-polynomial-fit-with-fixed-points by cc-by-sa and MIT license