다항식 보간
다항식 보간 (Polynomial Interpolation)
데이터 점들을 하나의 다항식이나 조각별 다항식(스플라인)으로 맞추고 싶을 때 쓰는 함수들을 소개할게요. 최소제곱 적합부터 조각별 다항식 구조까지 하나씩 다뤄볼게요.
출처: 문서
본문
Octave는 다양한 종류의 보간을 잘 지원해요. 대부분은 보간 항목에서 설명하고 있고요. 앞선 장에서 설명한 함수들에 대한 한 가지 간단한 대안은, 주어진 데이터 점에 하나의 다항식 또는 조각별 다항식(스플라인)을 맞추는 것이에요. 심하게 요동치는 다항식을 피하려면 보통 낮은 차수의 다항식을 데이터에 맞추는 걸 원해요. 이는 대개 최소제곱(least-squares) 의미로 다항식을 맞춰야 함을 뜻하는데, 바로 그 일을 하는 함수가 polyfit이에요.
polyfit
p = polyfit (x, y, n)
[p, s] = polyfit (x, y, n)
[p, s, mu] = polyfit (x, y, n)
점 [x(:), y(:)]에 대한 적합의 최소제곱 오차를 최소화하는 n차 다항식 p(x)의 계수를 반환해요.
n은 보통 근사 다항식의 차수를 지정하는 정수(≥ 0)예요. n이 논리 벡터라면 마스크로 사용되어 해당 다항식 계수를 선택적으로 사용하거나 무시하게 해요.
다항식 계수는 행 벡터 p로 반환돼요. 출력 p는 polyval과 함께 직접 사용해 맞춘 다항식으로 값을 추정할 수 있어요.
선택 출력 s는 다음 필드를 담은 구조체예요.
'yf'— 각 x 값에 대한 다항식의 값.'V'— 다항식 계수를 계산하는 데 사용된 반데르몽드(Vandermonde) 행렬.'X (deprecated, will be removed in Octave 13)'— 다항식 계수를 계산하는 데 사용된 반데르몽드 행렬.'R'— QR 분해에서 얻은 삼각 인자 R.'C'— (확장되지 않은) 공분산 행렬. 형식상 v'*v의 역과 같지만 반올림 오차 전파를 최소화하는 방식으로 계산돼요.'df'— 자유도.'normr'— 잔차의 노름.'rsquared'— 결정 계수(R-squared).
두 번째 출력은 polyval이 예측 값의 통계적 오차 한계를 계산하는 데 사용할 수 있어요. 특히 p 계수의 표준편차는 sqrt (diag (s.C)/s.df) * s.normr로 주어져요.
세 번째 출력 mu가 있으면 원본 데이터를 중심화·스케일링해서 적합의 수치적 안정성을 높일 수 있어요. 계수 p는 다음 다항식에 연관돼요.
xhat = (x - mu(1)) / mu(2)
여기서 mu(1) = mean (x), mu(2) = std (x)예요.
예시 1: 논리 n과 정수 n
f = @(x) x.^2 + 5; # data-generating function
x = 0:5;
y = f (x);
## Fit data to polynomial A*x^3 + B*x^1
p = polyfit (x, y, logical ([1, 0, 1, 0]))
⇒ p = [ 0.0680, 0, 4.2444, 0 ]
## Fit data to polynomial using all terms up to x^3
p = polyfit (x, y, 3)
⇒ p = [ -4.9608e-17, 1.0000e+00, -1.6906e-15, 5.0000e+00 ]
프로그래밍 참고: 원하는 다항식 차수 n ≥ 데이터 점 개수라면 해가 유일하지 않아요. 대신 적합은 m = numel (x) - 1개의 항을 계산해요. 반환된 다항식은 x^n, x^n-1, …, x^n-m과 상수항 x^0의 계수를 가지며, 나머지 계수는 0으로 설정돼요.
See also: polyval, polyaffine, roots, vander, zscore.
splinefit
하나의 다항식으로 부족한 상황에서는 여러 다항식을 이어 붙이는 해법이 있어요. splinefit 함수는 조각별 다항식(스플라인)을 데이터 집합에 맞춰요.
pp = splinefit (x, y, breaks)
pp = splinefit (x, y, p)
pp = splinefit (…, "periodic", periodic)
pp = splinefit (…, "robust", robust)
pp = splinefit (…, "beta", beta)
pp = splinefit (…, "order", order)
pp = splinefit (…, "constraints", constraints)
잡음이 있는 데이터 x와 y에 끊기는 지점(매듭, knots) breaks를 가진 조각별 3차 스플라인을 맞춰요.
x는 벡터이고 y는 벡터 또는 N차원 배열이에요. y가 N차원 배열이면 x(j)는 y(:,…,:,j)와 대응돼요.
p는 x를 따라 구간(interval)의 개수를 정의하는 양의 정수이고, p+1이 끊김 개수예요. 각 구간의 점 개수는 1 이하로 차이납니다.
선택 속성 periodic은 스플라인에 주기 경계 조건을 적용할지 지정하는 논리값이에요. 주기의 길이는 max (breaks) - min (breaks)예요. 기본값은 false예요.
선택 속성 robust는 이상치(outlying) 데이터 점의 영향을 줄이는 강건(robust) 적합을 적용할지 지정하는 논리값이에요. 가중 최소제곱을 세 번 반복 수행해요. 가중치는 이전 잔차에서 계산돼요. 이상치 식별의 민감도는 속성 beta로 제어돼요. beta의 값은 0 < beta < 1 범위로 제한돼요. 기본값은 beta = 1/2예요. 0에 가까운 값은 모든 데이터에 동등한 가중치를 주고, beta 값이 증가할수록 이상치 데이터의 영향이 줄어들어요. 1에 가까운 값은 불안정성이나 계수 결핍(rank deficiency)을 일으킬 수 있어요.
맞춘 스플라인은 조각별 다항식 pp로 반환되며 ppval로 평가할 수 있어요.
스플라인은 차수 order의 다항식으로 구성돼요. 기본값은 3차, order=3이에요. P개 조각을 가진 스플라인은 P+order 개의 자유도를 가져요. 주기 경계 조건에서는 자유도가 P로 줄어들어요.
선택 속성 constraints는 적합에 대한 선형 제약을 지정하는 구조체예요. 구조체는 "xc", "yc", "cc" 세 필드를 가져요.
"xc"— 제약의 x 위치 벡터."yc"— 위치 xc에서의 제약 값. 기본값은 0 배열."cc"— 계수(행렬). 기본값은 1 배열. 행 수는 조각별 다항식의 차수order로 제한돼요.
제약은 0차부터 order-1차까지의 도함수의 선형 결합이며,
cc(1,j) * y(xc(j)) + cc(2,j) * y'(xc(j)) + ... = yc(:,...,:,j).
예요.
See also: interp1, unmkpp, ppval, spline, pchip, ppder, ppint, ppjumps.
조각별 다항식을 만드는 데 사용되는 끊김(또는 매듭) 수는 입력 데이터 x, y에 존재하는 잡음을 억제하는 데 중요한 요소예요. 아래 예시가 이를 보여줘요.
x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.2 * randn (size (x));
## Uniform breaks
breaks = linspace (0, 2 * pi, 41); % 41 breaks, 40 pieces
pp1 = splinefit (x, y, breaks);
## Breaks interpolated from data
pp2 = splinefit (x, y, 10); % 11 breaks, 10 pieces
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "41 breaks, 40 pieces", "11 breaks, 10 pieces"})
그 결과는 Figure 28.1에서 볼 수 있어요. 끊김이 많은 적합에서는 밑바탕 함수에는 없는 빠른 리플(ripple)이 나타나요.
splinefit가 제공하는 조각별 다항식 적합은 order-1차까지 연속 도함수를 가져요. 예를 들어 3차 적합은 연속 1·2차 도함수를 가져요. 다음 코드가 이를 보여줍니다.
## Data (200 points)
x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.1 * randn (size (x));
## Piecewise constant
pp1 = splinefit (x, y, 8, "order", 0);
## Piecewise linear
pp2 = splinefit (x, y, 8, "order", 1);
## Piecewise quadratic
pp3 = splinefit (x, y, 8, "order", 2);
## Piecewise cubic
pp4 = splinefit (x, y, 8, "order", 3);
## Piecewise quartic
pp5 = splinefit (x, y, 8, "order", 4);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
y3 = ppval (pp3, xx);
y4 = ppval (pp4, xx);
y5 = ppval (pp5, xx);
plot (x, y, ".", xx, [y1; y2; y3; y4; y5])
axis tight
ylim auto
legend ({"data", "order 0", "order 1", "order 2", "order 3", "order 4"})
따라서 Figure 28.2에서 볼 수 있듯 고차 해가 밑바탕 함수를 더 정확히 나타내지만 계산 복잡도가 높아져요.
적합하려는 밑바탕 함수가 주기적일 때, splinefit은 주기적 적합을 나타내는 데 필요한 경계 조건을 적용할 수 있어요.
## Data (100 points)
x = 2 * pi * [0, (rand (1, 98)), 1];
y = sin (x) - cos (2 * x) + 0.2 * randn (size (x));
## No constraints
pp1 = splinefit (x, y, 10, "order", 5);
## Periodic boundaries
pp2 = splinefit (x, y, 10, "order", 5, "periodic", true);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "no constraints", "periodic"})
더 복잡한 제약도 추가할 수 있어요. 예를 들어 아래 코드는 끝점에서 값이 고정(clamped)된 주기적 적합과, 끝점에서 힌지(hinged)된 두 번째 주기적 적합을 보여줘요.
## Data (200 points)
x = 2 * pi * rand (1, 200);
y = sin (2 * x) + 0.1 * randn (size (x));
## Breaks
breaks = linspace (0, 2 * pi, 10);
## Clamped endpoints, y = y' = 0
xc = [0, 0, 2*pi, 2*pi];
cc = [(eye (2)), (eye (2))];
con = struct ("xc", xc, "cc", cc);
pp1 = splinefit (x, y, breaks, "constraints", con);
## Hinged periodic endpoints, y = 0
con = struct ("xc", 0);
pp2 = splinefit (x, y, breaks, "constraints", con, "periodic", true);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "clamped", "hinged periodic"})
splinefit 함수는 강건 적합의 편의 기능도 제공하는데, 이상치 데이터의 영향을 줄여줘요. 아래 예시에서는 세 가지 다른 적합을 보여줘요. 둘은 서로 다른 수준의 이상치 억제를, 셋째는 강건하지 않은 해를 보여줍니다.
## Data
x = linspace (0, 2*pi, 200);
y = sin (x) + sin (2 * x) + 0.05 * randn (size (x));
## Add outliers
x = [x, linspace(0,2*pi,60)];
y = [y, -ones(1,60)];
## Fit splines with hinged conditions
con = struct ("xc", [0, 2*pi]);
## Robust fitting, beta = 0.25
pp1 = splinefit (x, y, 8, "constraints", con, "beta", 0.25);
## Robust fitting, beta = 0.75
pp2 = splinefit (x, y, 8, "constraints", con, "beta", 0.75);
## No robust fitting
pp3 = splinefit (x, y, 8, "constraints", con);
## Plot
xx = linspace (0, 2*pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
y3 = ppval (pp3, xx);
plot (x, y, ".", xx, [y1; y2; y3])
legend ({"data with outliers","robust, beta = 0.25", ...
"robust, beta = 0.75", "no robust fitting"})
axis tight
ylim auto
padecoef
아주 특수한 형태의 다항식 근사가 바로 파데 근사(Padé approximant)예요. 제어 시스템에서 연속 시간 지연을 이 근사로 매우 간단히 모델링할 수 있어요.
[num, den] = padecoef (T)
[num, den] = padecoef (T, N)
연속 시간 지연 T의 N차 파데 근사를 전달 함수 형태로 계산해요. exp (-sT)의 파데 근사는 다음 방정식으로 정의돼요.
Pn(s)
exp (-sT) ~ -------
Qn(s)
여기서 Pn(s)와 Qn(s)는 모두 N차 유리 함수이며 다음 식으로 정의돼요.
N (2N - k)!N! k
Pn(s) = SUM --------------- (-sT)
k=0 (2N)!k!(N - k)!
Qn(s) = Pn(-s)
입력 T와 N은 음수가 아닌 수치 스칼라여야 해요. N이 지정되지 않으면 기본값은 1이에요.
출력 행 벡터 num과 den은 s의 내림차수 거듭제곱으로 분자·분모 계수를 담아요. 둘 다 N차 다항식이에요.
예를 들어,
t = 0.1;
n = 4;
[num, den] = padecoef (t, n)
⇒ num =
1.0000e-04 -2.0000e-02 1.8000e+00 -8.4000e+01 1.6800e+03
⇒ den =
1.0000e-04 2.0000e-02 1.8000e+00 8.4000e+01 1.6800e+03
조각별 다항식 함수들
ppval 함수는 mkpp나 다른 방법으로 만든 조각별 다항식을 평가하고, unmkpp는 조각별 다항식에 대한 상세 정보를 반환해요.
다음 예시는 두 개의 선형 함수와 하나의 2차 함수를 하나의 함수로 결합하는 방법을 보여줘요. 각 함수는 인접한 구간에서 표현돼요.
x = [-2, -1, 1, 2];
p = [ 0, 1, 0;
1, -2, 1;
0, -1, 1 ];
pp = mkpp (x, p);
xi = linspace (-2, 2, 50);
yi = ppval (pp, xi);
plot (xi, yi);
mkpp
pp = mkpp (breaks, coefs)
pp = mkpp (breaks, coefs, d)
표본 점 breaks와 계수 coefs에서 조각별 다항식(pp) 구조체를 구성해요.
breaks는 엄격히 증가하는 값의 벡터여야 해요. 구간 수는 ni = length (breaks) - 1로 주어져요.
m이 다항식 차수일 때 coefs의 크기는 ni-by-(m + 1)이어야 해요.
coefs의 i번째 행 coefs(i,:)는 i번째 구간의 다항식 계수를 담으며, 최고차(m)에서 최저차(0) 순서로 정렬돼요.
coefs는 벡터 값 또는 배열 값 다항식을 지정하는 다차원 배열일 수도 있어요. 그 경우 다항식 차수 m은 coefs의 마지막 차원 길이로 정의돼요. 첫 번째 차원의 크기는 스칼라 또는 벡터 d로 주어져요. d가 주어지지 않으면 1로 설정돼요. 이 경우 p(r, i, :)는 구간 i에 정의된 r번째 다항식의 계수를 담아요. 어느 경우든 coefs는 크기 [ni*prod(d) m]의 2-D 행렬로 재구성돼요.
프로그래밍 참고: ppval은 xi - breaks(i), 즉 현재 구간의 하한 끝점을 xi에서 뺀 값에서 다항식을 평가해요. mkpp로 조각별 다항식 객체를 만들 때 이 점을 반드시 고려해야 해요.
See also: unmkpp, ppval, spline, pchip, ppder, ppint, ppjumps.
unmkpp
[x, p, n, k, d] = unmkpp (pp)
조각별 다항식 구조체 pp의 구성 요소를 추출해요. 이 함수는 mkpp의 역으로, pp를 만드는 데 필요한 mkpp의 입력을 추출해요. 아래 코드가 이 관계를 명확히 보여줘요.
[breaks, coefs, numinter, order, dim] = unmkpp (pp);
pp2 = mkpp (breaks, coefs, dim);
이렇게 얻은 조각별 다항식 구조체 pp2는 원래 pp와 동일해요. 구조체 pp의 필드에 직접 접근해서도 같은 결과를 얻을 수 있어요.
구성 요소는 다음과 같아요.
x— 표본 점 또는 끊김(breaks).p— 표본 구간의 점들에 대한 다항식 계수.p(i, :)는 구간 i의 다항식 계수를 최고차에서 최저차 순서로 담아요.d > 1이면 p는 크기[n*prod(d) m]의 행렬이고,i + (1:d)행은 구간 i의 모든 d개 다항식의 계수예요.n— 다항식 조각 또는 구간 수,n = length (x) - 1.k— 다항식의 차수 + 1.d— 각 구간에 정의된 다항식 수.
See also: mkpp, ppval, spline, pchip.
ppval
yi = ppval (pp, xi)
점 xi에서 조각별 다항식 구조체 pp를 평가해요.
pp가 스칼라 다항식 함수를 설명하면 결과는 xi와 같은 형태의 배열이에요. 그렇지 않으면 결과의 크기는 xi가 벡터일 때 [pp.dim, length(xi)], 다차원 배열일 때 [pp.dim, size(xi)]예요.
See also: mkpp, unmkpp, spline, pchip.
ppder
ppd = ppder (pp)
ppd = ppder (pp, m)
조각별 다항식 구조체 pp의 m차 도함수를 계산해요. m이 생략되면 1차 도함수를 계산해요.
See also: mkpp, ppval, ppint.
ppint
ppi = ppint (pp)
ppi = ppint (pp, c)
조각별 다항식 구조체 pp의 적분을 계산해요. c가 주어지면 적분 상수예요.
See also: mkpp, ppval, ppder.
ppjumps
jumps = ppjumps (pp)
조각별 다항식의 경계 점프(boundary jumps)를 평가해요. 구간이 n개이고 pp의 차원이 d라면 결과 배열의 크기는 [d, n-1]이에요.
See also: mkpp.
더 알아보기
- 일반 보간 함수는 보간 항목을 참고하세요.
- 다항식 평가와 근은 다항식 처리 항목도 함께 보면 좋아요.