복붙노트

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

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.])
    
  4. from https://stackoverflow.com/questions/15191088/how-to-do-a-polynomial-fit-with-fixed-points by cc-by-sa and MIT license