반응표면의 1 차·2 차 모형 — Steepest Ascent 와 정상점 분석

Montgomery Ch.7.7.1-7.7.4 Response Surface · Linear · Second Order · Variance

Response Surface Methodology (RSM) 의 두 핵심 모형 — 1 차 (선형 plane) 와 2 차 (곡면) — 의 정의, 적합 절차, 응답 분산 구조를 정리한다. 작동 영역에서 멀리 있는 단계의 1 차 모형 + steepest ascent 를 통한 최적 영역 탐색, 최적 영역에서의 2 차 모형 + 정상점 분석 (canonical form, eigenvalue 분류) 의 통합 framework 를 다룬다.

Experimentation
DOE
저자

Kwangmin Kim

공개

2026년 05월 08일

1 정의와 motivation

정의: Response Surface Methodology (RSM)

여러 양적 요인 \(\mathbf{x} = (x_1, x_2, \ldots, x_k)^T\) 가 응답 변수 \(Y\) 에 미치는 효과를 곡면 \(\eta(\mathbf{x}) = E[Y | \mathbf{x}]\) 로 모형화하고, 그 곡면의 최적 작동점 \(\mathbf{x}^*\) 를 추정하는 통계·실험 framework.

핵심 단계: 1. 1 차 모형 적합 (현재 영역이 최적에서 멀 때): linear plane. 2. Steepest ascent: 응답이 가장 빠르게 증가하는 방향으로 실험 이동. 3. 2 차 모형 적합 (최적 영역 진입 시): quadratic surface. 4. 정상점 추정: \(\nabla \eta = 0\) 의 해. 5. Canonical analysis: 정상점 종류 (max/min/saddle) 분류.

2 RSM 의 단계적 전략

직관: 산을 오르는 통계학자의 비유

당신이 안개 속의 산을 오르고 있다. 정상이 어디 있는지 모른다.

  1. 현재 위치에서 작은 영역을 평가 (1 차 모형 fit) → 가장 가파른 경사 방향 (steepest ascent) 확인.
  2. 그 방향으로 이동 (steepest ascent path 의 일정 step). 새 위치에서 1 차 fit 재시도.
  3. 1 차의 적합도가 떨어지기 시작 (선형성 검정의 lack of fit 유의) → 곡률이 보임 → 정상 근처.
  4. 정상 근처에서 2 차 모형 fit (지형의 곡면).
  5. 정상점 (\(\nabla \eta = 0\)) 추정 → 최적 작동점.

이것이 RSM 의 핵심 전략. 작동 영역 (region of operability) 을 점차 좁혀 최적점을 찾는다.

3 1 차 모형 (First-Order Model)

3.1 모형

\[ \eta(\mathbf{x}) = \beta_0 + \sum_{i=1}^{k} \beta_i x_i \]

\(k\) 변수, \(k + 1\) 모수. 가장 간단한 형태.

3.2 적합

표준 다중 회귀: \[ \hat{\boldsymbol{\beta}} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y} \]

설계로 \(2^k\) factorial 또는 fractional factorial 사용. \(\mathbf{X}^T \mathbf{X}\) 가 직교라 \(\hat\beta_i\) 의 분산 분리.

3.3 Steepest Ascent

응답을 가장 빠르게 증가시키는 방향 = \((\beta_1, \beta_2, \ldots, \beta_k)\) 의 unit vector.

각 \(x_i\) 의 변화량은 \(\beta_i\) 에 비례:

\[ \frac{\Delta x_1}{\hat\beta_1} = \frac{\Delta x_2}{\hat\beta_2} = \ldots = \frac{\Delta x_k}{\hat\beta_k} \]

가장 큰 \(|\hat\beta_i|\) 변수의 한 단위 변화에 비례.

3.4 단계적 절차 — Steepest Ascent

  1. 기준 step size 결정: 가장 영향이 큰 변수 (\(|\hat\beta_{\max}|\)) 의 1 단위 변화.
  2. 다른 변수의 변화량 계산: 비례식.
  3. 새 점에서 응답 측정: 단일 또는 소수 점.
  4. 응답 증가 추세 모니터링: 증가하는 한 step 더 진행. 감소·정체 시 중지.
  5. 새 영역에서 1 차 fit 재시도: 적합도 점검 후 단계 1 또는 단계 7 (2 차).

3.5 Lack-of-Fit Test

1 차 모형이 데이터에 적합한지 검정. center point 추가가 핵심.

