SymPy 검산: 손으로 푼 것을 기계로 채점하기

초급 트랙의 마지막 글입니다. 새 수학은 없고, 앞의 47편을 채점하는 도구를 만듭니다.

3번 · 손으로 푸는 습관에서 검산 셋을 세웠습니다 — 숫자 넣기, 극단값 보기, 두 방법으로 세기. 그 셋을 그대로 코드로 옮기는 것이 이 글이고, 새로 배우는 것은 도구 사용법뿐입니다.

채점하는 길은 둘입니다. 기호 계산은 식 자체를 견주고, 수치 검산은 값을 견줍니다. 기호로 되면 기호가 확실하고, 안 되면 수치로 갑니다.

기호 계산 — 식과 식을 견준다

SymPy는 수를 계산하는 대신 식을 기호 그대로 다루는 파이썬 도구입니다. 먼저 문자를 만들어야 식을 적을 수 있습니다.

import sympy as sp

x = sp.symbols('x')
sp.diff(x**2 * sp.sin(x), x)     # 2*x*sin(x) + x**2*cos(x)

$x^2 \sin x$ 를 미분한 결과입니다. 32번 · 도함수 계산 규칙의 곱의 법칙으로 손으로 얻은 것과 같습니다.

그런데 눈으로 견주는 것은 믿을 수 없습니다. 같은 식도 적는 방법이 여럿이라 글자가 달라 보일 수 있습니다. 그래서 두 식을 빼서 0인지 봅니다.

mine = 2*x*sp.sin(x) + x**2*sp.cos(x)          # 손으로 구한 답
sp.simplify(sp.diff(x**2*sp.sin(x), x) - mine) == 0    # True

simplify는 식을 최대한 정리하는 명령입니다. 차가 0으로 정리되면 두 식은 같습니다. 이것이 기호 검산의 관용구이고, 아래 도구들과 늘 짝을 이룹니다.

앞의 글들을 그대로 채점해 봅니다.

k = sp.symbols('k')
sp.integrate(x*sp.exp(x), (x, 0, 1))               # 1
sp.series(sp.exp(x), x, 0, 4)                      # 1 + x + x**2/2 + x**3/6 + O(x**4)
sp.summation(sp.Rational(1, 2)**k, (k, 0, sp.oo))  # 2
sp.factor(x**2 - 5*x + 6)                          # (x - 3)*(x - 2)

셋째 줄이 29번의 첫째항 1, 공비 $\tfrac12$ 인 무한등비급수이고, 손으로 얻은 $\dfrac{1}{1 - 1/2} = 2$ 와 같습니다.

두 답이 갈리면 무엇부터 의심하나

대조가 False 를 뱉었을 때 아무 데나 뒤지면 오래 걸립니다. 자주 틀리는 자리부터 봅니다.

1. 부호. 가장 흔합니다. 빼기를 옮기다 하나 놓치거나, $\cos$ 을 미분하고 마이너스를 안 붙인 자리입니다.
2. 연쇄법칙의 속미분. $\sin(3x)$ 를 미분하고 $\cos(3x)$ 라고 적으면 안쪽 $3x$ 를 미분한 3이 빠진 것입니다. 옳은 답은 $3\cos(3x)$ 입니다.
3. 괄호. 종이에서는 보이던 묶음이 한 줄로 적히며 풀립니다. $1/(2x)$ 를 1/2*x 로 적으면 다른 식입니다.
4. 정의역 조건. 여기까지 와도 안 맞으면 두 식이 어떤 $x$ 에서 같은가를 봅니다. $\dfrac{x}{x}$ 와 $1$ 은 $x = 0$ 을 빼면 같습니다. simplify 는 이런 자리를 조용히 지나가므로, 기호 대조가 True 여도 정의역까지 같다는 뜻은 아닙니다.

수치 검산 — 값과 값을 견준다

기호로 안 풀리거나 결과가 너무 복잡하면 값으로 갑니다. 3번의 숫자 넣기가 여기 그대로 옵니다.

가장 싼 방법은 특수값입니다. 양변에 $x = 0$, $1$, $-1$ 을 넣어 봅니다. 값이 하나라도 갈리면 그 자리에서 틀린 것이 확정됩니다. 다만 모두 맞아도 참이라는 보장은 없습니다 — 세 점에서만 겹치는 다른 식도 있기 때문입니다. 그래서 여러 점에서 재어 그래프로 겹쳐 보는 편이 낫습니다.

도함수는 중심차분으로

31번 · 도함수의 정의는 $h$ 를 한쪽으로만 보냈지만, 검산에는 양쪽을 쓰는 편이 훨씬 정확합니다.

$f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$

import math

def central(f, x, h):
    return (f(x + h) - f(x - h)) / (2 * h)

central(math.sin, 1, 0.01)   # 0.5402933008747335
math.cos(1)                  # 0.5403023058681398

