Python 완전정복 시리즈 #27: SciPy 완벽 이해하기

지난 시간에는 NumPy를 이용해 수백만 개의 숫자를 빠르게 계산하는 방법을 배웠습니다.

배열을 만들고, 형태를 변경하고, 평균과 표준편차를 계산했으며, 행렬 곱셈과 브로드캐스팅까지 살펴보았습니다.

그런데 실무에서는 단순한 계산만으로 해결되지 않는 문제가 등장합니다.

비용이 가장 낮아지는 생산량은 얼마일까?

측정하지 않은 중간 지점의 값은 얼마일까?

복잡한 함수 아래의 넓이는 얼마일까?

두 그룹의 평균 차이는 우연일까?

센서 데이터에 섞인 잡음을 어떻게 제거할까?

신호에서 가장 중요한 봉우리는 어디일까?

거대한 행렬에 0이 대부분이라면 어떻게 저장할까?

NumPy는 숫자 군단을 빠르게 움직이는 훌륭한 지휘관입니다.

하지만 위와 같은 문제를 해결하려면 분야별 전문 기술자들이 필요합니다.

최적화 전문가
적분 전문가
통계 분석가
신호 처리 엔지니어
공간 계산 전문가
희소행렬 관리자

이 전문가들이 한 연구소에 모여 있는 라이브러리가 바로 SciPy입니다. 🧪🔬

NumPy가 튼튼한 과학 계산 엔진이라면 SciPy는 그 엔진에 다양한 실험 장비를 연결한 종합 연구소입니다.

1. SciPy란 무엇인가?

SciPy는 Python에서 과학·공학 계산에 필요한 알고리즘과 고수준 함수를 제공하는 오픈소스 라이브러리입니다.

NumPy 배열을 기반으로 동작하며 최적화, 적분, 보간, 선형대수, 통계, 신호 처리, 미분방정식, 푸리에 변환, 공간 계산, 희소배열 등 다양한 기능을 제공합니다. 2026년 7월 기준 공식 문서의 최신 안정 버전은 SciPy 1.18.0입니다.

SciPy가 활용되는 분야는 매우 넓습니다.

물리학
화학
생명과학
천문학
금융공학
기계공학
전자공학
통계 분석
머신러닝
이미지 처리
음향 분석
센서 데이터 분석

SciPy는 단순히 공식을 모아 놓은 라이브러리가 아닙니다.

수십 년 동안 연구되고 검증된 수치 알고리즘을 Python에서 비교적 간단한 코드로 사용할 수 있도록 연결해 줍니다. 일부 기능은 C, C++, Fortran으로 구현된 고성능 코드도 활용합니다.

2. NumPy와 SciPy의 차이

NumPy와 SciPy는 경쟁 관계가 아닙니다.

두 라이브러리는 역할이 다릅니다.

NumPy는 다음 작업에 강합니다.

다차원 배열 생성
배열 인덱싱과 슬라이싱
요소별 수학 연산
통계 계산
행렬 곱셈
난수 생성
배열 저장과 불러오기

SciPy는 NumPy 배열을 사용하면서 더욱 전문적인 문제를 해결합니다.

함수의 최솟값 찾기
비선형 방정식 풀기
정적분 계산
미분방정식 풀이
측정값 사이의 중간값 추정
통계적 가설 검정
신호 필터링
피크 탐지
주파수 분석
희소행렬 계산
공간 거리 계산

간단히 정리하면 다음과 같습니다.

NumPy
→ 숫자를 저장하고 빠르게 계산한다.

SciPy
→ 그 숫자를 이용해 전문 문제를 해결한다.

NumPy가 벽돌과 철근을 공급한다면 SciPy는 연구소, 다리, 발전소를 설계하는 전문 기술팀입니다.

3. SciPy 설치하기

터미널이나 명령 프롬프트에서 다음 명령어를 실행합니다.

python -m pip install scipy

SciPy 공식 초보자 설치 안내에서도 가상환경 안에서 `python -m pip install scipy`를 사용하는 방식을 안내합니다.

python -m venv venv
venv\Scripts\activate
python -m pip install scipy

python3 -m venv venv
source venv/bin/activate
python -m pip install scipy

이번 글의 실습에서는 NumPy와 Matplotlib도 사용합니다.

python -m pip install numpy scipy matplotlib

import scipy

print(scipy.__version__)

버전이 정상적으로 출력되면 설치가 완료된 것입니다.

4. SciPy의 주요 하위 패키지

SciPy는 기능별 하위 패키지로 구성되어 있습니다.

SciPy 공식 사용자 가이드도 기능 영역별 하위 패키지 구조를 중심으로 구성되어 있습니다.

모든 기능을 한꺼번에 외울 필요는 없습니다.

문제를 만났을 때 필요한 연구실 문을 찾아가면 됩니다.

최솟값을 찾는다
→ optimize 연구실

넓이를 계산한다
→ integrate 연구실

잡음을 제거한다
→ signal 연구실

평균 차이를 검정한다
→ stats 연구실

5. SciPy를 불러오는 권장 방식

SciPy 전체를 불러온 뒤 모든 함수를 찾기보다는 필요한 하위 패키지나 함수를 직접 불러오는 방식이 읽기 좋습니다.

from scipy import optimize
from scipy import integrate
from scipy import stats

사용:

result = optimize.minimize_scalar(
    lambda x: (x - 3) ** 2
)

from scipy.optimize import minimize_scalar
from scipy.integrate import quad

사용:

result = minimize_scalar(
    lambda x: (x - 3) ** 2
)

작은 예제에서는 필요한 함수를 직접 불러오는 것이 편리합니다.

from scipy.optimize import minimize

여러 함수를 함께 사용할 때는 하위 패키지 단위로 불러오면 출처가 명확합니다.

from scipy import optimize

result = optimize.minimize(...)

6. 물리·수학 상수 사용하기

SciPy는 자주 사용하는 물리 상수와 단위를 제공합니다.

from scipy import constants

print(constants.c)
print(constants.g)
print(constants.pi)

대표적인 값:

constants.c
→ 진공에서의 빛의 속도

constants.g
→ 표준 중력가속도

constants.pi
→ 원주율

from scipy.constants import value

planck_constant = value(
    "Planck constant"
)

print(planck_constant)

`scipy.constants.value()`는 물리 상수 데이터베이스에서 이름에 해당하는 값을 반환합니다.

물체가 정지 상태에서 낙하한다고 가정하면 다음 식으로 거리를 근사할 수 있습니다.

거리 = 1 / 2 × 중력가속도 × 시간²

from scipy import constants

time_seconds = 3

distance = (
    0.5
    * constants.g
    * time_seconds ** 2
)

print(
    f"{time_seconds}초 동안 "
    f"약 {distance:.2f}m 낙하합니다."
)

결과 예시:

3초 동안 약 44.13m 낙하합니다.

직접 숫자 `9.80665`를 입력하는 대신 이름이 있는 상수를 사용하면 코드의 의미가 명확해집니다.

7. 최적화란 무엇인가?

최적화는 주어진 조건 안에서 함수의 값이 가장 작거나 가장 커지는 지점을 찾는 과정입니다.

실무에서는 다음과 같은 문제로 등장합니다.

비용을 최소화한다.

수익을 최대화한다.

배송 거리를 최소화한다.

모델의 예측 오차를 최소화한다.

재료 사용량을 최소화한다.

에너지 소비를 최소화한다.

SciPy의 `optimize` 하위 패키지는 함수 최소화, 전역 최적화, 선형계획법, 비선형 최소제곱, 방정식의 해, 곡선 맞춤 등을 제공합니다.

최적화 문제에는 보통 세 가지 요소가 있습니다.

최소화하거나 최대화할 대상입니다.

생산 비용
예측 오차
이동 거리
연료 소비량

우리가 조절할 수 있는 값입니다.

생산량
가격
속도
온도
광고비

변수가 지켜야 할 범위나 규칙입니다.

생산량은 0 이상

광고비는 1,000만 원 이하

재료 사용량은 재고보다 작아야 함