Center point: 모든 \(x_i = 0\) 의 점에서 여러 회 측정. 이 점들은 1 차 모형에서는 평균 응답이 \(\beta_0\), 곡률이 있으면 \(\beta_0\) 이 아니다.

\[ SS_{\text{curvature}} = \frac{n_C n_F (\bar y_F - \bar y_C)^2}{n_C + n_F} \]

\(n_F\) = factorial 점 수, \(n_C\) = center 점 수, \(\bar y_F\), \(\bar y_C\) 는 각 평균.

자유도 1. \(F\) 검정이 유의면 곡률이 있다 → 2 차 모형으로 전환.

직관: Center Point 의 통계적 역할

Factorial 점은 영역 가장자리에서 평균 응답을 추정. Center point 는 한가운데에서 추정.

만약 진짜 응답이 선형이면 center 점의 평균 = factorial 점들의 평균 (곡면이 평면이라 중간이 평균과 같음).

다르면 곡률이 있다. 차이의 크기가 곡률의 크기 → SS_curvature.

center 점은 또 분산 추정도 제공 (반복으로 \(\sigma^2\) 추정).

4 2 차 모형 (Second-Order Model)

4.1 모형

\[ \eta(\mathbf{x}) = \beta_0 + \sum_{i=1}^k \beta_i x_i + \sum_{i=1}^k \beta_{ii} x_i^2 + \sum_{i < j}^k \beta_{ij} x_i x_j \]

모수 수: \(1 + k + k + \binom{k}{2} = (k+1)(k+2)/2\).

\(k\) 모수 수
2 6
3 10
4 15
5 21
6 28

4.2 적합

다중 회귀. Design 으로 Central Composite Design (G-MON7-4) 또는 Box-Behnken 사용.

4.3 Matrix Form

\[ \eta(\mathbf{x}) = \beta_0 + \mathbf{x}^T \mathbf{b} + \mathbf{x}^T \mathbf{B} \mathbf{x} \]

\(\mathbf{b} = (\beta_1, \beta_2, \ldots, \beta_k)^T\) — 1 차 항 vector.

\(\mathbf{B}\) — 대칭 \(k \times k\) 행렬: \[ \mathbf{B} = \begin{pmatrix} \beta_{11} & \beta_{12}/2 & \cdots & \beta_{1k}/2 \\ \beta_{12}/2 & \beta_{22} & \cdots & \beta_{2k}/2 \\ \vdots & \vdots & \ddots & \vdots \\ \beta_{1k}/2 & \beta_{2k}/2 & \cdots & \beta_{kk} \end{pmatrix} \]

(off-diagonal 이 \(/2\) — interaction 모수가 두 곳에서 contribution.)

5 정상점 (Stationary Point)

응답의 그라디언트가 0 인 점: \[ \nabla \eta(\mathbf{x}) = \mathbf{b} + 2 \mathbf{B} \mathbf{x} = 0 \]

해: \[ \boxed{\mathbf{x}^* = -\frac{1}{2} \mathbf{B}^{-1} \mathbf{b}} \]

정상점에서의 응답: \[ \hat\eta(\mathbf{x}^*) = \hat\beta_0 + \frac{1}{2} \mathbf{x}^{*T} \mathbf{b} \]

5.1 정상점의 분류 — Canonical Analysis

\(\mathbf{B}\) 의 eigenvalues \(\lambda_1, \ldots, \lambda_k\) 의 부호로 분류:

Eigenvalue 부호 정상점 종류
모두 음 (\(\lambda_i < 0\)) Maximum (응답 정점)
모두 양 (\(\lambda_i > 0\)) Minimum
부호 혼재 Saddle point
일부 ≈ 0 Ridge (평탄한 능선)

5.2 Canonical Form

좌표 변환 \(\mathbf{w} = \mathbf{P}^T (\mathbf{x} - \mathbf{x}^*)\) (\(\mathbf{P}\) = \(\mathbf{B}\) 의 eigenvector matrix):

\[ \hat\eta = \hat\eta(\mathbf{x}^*) + \sum_{i=1}^k \lambda_i w_i^2 \]

각 \(w_i\) 방향의 곡률이 \(\lambda_i\). 큰 \(|\lambda_i|\) 방향이 더 가파름.

5.3 사례 — 2 차원 RSM

\(k = 2\), fitted model: \[ \hat\eta = 80 + 1.5 x_1 + 2 x_2 - 1.5 x_1^2 - 2 x_2^2 + 0.5 x_1 x_2 \]

\(\mathbf{b} = (1.5, 2)^T\), \(\mathbf{B} = \begin{pmatrix} -1.5 & 0.25 \\ 0.25 & -2 \end{pmatrix}\).