앞의 네 자리 $0.5402$ 까지 같습니다. 47번 · 자릿수와 오차의 상대오차로 재면 $1.7 \times 10^{-5}$ 이라 약 4~5자리가 맞았습니다. $h$ 를 10분의 1인 $0.001$ 로 줄이면 상대오차가 $1.7 \times 10^{-7}$ 로 100분의 1이 됩니다 — 중심차분의 오차가 $O(h^2)$ 이라는 뜻입니다.

정적분은 리만합으로

33번 · 정적분에서 넓이를 직사각형으로 잘라 더했습니다. 그것을 그대로 적습니다. 각 조각의 가운데 높이를 쓰는 것이 중점법입니다.

def riemann(f, a, b, n):
    w = (b - a) / n
    return sum(f(a + (i + 0.5) * w) for i in range(n)) * w

riemann(lambda t: t * math.exp(t), 0, 1, 100)    # 0.999981514413776
riemann(lambda t: t * math.exp(t), 0, 1, 1000)   # 0.9999998151431911

위에서 sp.integrate 가 준 $1$ 과 맞습니다. 조각을 10배로 늘렸더니 오차가 100분의 1이 되었으니 중점법의 오차도 $O(h^2)$ 입니다.

확률은 난수로

39번 · 세는 법의 생일 문제입니다. 서른 명의 생일을 실제로 10만 번 뽑아 봅니다.

import random

hit = 0
for _ in range(100000):
    days = [random.randrange(365) for _ in range(30)]
    if len(set(days)) < 30:
        hit += 1
hit / 100000     # 0.70622

손계산 값은 $0.70632$ 였습니다. 상대오차가 $1.7 \times 10^{-4}$ 이라 약 3~4자리가 맞았습니다.

38번 · 조건부확률의 검사 문제도 같은 방법으로 잽니다. 10만 명을 만들어 양성인 사람만 세면 $0.171$ 이 나오고, 베이즈 정리로 얻은 $1/6 \approx 0.1667$ 과 앞의 두 자리가 맞습니다.

난수 검산이 사 주는 자릿수에는 한계가 있습니다. $N$ 번 돌려 얻는 상대오차는 대략 $1/\sqrt{N}$ 이라, 10만 번이면 $0.3$ % 남짓 — 두세 자리입니다. 자릿수를 하나 더 얻으려면 회수를 100배로 늘려야 합니다.

그래서 안 되는 것도 있습니다. 39번의 로또 1등 확률은 $1 / 8{,}145{,}060 \approx 1.2 \times 10^{-7}$ 입니다. 10만 번을 뽑으면 1등이 나올 기댓값이 $0.012$ 회라 거의 언제나 0회가 나오고, 추정값은 0이 됩니다. 대신 확률이 큰 사건으로 옮겨 재면 됩니다 — 6개 중 3개가 맞을 확률은 손계산으로

$\frac{C(6,3)\,C(39,3)}{C(45,6)} = \frac{20 \times 9139}{8145060} \approx 0.02244$

이고, 10만 번을 뽑으니 $0.02305$ 가 나와 두 자리가 맞았습니다. 재려는 확률의 역수보다 훨씬 많이 뽑아야 잴 수 있다 — 이것이 난수 검산의 조건입니다.

Σ를 코드로 옮길 때의 습관 하나

27번 · 시그마 표기의 이중합을 코드로 옮기면 첨자를 놓치기 쉽습니다. 같은 합을 두 축 순서로 더해 값이 같은지 봅니다.

a = [[1, 2, 3], [4, 5, 6]]
by_row = sum(sum(row) for row in a)                       # 21
by_col = sum(sum(row[j] for row in a) for j in range(3))  # 21

3번의 두 방법으로 세기를 코드로 옮긴 것입니다. 값이 갈리면 첨자를 잘못 옮긴 것입니다. 어떤 조건에서 합의 순서를 바꿔도 되는가라는 규칙 자체는 중급 3번의 몫이고, 여기서는 두 값이 같은지 보는 데까지입니다.

연습 문제

연습 1 — 기호 검산 읽기

1. 두 식을 sp.simplify(a - b) == 0 으로 대조해 True 를 얻었는데도 두 식이 같지 않을 수 있습니다. 어떤 경우인지 예를 하나 들어 답하세요.

   정의역이 다를 때입니다. $\dfrac{x}{x}$ 와 $1$ 은 차가 0으로 정리되지만 $x = 0$ 에서는 앞의 식이 정의되지 않습니다. 대조가 True 라도 "어떤 $x$ 에서 같은가"는 따로 봐야 합니다.
2. 손으로 $\dfrac{d}{dx}\sin(3x) = \cos(3x)$ 를 얻어 대조했더니 어긋났습니다. 의심 순서의 몇 번째 항목인지 답하고 옳은 답을 적으세요.

   두 번째, 연쇄법칙의 속미분입니다. 안쪽 $3x$ 를 미분한 3이 빠졌으므로 옳은 답은 $3\cos(3x)$ 입니다.
3. 손으로 $\displaystyle\int_0^1 6x^2\,dx = 6$ 을 얻어 sp.integrate 와 대조했더니 어긋났습니다. 옳은 값을 구하고 어디가 빠졌는지 적으세요.

   $\displaystyle\int_0^1 6x^2\,dx = \left[2x^3\right]_0^1 = 2$ 입니다. $x^2$ 의 부정적분이 $\dfrac{x^3}{3}$ 인데 3으로 나누는 것을 빠뜨렸습니다.

연습 2 — 수치로 확인하기

1. 중심차분으로 $f(x) = x^3$ 의 $x = 2$ 에서의 도함수를 $h = 0.1$ 로 손계산하고, 참값과의 절대오차를 구하세요.

   $\dfrac{2.1^3 - 1.9^3}{0.2} = \dfrac{9.261 - 6.859}{0.2} = \dfrac{2.402}{0.2} = 12.01$ 입니다. 참값은 $3x^2 = 12$ 이므로 절대오차가 $0.01$ 입니다.
2. 1번에서 $h$ 를 $0.05$ 로 줄이면 오차가 대략 얼마가 될지 47번의 차수 표기로 근거를 적고 예상 값을 답하세요.

   중심차분의 오차는 $O(h^2)$ 이므로 $h$ 를 반으로 줄이면 오차는 넷으로 나뉘어 약 $0.0025$ 입니다. 실제로 $\dfrac{2.05^3 - 1.95^3}{0.1} = \dfrac{1.20025}{0.1} = 12.0025$ 로 오차가 정확히 $0.0025$ 입니다.
3. 중점 리만합으로 $\displaystyle\int_0^1 x\,dx$ 를 $n = 2$ 로 손계산하고 참값과 견주세요.

   폭이 $0.5$ 이고 두 조각의 가운데는 $0.25$ 와 $0.75$ 이므로 $(0.25 + 0.75) \times 0.5 = 0.5$ 입니다. 참값 $\dfrac12$ 와 정확히 같습니다 — 일차함수에서는 가운데 높이가 조각의 평균 높이와 같기 때문입니다.

연습 3 — 난수 검산과 이중합

1. 난수 10만 번으로 잰 값이 $0.3327$ 이고 손계산이 $0.33333$ 이었습니다. 상대오차를 구하고 몇 자리가 맞았는지 답하세요.

   $\dfrac{|0.3327 - 0.33333|}{0.33333} = \dfrac{0.00063}{0.33333} \approx 1.9 \times 10^{-3}$ 입니다. $-\log_{10}(1.9 \times 10^{-3}) \approx 2.7$ 이므로 약 2~3자리가 맞았습니다.
2. 확률이 $1.2 \times 10^{-7}$ 인 사건을 난수 10만 번으로 재려 하면 어떤 값이 나오는지와 그 까닭을 적으세요.

   10만 번에서 그 사건이 나올 기댓값이 $100000 \times 1.2 \times 10^{-7} = 0.012$ 회라 거의 언제나 0회가 나오고 추정값이 0이 됩니다. 재려면 뽑는 횟수가 확률의 역수인 약 800만보다 훨씬 커야 합니다.
3. 이중합을 코드로 옮겨 놓고 첨자를 제대로 옮겼는지 코드만으로 한 번 더 확인하는 방법을 적으세요.

   더하는 축의 순서를 바꿔 행부터 더한 값과 열부터 더한 값을 견줍니다. 두 값이 갈리면 첨자를 잘못 옮긴 것입니다.

정리

• 채점하는 길은 둘입니다. 기호 계산은 식을, 수치 검산은 값을 견줍니다.
• 기호 검산의 관용구는 sp.simplify(a - b) == 0 입니다. 눈으로 식을 견주지 않습니다.
• 두 답이 갈리면 부호 → 속미분 → 괄호 → 정의역 순으로 의심합니다.
• 중심차분과 중점 리만합은 오차가 $O(h^2)$ 이라, 조각을 10배로 늘리면 오차가 100분의 1이 됩니다.
• 난수 검산이 사 주는 상대오차는 대략 $1/\sqrt{N}$ 입니다. 10만 번이면 두세 자리이고, 재려는 확률의 역수보다 훨씬 많이 뽑아야 합니다.
• 이중합은 두 축 순서로 더해 값을 견줍니다. 3번의 "두 방법으로 세기"를 코드로 옮긴 것입니다.

여기까지가 초급 48편입니다. 이제 앞의 47편 어느 자리로 돌아가도 손으로 낸 답을 스스로 채점할 수 있습니다.

읽어주셔서 감사합니다. 😊