최적화 알고리즘은 산속에서 가장 낮은 골짜기를 찾는 탐험가와 비슷합니다.

문제는 안개가 심해 산 전체를 한눈에 볼 수 없다는 점입니다. ⛰️

8. 한 변수 함수의 최솟값 찾기

다음 함수를 생각해 보겠습니다.

f(x) = (x - 3)² + 2

이 함수는 `x=3`일 때 가장 작은 값을 가집니다.

SciPy로 확인해 보겠습니다.

from scipy.optimize import minimize_scalar


def objective(x: float) -> float:
    return (x - 3) ** 2 + 2


result = minimize_scalar(
    objective
)

print(result)

필요한 값만 출력합니다.

print(
    f"최적의 x: {result.x:.6f}"
)

print(
    f"함수 최솟값: {result.fun:.6f}"
)

print(
    f"성공 여부: {result.success}"
)

결과 예시:

최적의 x: 3.000000
함수 최솟값: 2.000000
성공 여부: True

`minimize_scalar()`는 한 변수 스칼라 함수의 국소 최솟값을 찾는 함수이며 Brent, bounded, golden 등의 방법을 지원합니다.

result.x
→ 최적의 변수 값

result.fun
→ 목적함수의 최솟값

result.success
→ 계산 성공 여부

result.message
→ 종료 이유

result.nfev
→ 함수 평가 횟수

최적화 결과는 숫자만 확인하지 말고 `success`와 `message`도 확인하는 습관이 좋습니다.

탐험가가 골짜기를 찾았다고 무전했는데 실제로는 배터리가 꺼진 것일 수도 있기 때문입니다.

9. 여러 변수 함수 최적화하기

변수가 여러 개라면 `minimize()`를 사용할 수 있습니다.

다음 함수의 최솟값을 찾아 보겠습니다.

f(x, y) = (x - 2)² + (y + 1)²

최솟값은 다음 위치에 있습니다.

x = 2
y = -1

import numpy as np
from scipy.optimize import minimize


def objective(
    variables: np.ndarray
) -> float:
    x, y = variables

    return (
        (x - 2) ** 2
        + (y + 1) ** 2
    )


initial_guess = np.array([
    0.0,
    0.0
])

result = minimize(
    objective,
    initial_guess
)

print(
    f"최적의 변수: {result.x}"
)

print(
    f"최솟값: {result.fun:.6f}"
)

print(
    f"성공 여부: {result.success}"
)

결과 예시:

최적의 변수: [ 1.99999999 -1.00000001]
최솟값: 0.000000
성공 여부: True

`minimize()`는 하나 이상의 변수를 받는 스칼라 목적함수의 최소화를 수행하며 BFGS, Nelder-Mead, L-BFGS-B, SLSQP 등 다양한 알고리즘을 지원합니다.

initial_guess = np.array([
    0.0,
    0.0
])

최적화 알고리즘은 일반적으로 초기 위치에서 탐색을 시작합니다.

복잡한 함수에서는 초기값에 따라 서로 다른 국소 최솟값에 도착할 수도 있습니다.

첫 번째 등산로
→ 작은 골짜기 발견

두 번째 등산로
→ 더 깊은 골짜기 발견

따라서 중요한 문제에서는 여러 초기값으로 계산해 결과를 비교하는 방법도 검토해야 합니다.

10. 범위가 있는 최적화

실무 변수에는 허용 범위가 있는 경우가 많습니다.

생산량은 0개 이상

온도는 10℃ 이상 100℃ 이하

광고비는 0원 이상 1억 원 이하

다음 함수에서 `x`를 0부터 5 사이로 제한해 보겠습니다.

from scipy.optimize import minimize_scalar


def objective(x: float) -> float:
    return (x - 8) ** 2


result = minimize_scalar(
    objective,
    bounds=(0, 5),
    method="bounded"
)

print(
    f"범위 안의 최적값: {result.x:.6f}"
)

제약이 없다면 최솟값은 `x=8`이지만 허용 범위가 0부터 5이므로 경계에 가까운 `x=5`가 선택됩니다.

import numpy as np
from scipy.optimize import minimize


def objective(
    variables: np.ndarray
) -> float:
    x, y = variables

    return (
        (x - 8) ** 2
        + (y - 3) ** 2
    )


bounds = [
    (0, 5),
    (0, 10)
]

result = minimize(
    objective,
    x0=np.array([1.0, 1.0]),
    bounds=bounds,
    method="L-BFGS-B"
)

print(result.x)

결과는 대략 다음과 같습니다.

[5. 3.]

x는 원래 8이 최적
→ 최대 허용값 5 선택

y는 3이 최적
→ 범위 안이므로 3 선택

L-BFGS-B는 변수의 상한과 하한이 있는 최적화에 사용할 수 있는 방법입니다.

11. 최대값을 찾는 방법

SciPy의 주요 최적화 함수는 기본적으로 최솟값을 찾습니다.

최대값을 찾고 싶다면 함수에 음수를 붙여 최소화할 수 있습니다.

다음 수익 함수를 생각해 보겠습니다.

수익(x) = -2x² + 40x - 50

from scipy.optimize import minimize_scalar


def profit(x: float) -> float:
    return (
        -2 * x ** 2
        + 40 * x
        - 50
    )


def negative_profit(x: float) -> float:
    return -profit(x)


result = minimize_scalar(
    negative_profit,
    bounds=(0, 20),
    method="bounded"
)

optimal_x = result.x
maximum_profit = profit(optimal_x)

print(
    f"최적 생산량: {optimal_x:.2f}"
)

print(
    f"최대 수익: {maximum_profit:.2f}"
)

원리:

가장 큰 수익
=
음수로 바꿨을 때 가장 작은 값

수익 100 → -100
수익 200 → -200

-200이 더 작으므로
원래 수익 200이 선택됨

최대값을 찾기 위해 그래프를 뒤집어 놓는 셈입니다.

산봉우리를 찾기 어렵다면 지도를 뒤집어 골짜기로 만들어 버립니다. 🗺️

12. 방정식의 해 찾기

다음 방정식을 풀어 보겠습니다.

x² - 4 = 0

해는 다음과 같습니다.

x = -2
x = 2

`root_scalar()`를 이용해 양수 해를 찾겠습니다.

from scipy.optimize import root_scalar


def equation(x: float) -> float:
    return x ** 2 - 4


result = root_scalar(
    equation,
    bracket=(0, 5),
    method="brentq"
)

print(
    f"해: {result.root:.6f}"
)

print(
    f"수렴 여부: {result.converged}"
)

print(
    f"반복 횟수: {result.iterations}"
)

`root_scalar()`는 한 변수 함수가 0이 되는 지점을 찾으며 이분법, Brent 방법, Newton 방법, 할선법 등 여러 알고리즘을 지원합니다.

bracket=(0, 5)

해가 존재한다고 예상되는 구간입니다.

f(0) = -4
f(5) = 21

함수값의 부호가 바뀌므로 그 사이에 해가 있다고 판단할 수 있습니다.

음수에서 양수로 넘어가는 길목을 검문해 정확히 0이 되는 지점을 찾는 과정입니다.

result = root_scalar(
    equation,
    bracket=(-5, 0),
    method="brentq"
)

print(result.root)

결과:

-2.0

13. 곡선 맞춤과 매개변수 추정

실험 데이터는 완벽한 직선이나 곡선 위에 놓이지 않습니다.

측정 오차와 잡음이 섞여 있기 때문입니다.

다음 모델을 데이터에 맞춰 보겠습니다.

y = ax + b

우리가 찾을 값은 기울기 `a`와 절편 `b`입니다.

import numpy as np

x_data = np.array([
    0,
    1,
    2,
    3,
    4,
    5
], dtype=float)

y_data = np.array([
    2.1,
    5.0,
    8.2,
    11.1,
    13.9,
    17.2
])

대략 다음 관계처럼 보입니다.

y ≈ 3x + 2

import numpy as np
from scipy.optimize import curve_fit


def linear_model(
    x: np.ndarray,
    slope: float,
    intercept: float
) -> np.ndarray:
    return (
        slope * x
        + intercept
    )


parameters, covariance = curve_fit(
    linear_model,
    x_data,
    y_data
)

slope, intercept = parameters

print(
    f"추정 기울기: {slope:.4f}"
)

print(
    f"추정 절편: {intercept:.4f}"
)

`curve_fit()`은 비선형 최소제곱법을 사용해 주어진 모델 함수가 관측 데이터에 맞도록 매개변수를 추정합니다. 모델 함수는 첫 번째 인수로 독립변수를 받고 나머지 인수로 추정할 매개변수를 받습니다.

import matplotlib.pyplot as plt
import numpy as np

x_fine = np.linspace(
    x_data.min(),
    x_data.max(),
    200
)

y_fitted = linear_model(
    x_fine,
    slope,
    intercept
)

fig, ax = plt.subplots()

ax.scatter(
    x_data,
    y_data,
    label="측정 데이터"
)

ax.plot(
    x_fine,
    y_fitted,
    label="추정 모델"
)

ax.set_title(
    "측정 데이터 곡선 맞춤"
)

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.legend()
ax.grid(alpha=0.3)

plt.show()

`curve_fit()`은 측정값을 억지로 한 줄에 세우는 것이 아닙니다.

전체 오차가 작아지도록 가장 설득력 있는 곡선을 찾는 중재자입니다.

14. 수치 적분 이해하기

적분은 함수 아래의 넓이를 계산하는 방법으로 자주 설명됩니다.

다음 함수를 생각해 보겠습니다.

f(x) = x²

0부터 2까지의 정적분은 다음과 같습니다.

∫₀² x² dx = 8 / 3

정확한 결과:

2.666666...

함수가 복잡하거나 측정 데이터로 주어지면 손으로 적분하기 어렵습니다.

이럴 때 수치 적분을 사용합니다.

수치 적분은 구간을 작은 조각으로 나누고 각 조각의 넓이를 합산해 전체 값을 근사합니다.

커다란 모양의 넓이
→ 작은 조각으로 분할
→ 각 조각의 넓이 계산
→ 모두 합산

SciPy의 `integrate` 하위 패키지는 정적분과 미분방정식 계산을 위한 여러 기법을 제공합니다.

15. 정적분 계산하기

`quad()`를 사용해 `x²`을 0부터 2까지 적분해 보겠습니다.

from scipy.integrate import quad


def function(x: float) -> float:
    return x ** 2


integral, error = quad(
    function,
    0,
    2
)

print(
    f"적분 결과: {integral:.10f}"
)

print(
    f"추정 오차: {error:.10e}"
)

결과 예시:

적분 결과: 2.6666666667
추정 오차: 2.9605947323e-14

`quad()`는 한 변수 함수의 정적분을 계산하고 적분 결과와 절대오차 추정값을 반환합니다. 내부적으로 QUADPACK 계열 알고리즘을 활용합니다.

integral, error = quad(...)

integral
→ 계산된 적분값

error
→ 절대오차 추정값

오차가 0이라고 단정하지 않고 추정 오차까지 함께 제공한다는 점이 중요합니다.

수치 계산은 “정답입니다”보다 “이 정도 오차 범위로 계산되었습니다”라는 태도를 가집니다.

16. 무한 구간 적분하기

SciPy는 무한 구간 적분도 지원합니다.

다음 가우스 함수의 전체 적분을 계산해 보겠습니다.

e^(-x²)

마이너스 무한대부터 플러스 무한대까지 적분한 값은 `√π`입니다.

import numpy as np
from scipy.integrate import quad


def gaussian(x: float) -> float:
    return np.exp(
        -(x ** 2)
    )


integral, error = quad(
    gaussian,
    -np.inf,
    np.inf
)

print(
    f"적분 결과: {integral:.10f}"
)

print(
    f"sqrt(pi): {np.sqrt(np.pi):.10f}"
)

`quad()`에서는 `-np.inf`와 `np.inf`를 적분 구간으로 지정할 수 있습니다.

무한대를 직접 잘게 나누는 것이 아니라 무한 구간을 처리할 수 있는 수치 기법을 사용합니다.

물론 모든 함수가 무한 구간에서 적분 가능한 것은 아닙니다.

함수가 발산하거나 특이점이 심하면 계산 결과와 경고를 꼼꼼하게 확인해야 합니다.

17. 미분방정식 풀기

미분방정식은 어떤 값 자체보다 그 값이 시간에 따라 어떻게 변하는지를 표현합니다.

예를 들어 물체의 온도가 주변 온도에 가까워지는 과정은 다음과 같이 모델링할 수 있습니다.

dT / dt = -k(T - T환경)

T
→ 물체의 현재 온도

T환경
→ 주변 온도

k
→ 냉각 속도 계수

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp


ROOM_TEMPERATURE = 25.0
COOLING_RATE = 0.15


def cooling_model(
    time: float,
    temperature: np.ndarray
) -> np.ndarray:
    return (
        -COOLING_RATE
        * (
            temperature
            - ROOM_TEMPERATURE
        )
    )


time_points = np.linspace(
    0,
    30,
    200
)

solution = solve_ivp(
    cooling_model,
    t_span=(0, 30),
    y0=[90.0],
    t_eval=time_points
)

print(
    f"계산 성공: {solution.success}"
)

print(
    f"30분 후 온도: "
    f"{solution.y[0, -1]:.2f}℃"
)

그래프로 표현합니다.

fig, ax = plt.subplots()

ax.plot(
    solution.t,
    solution.y[0]
)

ax.axhline(
    ROOM_TEMPERATURE,
    linestyle="--",
    label="주변 온도"
)

ax.set_title(
    "시간에 따른 물체의 냉각"
)

ax.set_xlabel("시간")
ax.set_ylabel("온도")
ax.legend()
ax.grid(alpha=0.3)

plt.show()

`solve_ivp()`는 초기값이 주어진 상미분방정식 시스템을 수치적으로 풉니다. 기본 방법은 RK45이며 비강성 문제에는 명시적 Runge-Kutta 계열, 강성 문제에는 Radau나 BDF 같은 방법을 검토할 수 있습니다.

solve_ivp(
    fun,
    t_span,
    y0,
    t_eval
)

fun
→ 미분방정식 함수

t_span
→ 계산할 시작 시간과 종료 시간

y0
→ 초기 상태

t_eval
→ 결과를 저장할 시간 지점

미분방정식은 미래를 정확히 예언하는 수정구슬이 아닙니다.

우리가 만든 변화 규칙을 따라 시스템이 어떻게 움직일지 계산하는 시뮬레이터입니다. 🔮

18. 보간법이란?

다음 측정 데이터가 있다고 가정해 보겠습니다.

그런데 4시의 온도는 측정하지 못했습니다.

0시  3시  ?  6시  9시
10℃  13℃  ?  19℃  22℃

보간은 알고 있는 데이터 사이의 값을 추정하는 방법입니다.

보간
→ 측정 범위 안의 빈 값 추정

외삽
→ 측정 범위 바깥의 값 추정

외삽은 관측 범위를 넘어가기 때문에 일반적으로 보간보다 위험합니다.

측정한 길 안에서 빈 발자국을 연결하는 것과, 지도 밖으로 길을 상상하는 것은 위험도가 다릅니다.

19. PCHIP 보간으로 중간값 추정하기

PCHIP는 데이터의 형태와 단조성을 비교적 잘 보존하는 조각별 3차 보간법입니다.

SciPy의 `PchipInterpolator`는 단조 3차 스플라인을 이용하는 형태 보존 보간기를 제공합니다. 입력 x는 증가하는 값이어야 하며 중복값을 포함할 수 없습니다.

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import PchipInterpolator


hours = np.array([
    0,
    3,
    6,
    9,
    12
], dtype=float)

temperatures = np.array([
    10,
    13,
    19,
    22,
    24
], dtype=float)

interpolator = PchipInterpolator(
    hours,
    temperatures
)

estimated_temperature = interpolator(
    4.0
)

print(
    f"4시 추정 온도: "
    f"{estimated_temperature:.2f}℃"
)

fine_hours = np.linspace(
    hours.min(),
    hours.max(),
    300
)

fine_temperatures = interpolator(
    fine_hours
)

fig, ax = plt.subplots()

ax.scatter(
    hours,
    temperatures,
    label="측정값"
)

ax.plot(
    fine_hours,
    fine_temperatures,
    label="PCHIP 보간"
)

ax.scatter(
    [4],
    [estimated_temperature],
    marker="x",
    s=100,
    label="4시 추정값"
)

ax.set_title(
    "시간별 온도 보간"
)

ax.set_xlabel("시간")
ax.set_ylabel("온도")
ax.legend()
ax.grid(alpha=0.3)

plt.show()

PCHIP에 전달하는 x 값은 증가 순서로 정렬되어야 합니다.

hours = np.array([
    0,
    3,
    6,
    9
])

다음과 같은 중복 x는 사용할 수 없습니다.

hours = np.array([
    0,
    3,
    3,
    6
])

같은 시간에 서로 다른 두 온도를 전달하면 함수가 어느 값을 선택해야 할지 결정할 수 없기 때문입니다.

20. 통계 요약하기

SciPy의 `stats` 하위 패키지는 확률분포, 요약 통계, 상관관계, 가설 검정 등 다양한 통계 기능을 제공합니다.

import numpy as np

scores = np.array([
    72,
    85,
    91,
    68,
    77,
    88,
    95,
    81
], dtype=float)

from scipy.stats import describe

summary = describe(
    scores
)

print(summary)

주요 속성을 확인합니다.

print(
    f"데이터 개수: {summary.nobs}"
)

print(
    f"최솟값: {summary.minmax[0]}"
)

print(
    f"최댓값: {summary.minmax[1]}"
)

print(
    f"평균: {summary.mean:.2f}"
)

print(
    f"분산: {summary.variance:.2f}"
)

print(
    f"왜도: {summary.skewness:.2f}"
)

print(
    f"첨도: {summary.kurtosis:.2f}"
)

`scipy.stats.describe()`는 관측치 수, 최소·최대, 평균, 분산, 왜도, 첨도 등의 기술통계를 계산합니다.

분포가 한쪽으로 얼마나 치우쳐 있는지를 나타냅니다.

왜도 ≈ 0
→ 좌우가 비교적 대칭

왜도 > 0
→ 오른쪽 꼬리가 긴 경향

왜도 < 0
→ 왼쪽 꼬리가 긴 경향

분포의 꼬리와 중심 집중 형태를 설명하는 통계량입니다.

다만 왜도와 첨도 하나만으로 데이터의 성격을 단정해서는 안 됩니다.

히스토그램, 상자 수염 그래프, 표본 크기 등을 함께 살펴보는 것이 좋습니다.

21. 독립표본 t 검정

두 개의 서로 독립적인 그룹 평균에 차이가 있는지 분석해 보겠습니다.

예를 들어 새로운 학습법을 사용한 A그룹과 기존 학습법을 사용한 B그룹의 점수가 있습니다.

import numpy as np

group_a = np.array([
    82,
    85,
    88,
    90,
    86,
    84,
    91,
    87
], dtype=float)

group_b = np.array([
    76,
    78,
    80,
    77,
    81,
    79,
    75,
    82
], dtype=float)

from scipy.stats import ttest_ind

result = ttest_ind(
    group_a,
    group_b,
    equal_var=False
)

print(
    f"t 통계량: "
    f"{result.statistic:.4f}"
)

print(
    f"p-value: "
    f"{result.pvalue:.6f}"
)

`ttest_ind()`는 독립된 두 표본의 평균에 관한 t 검정을 수행합니다. `equal_var=False`를 지정하면 두 모집단의 분산이 같다고 가정하지 않는 Welch t 검정을 수행합니다.

alpha = 0.05

if result.pvalue < alpha:
    print(
        "설정한 유의수준에서 "
        "평균 차이가 통계적으로 "
        "유의하다고 판단합니다."
    )
else:
    print(
        "설정한 유의수준에서 "
        "평균 차이가 유의하다고 "
        "판단할 근거가 부족합니다."
    )

22. p-value를 해석할 때 주의할 점

p-value는 다음 의미가 아닙니다.

❌ 가설이 틀릴 확률

❌ 결과가 우연일 확률

❌ 효과가 클 확률

❌ 연구가 성공할 확률

p-value는 귀무가설과 검정의 가정이 성립한다고 할 때, 관측한 결과만큼 또는 그보다 더 극단적인 결과가 나타날 가능성을 계산한 값입니다.

귀무가설 아래에서는
현재와 같은 결과가 비교적 드물다.

귀무가설을 기각할 만큼
충분한 증거를 얻지 못했다.

큰 p-value가 두 그룹이 완전히 같다는 사실을 증명하는 것은 아닙니다.

표본 크기
효과 크기
신뢰구간
데이터 분포
이상치
표본 추출 방식
검정의 가정
다중 검정 여부

표본이 매우 크면 실무적으로 작은 차이도 통계적으로 유의하게 나타날 수 있습니다.

따라서 다음 두 질문을 함께 해야 합니다.

통계적으로 차이가 있는가?

실무적으로 의미 있는 크기의 차이인가?

p-value는 판결문이 아니라 증거물 하나입니다. ⚖️

23. 신호 처리란?

신호는 시간이나 공간에 따라 변하는 값을 의미합니다.

마이크로 녹음한 음성
심전도
기계 진동
온도 센서
네트워크 트래픽
주가 변화
빛의 세기

현실의 신호에는 잡음이 포함되는 경우가 많습니다.

원래 신호 + 잡음 = 측정 신호

신호 처리에서는 다음 작업을 수행할 수 있습니다.

잡음 제거
특정 주파수 제거
피크 탐지
패턴 검출
주파수 분석
신호 평활화

SciPy의 `signal` 하위 패키지는 디지털·아날로그 필터 설계와 신호 처리 함수를 제공합니다.

24. 저주파 통과 필터로 잡음 제거하기

저주파 통과 필터는 낮은 주파수 성분은 통과시키고 높은 주파수 성분을 줄입니다.

천천히 변하는 온도 신호에 빠르게 흔들리는 잡음이 섞여 있다면 저주파 통과 필터를 사용할 수 있습니다.

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import (
    butter,
    sosfiltfilt
)


sampling_rate = 100
duration = 5

time = np.linspace(
    0,
    duration,
    sampling_rate * duration,
    endpoint=False
)

clean_signal = np.sin(
    2 * np.pi * 1 * time
)

high_frequency_noise = (
    0.35
    * np.sin(
        2 * np.pi * 20 * time
    )
)

measured_signal = (
    clean_signal
    + high_frequency_noise
)

sos = butter(
    N=4,
    Wn=5,
    btype="lowpass",
    fs=sampling_rate,
    output="sos"
)

filtered_signal = sosfiltfilt(
    sos,
    measured_signal
)

`butter()`는 Butterworth 필터 계수를 설계하며 `fs`를 지정하면 임계주파수 `Wn`을 같은 주파수 단위로 전달할 수 있습니다. `sosfiltfilt()`는 2차 구간 계수로 앞뒤 방향 필터링을 수행합니다.

fig, ax = plt.subplots(
    figsize=(12, 6)
)

ax.plot(
    time,
    measured_signal,
    alpha=0.5,
    label="잡음 포함 신호"
)

ax.plot(
    time,
    filtered_signal,
    linewidth=2,
    label="필터 적용 신호"
)

ax.set_title(
    "저주파 통과 필터를 이용한 잡음 제거"
)

ax.set_xlabel("시간")
ax.set_ylabel("진폭")
ax.legend()
ax.grid(alpha=0.3)

plt.show()

N=4
→ 필터 차수

Wn=5
→ 차단주파수 5Hz

btype="lowpass"
→ 저주파 통과 필터

fs=100
→ 초당 100개 샘플

output="sos"
→ 2차 구간 형식으로 계수 반환

필터가 잡음을 지우는 마법 지우개는 아닙니다.

무엇을 신호로 보고 무엇을 잡음으로 볼 것인지 먼저 결정해야 합니다.

차단주파수를 잘못 선택하면 잡음뿐 아니라 중요한 신호까지 퇴장할 수 있습니다.

25. 신호의 피크 찾기

피크는 주변 값보다 높은 지점을 의미합니다.

다음과 같은 작업에 활용할 수 있습니다.

심박 신호의 박동 위치 탐지

기계 진동의 충격 지점 탐지

하루 중 트래픽이 높은 시간 탐지

매출이 급증한 시점 탐지

스펙트럼의 주요 주파수 탐지

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import find_peaks


signal = np.array([
    0,
    1,
    3,
    1,
    0,
    2,
    5,
    2,
    1,
    0,
    4,
    1,
    0
], dtype=float)

peaks, properties = find_peaks(
    signal,
    height=2.5,
    distance=2
)

print(
    f"피크 위치: {peaks}"
)

print(
    f"피크 높이: "
    f"{properties['peak_heights']}"
)

fig, ax = plt.subplots()

ax.plot(
    signal,
    marker="o"
)

ax.scatter(
    peaks,
    signal[peaks],
    marker="x",
    s=120,
    label="탐지된 피크"
)

ax.set_title(
    "신호 피크 탐지"
)

ax.set_xlabel("샘플 위치")
ax.set_ylabel("신호 값")
ax.legend()
ax.grid(alpha=0.3)

plt.show()

`find_peaks()`는 1차원 신호의 국소 최대값을 찾고 높이, 최소 간격, 돌출도, 너비 등의 조건으로 피크를 선별할 수 있습니다. NaN이 포함된 데이터에서는 예상과 다른 결과가 발생할 수 있어 사전 처리가 필요합니다.

height
→ 피크의 최소·최대 높이

distance
→ 피크 사이의 최소 샘플 거리

prominence
→ 주변 기준선보다 얼마나 두드러지는지

width
→ 피크의 너비

잡음이 심하면 작은 봉우리까지 모두 피크로 인식할 수 있습니다.

이럴 때는 필터링한 뒤 `height`, `distance`, `prominence`를 조절합니다.

피크 탐지는 산봉우리 개수를 세는 일과 비슷합니다.

동네 언덕까지 에베레스트라고 부르지 않도록 기준이 필요합니다. 🏔️

26. FFT로 주파수 분석하기

시간에 따라 변하는 신호를 보면 어떤 주파수가 포함되어 있는지 바로 알기 어렵습니다.

푸리에 변환은 시간 영역 신호를 주파수 영역으로 변환합니다.

시간 영역
→ 언제 얼마나 흔들렸는가?

주파수 영역
→ 어떤 속도의 진동이 얼마나 포함되었는가?

import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import (
    rfft,
    rfftfreq
)


sampling_rate = 200
duration = 2

time = np.linspace(
    0,
    duration,
    sampling_rate * duration,
    endpoint=False
)

signal = (
    np.sin(
        2 * np.pi * 5 * time
    )
    + 0.5
    * np.sin(
        2 * np.pi * 20 * time
    )
)

5Hz와 20Hz 신호가 섞여 있습니다.

spectrum = rfft(
    signal
)

frequencies = rfftfreq(
    signal.size,
    d=1 / sampling_rate
)

amplitude = (
    np.abs(spectrum)
    * 2
    / signal.size
)

실수 입력 신호에서는 `rfft()`를 사용해 음수 주파수의 중복 부분을 제외한 스펙트럼을 계산할 수 있습니다.

fig, ax = plt.subplots()

ax.plot(
    frequencies,
    amplitude
)

ax.set_xlim(
    0,
    50
)

ax.set_title(
    "신호 주파수 스펙트럼"
)

ax.set_xlabel("주파수")
ax.set_ylabel("진폭")
ax.grid(alpha=0.3)

plt.show()

그래프에서는 5Hz와 20Hz 부근에 큰 봉우리가 나타납니다.

시간 그래프에서는 복잡한 파도처럼 보였지만 주파수 그래프에서는 어떤 진동들이 숨어 있었는지 신분증 검사가 가능합니다. 🌊

27. 공간 거리 계산하기

`scipy.spatial`은 좌표, 거리, 최근접 이웃, 공간 분할과 관련된 알고리즘을 제공합니다.

여러 지점 사이의 거리를 한 번에 계산해 보겠습니다.

import numpy as np

warehouses = np.array([
    [0, 0],
    [5, 2],
    [8, 8]
], dtype=float)

customers = np.array([
    [1, 1],
    [6, 4],
    [10, 7]
], dtype=float)

from scipy.spatial.distance import cdist

distances = cdist(
    warehouses,
    customers,
    metric="euclidean"
)

print(distances)

결과의 형태:

행
→ 창고

열
→ 고객

distances[0, 0]
→ 첫 번째 창고와 첫 번째 고객의 거리

distances[2, 1]
→ 세 번째 창고와 두 번째 고객의 거리

`cdist()`는 두 관측 집합에 포함된 모든 점 쌍 사이의 거리를 계산하며 유클리드, 맨해튼, 코사인, 해밍 등 다양한 거리 측정 방식을 지원합니다.

nearest_warehouse = np.argmin(
    distances,
    axis=0
)

nearest_distance = np.min(
    distances,
    axis=0
)

for customer_index, (
    warehouse_index,
    distance
) in enumerate(
    zip(
        nearest_warehouse,
        nearest_distance
    ),
    start=1
):
    print(
        f"{customer_index}번 고객 → "
        f"{warehouse_index + 1}번 창고, "
        f"거리 {distance:.2f}"
    )

거리 계산은 물류뿐 아니라 추천 시스템, 군집화, 이미지 비교, 이상 탐지에도 사용됩니다.

28. 희소배열 이해하기

다음과 같은 거대한 배열이 있다고 가정해 보겠습니다.

[5, 0, 0, 0, 0, 0]
[0, 0, 3, 0, 0, 0]
[0, 0, 0, 0, 7, 0]
[0, 0, 0, 0, 0, 0]

값 대부분이 0입니다.

이런 배열을 희소배열이라고 합니다.

반대로 값이 대부분 채워져 있는 일반 배열은 밀집배열이라고 부릅니다.

사용자와 상품의 구매 여부

문서와 단어의 등장 횟수

대규모 네트워크 연결

그래프 인접 행렬

유한요소 해석 행렬

추천 시스템 평점

0을 모두 저장하면 메모리가 낭비될 수 있습니다.

희소 형식은 주로 0이 아닌 값과 그 위치를 저장합니다.

일반 배열:
모든 칸 저장

희소배열:
값이 있는 칸 중심으로 저장

import numpy as np
from scipy.sparse import csr_array


dense_array = np.array([
    [1, 0, 2],
    [0, 0, 3],
    [4, 5, 0]
])

sparse_array = csr_array(
    dense_array
)

print(sparse_array)

print(
    f"형태: {sparse_array.shape}"
)

print(
    f"저장된 값 개수: "
    f"{sparse_array.nnz}"
)

print(
    f"실제 데이터: "
    f"{sparse_array.data}"
)

restored = sparse_array.toarray()

print(restored)

CSR은 Compressed Sparse Row의 줄임말입니다.

행 단위 연산과 행 슬라이싱, 행렬·벡터 곱셈에 효율적인 희소 형식입니다. SciPy는 기존 희소행렬 인터페이스에서 희소배열 인터페이스로 이동하고 있으므로 신규 코드에서는 `csr_array` 같은 희소배열 API를 우선 검토할 수 있습니다.

0이 별로 없는 배열에서는 희소 형식의 위치 정보 관리 비용 때문에 오히려 비효율적일 수 있습니다.

0이 대부분
→ 희소배열 검토

값이 대부분 존재
→ 일반 NumPy 배열 검토

빈 좌석이 많은 경기장에서는 앉아 있는 사람의 좌석번호만 기록하는 것이 효율적입니다.

모든 좌석이 가득 찼다면 명단 관리가 더 복잡해질 수 있습니다. 🏟️

29. 희소 연립방정식 풀기

다음 연립방정식을 생각해 보겠습니다.

3x = 6
4y = 8
5z = 10

행렬로 표현하면 다음과 같습니다.

A × x = b

import numpy as np
from scipy.sparse import csr_array
from scipy.sparse.linalg import spsolve


coefficient_matrix = csr_array([
    [3.0, 0.0, 0.0],
    [0.0, 4.0, 0.0],
    [0.0, 0.0, 5.0]
])

right_hand_side = np.array([
    6.0,
    8.0,
    10.0
])

solution = spsolve(
    coefficient_matrix,
    right_hand_side
)

print(solution)

결과:

[2. 2. 2.]

`spsolve()`는 희소배열이나 희소행렬로 표현된 `A`에 대해 `Ax=b` 형태의 선형 시스템을 풉니다. 입력 A는 내부적으로 CSR 또는 CSC 형식으로 변환될 수 있습니다.

희소 선형 시스템은 공학 시뮬레이션, 그래프 분석, 편미분방정식 근사 등에서 중요합니다.

30. SciPy 선형대수 활용하기

NumPy에도 `numpy.linalg`이 있지만 SciPy의 `scipy.linalg`은 보다 폭넓은 선형대수 기능을 제공합니다.

import numpy as np
from scipy.linalg import solve


coefficient_matrix = np.array([
    [2.0, 1.0],
    [1.0, -1.0]
])

right_hand_side = np.array([
    5.0,
    1.0
])

solution = solve(
    coefficient_matrix,
    right_hand_side
)

print(solution)

결과:

[2. 1.]

from scipy.linalg import lu

permutation, lower, upper = lu(
    coefficient_matrix
)

print("P:")
print(permutation)

print("L:")
print(lower)

print("U:")
print(upper)

LU 분해는 행렬을 하삼각행렬과 상삼각행렬 등의 곱으로 분해합니다.

반복적으로 연립방정식을 풀거나 행렬 구조를 분석할 때 활용할 수 있습니다.

from scipy.linalg import eig

matrix = np.array([
    [2.0, 0.0],
    [0.0, 3.0]
])

eigenvalues, eigenvectors = eig(
    matrix
)

print(eigenvalues)
print(eigenvectors)

선형대수 계산에서는 수치 안정성이 중요합니다.

역행렬을 직접 구해 곱하기보다 문제에 맞는 `solve()` 계열 함수를 사용하는 편이 일반적으로 더 직접적입니다.

31. 특수 함수 사용하기

SciPy의 `special` 하위 패키지는 감마함수, 베타함수, 오차함수, 베셀함수 등 과학 계산에서 자주 등장하는 특수 함수를 제공합니다.

감마함수는 팩토리얼을 실수와 복소수 영역으로 확장한 함수입니다.

양의 정수 n에 대해 다음 관계가 있습니다.

Γ(n) = (n - 1)!

from scipy.special import gamma

print(
    gamma(1)
)

print(
    gamma(5)
)

print(
    gamma(0.5)
)

결과 해석:

gamma(1)
→ 0! = 1

gamma(5)
→ 4! = 24

gamma(0.5)
→ √π

import numpy as np
from scipy.special import gamma

print(
    np.sqrt(np.pi)
)

print(
    gamma(0.5)
)

특수 함수는 일반 애플리케이션에서 자주 보이지 않지만 확률분포, 물리학, 신호 처리, 미분방정식에서는 갑자기 주연 배우로 등장합니다.

32. 실전 미니 프로젝트: 카페 음료 냉각 분석 시스템

뜨거운 커피가 시간이 지나면서 식는 과정을 분석해 보겠습니다.

이번 프로젝트에서는 다음 기능을 함께 사용합니다.

curve_fit
→ 냉각 모델의 매개변수 추정

root_scalar
→ 목표 온도에 도달하는 시간 계산

quad
→ 일정 시간 동안 주변 온도보다 높은 정도 계산

Matplotlib
→ 측정값과 모델 시각화

카페에서 커피 온도를 2분 간격으로 측정했다고 가정합니다.

초기 온도는 몇 도였는가?

주변 온도는 몇 도로 추정되는가?

냉각 속도는 얼마인가?

마시기 좋은 45℃가 되는 시점은 언제인가?

처음 30분 동안 주변 온도보다 얼마나 뜨거웠는가?

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import quad
from scipy.optimize import (
    curve_fit,
    root_scalar
)


OUTPUT_FILE = Path(
    "coffee_cooling_analysis.png"
)


def cooling_model(
    time: np.ndarray | float,
    room_temperature: float,
    temperature_difference: float,
    cooling_rate: float
) -> np.ndarray | float:
    """
    뉴턴의 냉각 법칙 형태의 모델입니다.

    T(t) = T_room + A * exp(-k * t)
    """
    return (
        room_temperature
        + temperature_difference
        * np.exp(
            -cooling_rate * time
        )
    )


def fit_model(
    times: np.ndarray,
    temperatures: np.ndarray
) -> tuple[
    np.ndarray,
    np.ndarray
]:
    """측정 데이터에 냉각 모델을 맞춥니다."""
    initial_guess = [
        24.0,
        65.0,
        0.08
    ]

    lower_bounds = [
        10.0,
        0.0,
        0.0001
    ]

    upper_bounds = [
        35.0,
        100.0,
        1.0
    ]

    parameters, covariance = curve_fit(
        cooling_model,
        times,
        temperatures,
        p0=initial_guess,
        bounds=(
            lower_bounds,
            upper_bounds
        )
    )

    return parameters, covariance


def find_target_time(
    parameters: np.ndarray,
    target_temperature: float
) -> float:
    """모델 온도가 목표값에 도달하는 시간을 찾습니다."""
    room_temperature, \
        temperature_difference, \
        cooling_rate = parameters

    def target_equation(
        time: float
    ) -> float:
        predicted = cooling_model(
            time,
            room_temperature,
            temperature_difference,
            cooling_rate
        )

        return (
            predicted
            - target_temperature
        )

    result = root_scalar(
        target_equation,
        bracket=(0, 120),
        method="brentq"
    )

    if not result.converged:
        raise RuntimeError(
            "목표 온도 도달 시간을 "
            "계산하지 못했습니다."
        )

    return float(result.root)


def calculate_heat_exposure(
    parameters: np.ndarray,
    end_time: float
) -> tuple[
    float,
    float
]:
    """
    주변 온도를 초과한 온도 차이를
    지정 시간까지 적분합니다.
    """
    room_temperature, \
        temperature_difference, \
        cooling_rate = parameters

    def excess_temperature(
        time: float
    ) -> float:
        predicted = cooling_model(
            time,
            room_temperature,
            temperature_difference,
            cooling_rate
        )

        return (
            predicted
            - room_temperature
        )

    integral, error = quad(
        excess_temperature,
        0,
        end_time
    )

    return float(integral), float(error)


def create_chart(
    times: np.ndarray,
    temperatures: np.ndarray,
    parameters: np.ndarray,
    target_temperature: float,
    target_time: float
) -> None:
    """측정 데이터와 추정 모델을 그래프로 저장합니다."""
    fine_times = np.linspace(
        0,
        40,
        400
    )

    fitted_temperatures = cooling_model(
        fine_times,
        *parameters
    )

    fig, ax = plt.subplots(
        figsize=(11, 6)
    )

    ax.scatter(
        times,
        temperatures,
        label="측정 온도",
        zorder=3
    )

    ax.plot(
        fine_times,
        fitted_temperatures,
        linewidth=2.5,
        label="추정 냉각 모델"
    )

    ax.axhline(
        target_temperature,
        linestyle="--",
        label=(
            f"목표 온도 "
            f"{target_temperature:.0f}℃"
        )
    )

    ax.axvline(
        target_time,
        linestyle=":",
        label=(
            f"도달 시간 "
            f"{target_time:.1f}분"
        )
    )

    ax.scatter(
        [target_time],
        [target_temperature],
        marker="x",
        s=130,
        zorder=4
    )

    ax.set_title(
        "카페 음료 냉각 분석"
    )

    ax.set_xlabel(
        "경과 시간(분)"
    )

    ax.set_ylabel(
        "온도(℃)"
    )

    ax.grid(
        alpha=0.3
    )

    ax.legend()

    fig.savefig(
        OUTPUT_FILE,
        dpi=250,
        bbox_inches="tight"
    )

    plt.show()


def main() -> None:
    times = np.array([
        0,
        2,
        4,
        6,
        8,
        10,
        12,
        15,
        18,
        22,
        26,
        30
    ], dtype=float)

    temperatures = np.array([
        89.5,
        80.7,
        73.2,
        66.8,
        61.1,
        56.4,
        52.1,
        47.4,
        43.5,
        39.5,
        36.6,
        34.4
    ], dtype=float)

    parameters, covariance = fit_model(
        times,
        temperatures
    )

    room_temperature, \
        temperature_difference, \
        cooling_rate = parameters

    estimated_initial_temperature = (
        room_temperature
        + temperature_difference
    )

    target_temperature = 45.0

    target_time = find_target_time(
        parameters,
        target_temperature
    )

    heat_exposure, error = (
        calculate_heat_exposure(
            parameters,
            end_time=30
        )
    )

    print(
        "===== 커피 냉각 분석 ====="
    )

    print(
        "추정 주변 온도: "
        f"{room_temperature:.2f}℃"
    )

    print(
        "추정 초기 온도: "
        f"{estimated_initial_temperature:.2f}℃"
    )

    print(
        "추정 냉각 계수: "
        f"{cooling_rate:.5f}"
    )

    print(
        f"{target_temperature:.0f}℃ "
        "도달 예상 시간: "
        f"{target_time:.2f}분"
    )

    print(
        "30분간 누적 초과 온도: "
        f"{heat_exposure:.2f}℃·분"
    )

    print(
        "적분 추정 오차: "
        f"{error:.3e}"
    )

    standard_errors = np.sqrt(
        np.diag(covariance)
    )

    print(
        "매개변수 표준오차: "
        f"{standard_errors}"
    )

    create_chart(
        times,
        temperatures,
        parameters,
        target_temperature,
        target_time
    )

    print(
        f"그래프 저장 위치: "
        f"{OUTPUT_FILE.resolve()}"
    )


if __name__ == "__main__":
    main()

프로젝트에서 사용한 함수는 다음과 같습니다.

T(t) = T환경 + A × e^(-kt)

T환경
→ 주변 온도

A
→ 초기 온도와 주변 온도의 차이

k
→ 냉각 속도 계수

parameters, covariance = curve_fit(
    cooling_model,
    times,
    temperatures,
    p0=initial_guess,
    bounds=(
        lower_bounds,
        upper_bounds
    )
)

측정값에 가장 잘 맞는 다음 값을 찾습니다.

주변 온도
초기 온도 차이
냉각 계수

predicted_temperature - 45 = 0

이 방정식의 해를 `root_scalar()`로 찾습니다.

result = root_scalar(
    target_equation,
    bracket=(0, 120),
    method="brentq"
)

주변 온도보다 높은 온도 차이를 0분부터 30분까지 적분합니다.

integral, error = quad(
    excess_temperature,
    0,
    30
)

단위는 대략 다음과 같습니다.

℃ × 분
→ ℃·분

이 값은 30분 동안 음료가 주변 온도보다 얼마나 높은 상태를 유지했는지 나타내는 하나의 요약 지표로 볼 수 있습니다.

이 프로젝트는 SciPy 사용법을 설명하기 위한 단순화된 예제입니다.

실제 냉각 과정은 용기 재질, 뚜껑, 대류, 증발, 음료량, 주변 공기 흐름 등 여러 요인의 영향을 받을 수 있습니다.

커피 한 잔을 분석했을 뿐인데 최적화팀, 방정식팀, 적분팀이 모두 출동했습니다. ☕🔬

33. SciPy 사용 시 주의할 점

최적화 결과는 반드시 성공 여부를 확인합니다.

if not result.success:
    raise RuntimeError(
        result.message
    )

방정식 풀이에서는 다음을 확인합니다.

if not result.converged:
    raise RuntimeError(
        "해를 찾지 못했습니다."
    )

수치 알고리즘은 입력값, 초기값, 허용오차, 함수 형태에 영향을 받습니다.

잘못된 모델
+ 정확한 계산
= 정확하게 계산된 잘못된 결론

다음 단위가 섞이면 결과가 크게 잘못될 수 있습니다.

초와 분
미터와 킬로미터
Hz와 kHz
섭씨와 절대온도
원과 천 원

변수명에 단위를 포함하면 실수를 줄일 수 있습니다.

duration_seconds = 30
distance_kilometers = 12.5
sampling_rate_hz = 100

로그 함수
→ 0 이하에서 주의

제곱근
→ 음수 입력 주의

나눗셈
→ 분모 0 주의

보간
→ x 중복값 주의

SciPy는 수렴 실패, 부정확한 적분, 잘못된 희소구조 등에 대해 경고를 출력할 수 있습니다.

경고를 숨기기 전에 원인을 확인해야 합니다.

곡선 맞춤과 보간 결과는 숫자만 보지 말고 그래프로 확인합니다.

모델이 데이터 흐름을 잘 따라가는가?

일부 구간에서 비정상적으로 튀는가?

외삽 결과가 현실적인가?

이상치가 모델을 지배하는가?

통계 함수가 실행되었다고 해서 데이터가 검정 조건을 만족한다는 뜻은 아닙니다.

독립성, 분포, 분산, 표본 추출 방식 등을 별도로 검토해야 합니다.

한 변수 최적화
→ minimize_scalar

여러 변수 최적화
→ minimize

한 변수 방정식
→ root_scalar

정적분
→ quad

초기값 미분방정식
→ solve_ivp

곡선 매개변수 추정
→ curve_fit

SciPy는 도구 상자입니다.

망치가 훌륭하다고 모든 문제를 못으로 보면 연구소 벽에 구멍만 늘어납니다. 🔨

34. 자주 발생하는 오류

ModuleNotFoundError:
No module named 'scipy'

해결:

python -m pip install scipy

현재 실행 환경을 확인합니다.

python -m pip show scipy

`minimize()`의 목적함수는 일반적으로 하나의 스칼라 값을 반환해야 합니다.

잘못된 예:

def objective(x):
    return np.array([
        x[0] ** 2,
        x[1] ** 2
    ])

수정:

def objective(x):
    return (
        x[0] ** 2
        + x[1] ** 2
    )

def objective(values):
    x, y = values

    return x ** 2 + y ** 2

변수가 두 개이므로 초기값도 두 개가 필요합니다.

initial_guess = np.array([
    0.0,
    0.0
])

def function(x):
    return x ** 2 + 1

이 함수는 실수 영역에서 항상 양수입니다.

root_scalar(
    function,
    bracket=(-2, 2)
)

부호가 바뀌지 않으므로 bracket 기반 방법으로 실수 해를 찾을 수 없습니다.

f(-2) = 5
f(2) = 5

먼저 구간 양 끝의 함수값을 확인합니다.

print(function(-2))
print(function(2))

x = np.array([
    0,
    1,
    1,
    2
])

PCHIP 입력 x에는 중복값이 없어야 합니다.

중복 측정이 있다면 평균이나 다른 집계 기준으로 먼저 정리해야 합니다.

x = np.array([
    0,
    3,
    2,
    5
])

정렬 순서를 구합니다.

order = np.argsort(x)

x_sorted = x[order]
y_sorted = y[order]

샘플링 주파수가 100Hz라면 표현 가능한 최고 주파수는 나이퀴스트 주파수인 50Hz입니다.

sampling_rate = 100

다음 차단주파수는 잘못된 설정입니다.

cutoff = 60

`fs`와 차단주파수의 관계를 확인해야 합니다.

`sosfiltfilt()`는 신호 양끝을 확장해 앞뒤 방향으로 필터링합니다.

신호가 너무 짧으면 패딩 길이와 관련된 오류가 발생할 수 있습니다.

더 긴 신호를 사용하거나 필터 설정과 `padlen`을 검토합니다.

peaks, _ = find_peaks(
    noisy_signal
)

잡음의 작은 봉우리도 모두 탐지될 수 있습니다.

조건을 추가합니다.

peaks, properties = find_peaks(
    noisy_signal,
    height=1.0,
    distance=20,
    prominence=0.5
)

필요하면 신호를 먼저 평활화하거나 필터링합니다.

`ttest_ind()`는 `nan_policy`를 지정할 수 있습니다.

result = ttest_ind(
    group_a,
    group_b,
    equal_var=False,
    nan_policy="omit"
)

하지만 무조건 제거하기 전에 NaN이 발생한 이유를 확인해야 합니다.

dense = sparse_array.toarray()

거대한 희소배열을 밀집배열로 변환하면 메모리가 급격히 증가할 수 있습니다.

필요한 계산이 희소 형식을 지원하는지 먼저 확인합니다.

35. 연습 문제

다음 함수의 최솟값을 찾으세요.

f(x) = (x - 7)² + 3

from scipy.optimize import minimize_scalar

다음 함수의 최대값을 찾으세요.

f(x) = -x² + 10x + 5

범위는 0부터 10입니다.

함수에 음수를 붙인 뒤 최소화합니다.

다음 방정식의 양수 해를 구하세요.

x³ - 8 = 0

root_scalar(
    function,
    bracket=(0, 5),
    method="brentq"
)

다음 함수를 0부터 3까지 적분하세요.

f(x) = 2x + 1

quad(
    function,
    0,
    3
)

다음 측정값에서 x=2.5일 때의 값을 PCHIP 보간으로 추정하세요.

x = np.array([
    0,
    1,
    2,
    3,
    4
])

y = np.array([
    0,
    2,
    5,
    9,
    14
])

두 그룹의 평균 차이를 Welch t 검정으로 분석하세요.

group_a = np.array([
    10,
    12,
    11,
    13,
    12
])

group_b = np.array([
    8,
    9,
    10,
    9,
    8
])

ttest_ind(
    group_a,
    group_b,
    equal_var=False
)

초당 100개씩 수집한 신호에서 5Hz 이하 성분만 남기는 Butterworth 저주파 통과 필터를 설계하세요.

sos = butter(
    N=4,
    Wn=5,
    btype="lowpass",
    fs=100,
    output="sos"
)

다음 신호에서 높이가 4 이상인 피크를 찾으세요.

signal = np.array([
    0,
    2,
    5,
    1,
    0,
    4,
    7,
    2,
    0
])

다음 두 지점 집합 사이의 모든 유클리드 거리를 계산하세요.

group_a = np.array([
    [0, 0],
    [2, 2]
])

group_b = np.array([
    [1, 1],
    [5, 5]
])

cdist(
    group_a,
    group_b,
    metric="euclidean"
)

다음 배열을 CSR 희소배열로 변환하고 0이 아닌 저장값 개수를 출력하세요.

array = np.array([
    [1, 0, 0],
    [0, 0, 3],
    [0, 5, 0]
])

sparse = csr_array(array)

print(sparse.nnz)

36. 핵심 요약

python -m pip install scipy

from scipy.optimize import minimize_scalar

result = minimize_scalar(
    objective
)

from scipy.optimize import minimize

result = minimize(
    objective,
    initial_guess
)

from scipy.optimize import root_scalar

result = root_scalar(
    function,
    bracket=(start, end)
)

from scipy.optimize import curve_fit

parameters, covariance = curve_fit(
    model,
    x_data,
    y_data
)

from scipy.integrate import quad

integral, error = quad(
    function,
    start,
    end
)

from scipy.integrate import solve_ivp

solution = solve_ivp(
    model,
    t_span,
    initial_state,
    t_eval=time_points
)

from scipy.interpolate import PchipInterpolator

interpolator = PchipInterpolator(
    x,
    y
)

estimated = interpolator(
    new_x
)

from scipy.stats import ttest_ind

result = ttest_ind(
    group_a,
    group_b,
    equal_var=False
)

from scipy.signal import (
    butter,
    sosfiltfilt
)

sos = butter(
    N=4,
    Wn=5,
    btype="lowpass",
    fs=100,
    output="sos"
)

filtered = sosfiltfilt(
    sos,
    signal
)

from scipy.signal import find_peaks

peaks, properties = find_peaks(
    signal,
    height=1,
    prominence=0.5
)

from scipy.spatial.distance import cdist

distances = cdist(
    group_a,
    group_b
)

from scipy.sparse import csr_array

sparse = csr_array(
    dense_array
)

from scipy.sparse.linalg import spsolve

solution = spsolve(
    sparse_matrix,
    right_hand_side
)

37. 마무리

이번 시간에는 NumPy를 넘어 전문적인 과학 계산 문제를 해결하는 SciPy를 알아보았습니다.

SciPy를 사용하면 다음과 같은 질문에 Python 코드로 답할 수 있습니다.

어떤 값에서 비용이 가장 작아지는가?

함수가 0이 되는 위치는 어디인가?

관측 데이터에 가장 잘 맞는 모델은 무엇인가?

복잡한 함수 아래의 넓이는 얼마인가?

시간에 따라 시스템은 어떻게 변화하는가?

측정하지 않은 중간값은 얼마인가?

두 집단의 차이가 통계적으로 유의한가?

신호의 잡음을 어떻게 줄일 수 있는가?

중요한 피크는 어디에 있는가?

거대한 희소 시스템을 어떻게 계산할 것인가?

NumPy에서는 배열을 이렇게 계산했습니다.

result = array * 2 + 10

SciPy에서는 한 단계 더 나아갑니다.

result = minimize(
    objective,
    initial_guess
)

integral, error = quad(
    function,
    start,
    end
)

solution = solve_ivp(
    model,
    time_range,
    initial_state
)

SciPy의 핵심은 어려운 수학 공식을 외우지 않아도 된다는 뜻이 아닙니다.

어떤 문제를 풀고 있는지, 어떤 가정을 사용했는지, 결과의 오차와 한계가 무엇인지 이해해야 올바르게 사용할 수 있습니다.

코드가 실행되었다.
≠
현실적으로 올바른 답을 얻었다.

SciPy는 정답 자동판매기가 아닙니다.

우리가 문제를 올바르게 정의하면 강력한 수치 알고리즘으로 계산을 도와주는 연구 파트너입니다.

NumPy가 숫자 군단을 정렬했다면, SciPy는 그 군단을 데리고 산의 최저점을 찾고, 곡선의 넓이를 재고, 잡음 속에서 중요한 신호를 찾아냅니다.

이제 복잡한 함수나 센서 데이터가 등장해도 겁먹지 마세요.

SciPy 연구소의 어느 문을 두드려야 할지부터 생각하면 됩니다. 🧪⚙️

[Python 완전정복 시리즈 #28] scikit-learn 완벽 이해하기 | 데이터로 미래를 예측하는 머신러닝 첫걸음

다음 시간에는 NumPy와 SciPy를 기반으로 동작하는 대표적인 머신러닝 라이브러리 scikit-learn을 알아봅니다.

데이터 준비, 학습용·평가용 데이터 분리, 모델 학습, 예측, 정확도 평가, 전처리, 파이프라인까지 컴퓨터가 데이터에서 규칙을 배우는 과정을 실습해 보겠습니다.

#Python #파이썬 #Python강좌 #파이썬기초 #SciPy #사이파이 #NumPy #과학계산 #수치계산 #최적화 #수치적분 #보간법 #통계분석 #가설검정 #신호처리 #FFT #희소행렬 #선형대수 #데이터분석 #머신러닝 #코딩공부 #프로그래밍