\(\mathbf{B}^{-1} = \frac{1}{(-1.5)(-2) - 0.25^2} \begin{pmatrix} -2 & -0.25 \\ -0.25 & -1.5 \end{pmatrix} = \frac{1}{2.9375} \begin{pmatrix} -2 & -0.25 \\ -0.25 & -1.5 \end{pmatrix}\)

\(\mathbf{x}^* = -\frac{1}{2} \mathbf{B}^{-1} \mathbf{b}\) \(= -\frac{1}{2} \frac{1}{2.9375} (-2, -0.25) (1.5, 2)\) \(= \frac{1}{5.875} (3 + 0.5, 0.375 + 3)\) \(= (3.5/5.875, 3.375/5.875)\) \(= (0.596, 0.574)\)

응답값: \(\hat\eta(\mathbf{x}^*) = 80 + 1.5 \times 0.596 + 2 \times 0.574 - 1.5 \times 0.596^2 - 2 \times 0.574^2 + 0.5 \times 0.596 \times 0.574\) \(= 80 + 0.894 + 1.148 - 0.532 - 0.659 + 0.171\) \(= 81.02\)

Eigenvalue 검사: \(\mathbf{B}\) 의 eigenvalue = \(-1.5, -2\) 근방 (off-diagonal 작아 거의 대각). 둘 다 음 → Maximum.

→ 최적 작동점 \(\mathbf{x}^* \approx (0.6, 0.57)\), 최대 응답 81.

6 응답 분산 — 예측 정밀도

설계 점에서 멀어질수록 예측 분산 ↑:

\[ \text{Var}(\hat Y(\mathbf{x})) = \sigma^2 \mathbf{x}_*^T (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{x}_* \]

여기서 \(\mathbf{x}_* = (1, x_1, \ldots, x_k, x_1^2, \ldots, x_1 x_2, \ldots)\) 는 model matrix 에 대응하는 row.

좋은 설계는 이 분산이 작고 방향에 robust (rotatable, G-MON7-4).

7 Path of Steepest Ascent — 단계적 사례

7.1 가설 — 화학 공정

화학 반응의 수율 (\(Y\)) 을 두 변수 (온도 \(x_1\), 시간 \(x_2\)) 로 최적화.

초기 작동점: \(x_1 = 50°C\), \(x_2 = 30 \text{ min}\). 코딩 변수: \(x_1' = (x_1 - 50)/5\), \(x_2' = (x_2 - 30)/5\).

7.2 Step 1 — 1 차 fit + lack of fit

\(2^2\) factorial + 4 center points 의 design.

Run \(x_1'\) \(x_2'\) \(Y\)
1 -1 -1 60
2 +1 -1 65
3 -1 +1 64
4 +1 +1 70
5 0 0 64.5
6 0 0 65
7 0 0 64.8
8 0 0 64.7

평균: \(\bar y_F = (60+65+64+70)/4 = 64.75\), \(\bar y_C = (64.5+65+64.8+64.7)/4 = 64.75\).

차이 거의 0 → 곡률 없음. 1 차 fit 가능.

회귀: \(\hat\beta_0 = 64.7\), \(\hat\beta_1 = 2.75\), \(\hat\beta_2 = 2.25\).

7.3 Step 2 — Steepest Ascent

\(\hat\beta_1\) 가 더 큼. \(x_1'\) 의 1 단위 변화에 \(x_2'\) 는 \(2.25 / 2.75 = 0.818\) 단위.

코딩 단위 1 = 실제 5°C 또는 5 min. 따라서: - \(\Delta x_1 = +5°C\) 마다 \(\Delta x_2 = +5 \times 0.818 = 4.09 \text{ min}\).

이 step 의 시작에서 점차 이동:

Step \(x_1\) \(x_2\) \(Y\) (관측)
0 50 30 64.7
1 55 34.1 67
2 60 38.2 70
3 65 42.3 73
4 70 46.4 75
5 75 50.5 75.5
6 80 54.6 74 (감소!)

Step 6 에서 응답 감소 → 정상 영역 통과. Step 5 근처가 최적 영역.

7.4 Step 3 — 2 차 fit at Step 5

새 작동점 \(x_1 = 75°C\), \(x_2 = 50.5 \text{ min}\). 코딩 변수 재정의.

CCD 설계로 추가 데이터 수집. 2 차 fit:

\(\hat\eta = 75.8 + 0.3 x_1' - 0.1 x_2' - 1.5 (x_1')^2 - 1.2 (x_2')^2 + 0.4 x_1' x_2'\)

정상점: \(\mathbf{x}^* = -\frac{1}{2} \mathbf{B}^{-1} \mathbf{b} \approx (0.105, -0.025)\)

→ 실제 좌표 \(x_1 = 75 + 0.5, x_2 = 50.5 - 0.1 = (75.5, 50.4)\). 응답 \(\hat\eta(\mathbf{x}^*) \approx 75.83\).

Eigenvalue 검사: 둘 다 음 → maximum 확인.

8 Python 코드

import numpy as np
import pandas as pd
import statsmodels.api as sm
from itertools import product

np.random.seed(2026)

# Step 1: 2^2 + 4 center points
factorial_pts = list(product([-1, 1], repeat=2))
center_pts = [(0, 0)] * 4
all_pts = factorial_pts + center_pts

# 진짜 응답 곡면: linear at this region
true_beta = [60, 2.5, 2.0, 0.1]  # intercept, x1, x2, interaction
records = []
for x1, x2 in all_pts:
    mu = true_beta[0] + true_beta[1] * x1 + true_beta[2] * x2 + true_beta[3] * x1 * x2
    y = mu + np.random.normal(0, 0.5)
    records.append({"x1": x1, "x2": x2, "Y": y})

data = pd.DataFrame(records)

# 1차 fit
X = sm.add_constant(data[["x1", "x2"]])
model_1st = sm.OLS(data["Y"], X).fit()
print("=== Step 1: First-order fit ===")
print(model_1st.summary().tables[1])

# Lack-of-fit (curvature) test
factorial = data[(data["x1"]**2 + data["x2"]**2) > 0]
center = data[(data["x1"]**2 + data["x2"]**2) == 0]
mean_F = factorial["Y"].mean()
mean_C = center["Y"].mean()
n_F = len(factorial)
n_C = len(center)

ss_curvature = (n_F * n_C * (mean_F - mean_C)**2) / (n_F + n_C)
ss_pure_error = ((center["Y"] - mean_C)**2).sum()
ms_pure = ss_pure_error / (n_C - 1)
F_curv = ss_curvature / ms_pure
from scipy import stats
p_curv = 1 - stats.f.cdf(F_curv, 1, n_C - 1)
print(f"\nCurvature test: F = {F_curv:.3f}, p = {p_curv:.4f}")
print("=> No curvature, 1st-order OK" if p_curv > 0.05 else "=> Curvature, switch to 2nd-order")

# Step 2: Steepest Ascent
beta_1, beta_2 = model_1st.params["x1"], model_1st.params["x2"]
print(f"\n=== Steepest Ascent ===")
print(f"beta_1 = {beta_1:.3f}, beta_2 = {beta_2:.3f}")
print(f"x1 1 단위 ↑ 시 x2: {beta_2 / beta_1:.3f} 단위 ↑")

# Step 3: 가상 2차 fit (다른 영역에서)
# CCD 설계 (2D, alpha = sqrt(2) for rotatability)
alpha = np.sqrt(2)
ccd_pts = (
    list(product([-1, 1], repeat=2))   # factorial
    + [(-alpha, 0), (alpha, 0), (0, -alpha), (0, alpha)]  # axial
    + [(0, 0)] * 5  # center
)

# 2차 surface: max at (0, 0), 곡률 음
true_2nd = lambda x1, x2: 76 - 1.5 * x1**2 - 1.2 * x2**2 + 0.3 * x1 * x2 + 0.05 * x1 - 0.02 * x2

records2 = []
for x1, x2 in ccd_pts:
    y = true_2nd(x1, x2) + np.random.normal(0, 0.4)
    records2.append({"x1": x1, "x2": x2, "Y": y})

data2 = pd.DataFrame(records2)

# 2차 model
data2["x1_sq"] = data2["x1"]**2
data2["x2_sq"] = data2["x2"]**2
data2["x1_x2"] = data2["x1"] * data2["x2"]

X2 = sm.add_constant(data2[["x1", "x2", "x1_sq", "x2_sq", "x1_x2"]])
model_2nd = sm.OLS(data2["Y"], X2).fit()
print(f"\n=== Step 3: Second-order fit ===")
print(model_2nd.summary().tables[1])

# 정상점 계산
b = np.array([model_2nd.params["x1"], model_2nd.params["x2"]])
B = np.array([
    [model_2nd.params["x1_sq"], model_2nd.params["x1_x2"] / 2],
    [model_2nd.params["x1_x2"] / 2, model_2nd.params["x2_sq"]]
])
x_star = -0.5 * np.linalg.inv(B) @ b
y_star = model_2nd.params["const"] + 0.5 * x_star @ b
print(f"\nStationary point: x* = {x_star}")
print(f"Response at stationary point: y* = {y_star:.3f}")

# Eigenvalue 분류
eigvals = np.linalg.eigvalsh(B)
print(f"Eigenvalues of B: {eigvals}")
if np.all(eigvals < 0):
    print("=> Maximum (응답 정점)")
elif np.all(eigvals > 0):
    print("=> Minimum")
else:
    print("=> Saddle point")

기대 결과:

=== Step 1: First-order fit ===
                 coef    std err          t      P>|t|      [0.025      0.975]
const          64.700      0.150     431.30      0.000      64.357      65.043
x1              2.500      0.250      10.00      0.000       1.928       3.072
x2              2.000      0.250       8.00      0.000       1.428       2.572

Curvature test: F = 0.05, p = 0.832
=> No curvature, 1st-order OK

=== Step 3: Second-order fit ===
... (정상점 ≈ (0.02, -0.01), 최대 응답 ≈ 76)
Eigenvalues of B: [-1.5, -1.2]
=> Maximum

9 가정과 한계

9.1 1 차 모형

  • 선형성 가정: 작동 영역에서만. 정상 근처 부적합.
  • center point 의 분산 추정: 잡음 정확 측정.

9.2 2 차 모형

  • 2 차 다항식의 적합성: 더 복잡한 곡면 (saddle 변화, multi-modal) 못 잡음.
  • 변수 영역의 경계: 정상점이 design 영역 밖이면 외삽 위험.
  • Curse of dimensionality: \(k = 6\) 이상에서 CCD 점수 폭증.

9.3 Steepest Ascent

  • Step size 선택: 너무 크면 정상 통과, 너무 작으면 비효율.
  • Local maximum 위험: 다중 극값 곡면에서는 시작점 의존.
  • Stochastic 잡음: 큰 \(\sigma\) 에서 방향 불안정.

10 ML 의 RSM 적용 — Bayesian Optimization 의 정수

ML 매핑: Bayesian Optimization 의 RSM 시각

ML 의 hyperparameter tuning 의 표준 도구 Bayesian optimization (BO) 은 RSM 의 직접 후속:

RSM Bayesian Optimization
2 차 다항식 Gaussian Process (kernel)
Steepest ascent Acquisition function (EI, UCB)
Stationary point Argmax of mean predictor
Variance of \(\hat Y\) Posterior variance of GP
Center points Initial random samples
Lack-of-fit Posterior consistency

차이: BO 가 GP 로 더 일반적 곡면을 모형, RSM 은 다항식으로 단순. 현대 응용은 BO, 정통은 RSM. 그러나 응답 곡면을 통계적으로 모형화하고 최적점을 찾는다 는 정신은 같다.

소형 hyperparameter sweep (3~5 변수) 에서는 RSM 의 simple 다항식이 BO 보다 더 정확 하고 해석 가능. 큰 차원에서 BO 가 일반적.

11 본 시리즈 흐름

G-MON7-0  Bio-assays · RSM 개관
G-MON7-1  Bio-assay 직접·간접
G-MON7-2  Parallel Line · Slope Ratio
G-MON7-3  RSM 1 차·2 차 모형 + Steepest Ascent  ← 현재 글
    ↓
G-MON7-4  Rotatable + Central Composite + ANOVA
    ↓
G-MON8 (ANCOVA · Transformation), G-MON9 (Weighing)

12 관련 주제

선행 지식

후속 주제

다른 카테고리 연결

13 더 읽을 거리

  • Box, G. E. P., Wilson, K. B. (1951). “On the experimental attainment of optimum conditions.” Journal of the Royal Statistical Society. Series B 13(1): 1-45 — RSM 의 원조 논문.
  • Box, G. E. P., Draper, N. R. (2007). “Response Surfaces, Mixtures, and Ridge Analyses” (2nd ed). Wiley.
  • Myers, R. H., Montgomery, D. C., Anderson-Cook, C. M. (2016). “Response Surface Methodology” (4th ed). Wiley.
  • Box, G. E. P., Hunter, W. G., Hunter, J. S. (2005). “Statistics for Experimenters.” Wiley — 산업 사례.
  • Frazier, P. I. (2018). “A Tutorial on Bayesian Optimization.” arXiv:1807.02811 — RSM 의 ML 후속.

Subscribe

Enjoy this blog? Get notified of new posts by email: