서술에서 스물넷을 옮겨 적었다. 이 노트북이 하는 일은 셋이다.
검산 — 세운 , , 가 정말 같은 이야기인지 잰다.
흔들기 — 계수를 밀어 답이 어디서 무너지는지 본다.
부수기 — 일부러 틀린 번역을 실행해 결과가 얼마나 어긋나는지 본다.
세 번째가 제일 중요하다. 전치를 빼먹으면 얼마나 틀리는지를 숫자로 봐야 다음에 안 빼먹는다.
import itertools
import numpy as np
np.set_printoptions(precision=4, suppress=True, linewidth=110)
def 세웠다(번호, A, x, b, 말=""):
A, x, b = map(np.asarray, (A, x, b))
잔차 = np.abs(A @ x - b).max()
맞음 = 잔차 < 1e-9
print(f"[{번호}] {'맞다 ' if 맞음 else '틀리다'} "
f"|Ax-b|_max = {잔차:.2e} {말}")
def 같다(번호, 얻은, 적힌, 말="", 허용=1e-9):
얻은, 적힌 = np.asarray(얻은, dtype=float), np.asarray(적힌, dtype=float)
차 = np.abs(얻은 - 적힌).max()
print(f"[{번호}] {'맞다 ' if 차 < 허용 else '틀리다'} "
f"차이 {차:.2e} {말}")A = np.array([[0.5, 0.1, 0.3],
[0.2, 0.5, 0.1],
[0.1, 0.2, 0.4]])
IA = np.eye(3) - A
print("열 합 :", A.sum(axis=0), " -> 전부 0.8 이면 부가가치가 0.2")
print(f"det(I-A) = {np.linalg.det(IA):.6f}")
print()
# 표를 옳게 옮겼다면 원래 총산출이 복원돼야 한다
x = np.linalg.solve(IA, [100, 50, 50])
세웠다("1", IA, x, [100, 50, 50])
같다("1", x, [420, 320, 260], "총산출이 복원되면 표를 옳게 읽은 것이다")
x2 = np.linalg.solve(IA, [110, 50, 50])
print()
print(f"농림 주문 +10 -> x' = {x2}")
print(f"증가분 = {x2 - x} <- 10 이 아니라 28 이다")열 합 : [0.8 0.8 0.8] -> 전부 0.8 이면 부가가치가 0.2
det(I-A) = 0.100000
[1] 맞다 |Ax-b|_max = 2.13e-14
[1] 맞다 차이 5.68e-14 총산출이 복원되면 표를 옳게 읽은 것이다
농림 주문 +10 -> x' = [448. 333. 269.]
증가분 = [28. 13. 9.] <- 10 이 아니라 28 이다
농림은 10이 아니라 28을 더 만들어야 한다. 그리고 아무도 주문하지 않은 제조가 13, 서비스가 9를 함께 더 만든다.
하나 = np.ones((3, 3))
A2 = 4 * np.eye(3) - 하나 # 문제 2
p = np.linalg.solve(A2, [12, 16, 12])
세웠다("2", A2, p, [12, 16, 12])
같다("2", p, [13, 14, 13], "1300원, 1400원, 1300원")
같다("2", np.linalg.inv(A2), 0.25 * (np.eye(3) + 하나),
"닫힌 꼴 A^-1 = (I + 11^T)/4")
M쿠르노 = np.eye(3) + 하나 # 문제 3
q = np.linalg.solve(M쿠르노, [90, 80, 70])
세웠다("3", M쿠르노, q, [90, 80, 70])
같다("3", q, [30, 20, 10], "쿠르노 균형 생산량")
P = 100 - q.sum()
같다("3", q, P - np.array([10, 20, 30]), f"q_i = P - c_i 가 맞는가 (P={P:.0f})")
print()
print("두 행렬이 서로 물려 있다 : A2 @ M쿠르노 =")
print(A2 @ M쿠르노)[2] 맞다 |Ax-b|_max = 7.11e-15
[2] 맞다 차이 1.78e-15 1300원, 1400원, 1300원
[2] 맞다 차이 0.00e+00 닫힌 꼴 A^-1 = (I + 11^T)/4
[3] 맞다 |Ax-b|_max = 0.00e+00
[3] 맞다 차이 0.00e+00 쿠르노 균형 생산량
[3] 맞다 차이 0.00e+00 q_i = P - c_i 가 맞는가 (P=40)
두 행렬이 서로 물려 있다 : A2 @ M쿠르노 =
[[4. 0. 0.]
[0. 4. 0.]
[0. 0. 4.]]
다. 균형가격을 정하는 행렬과 균형수량을 정하는 행렬이 서로 역행렬 관계라는 것이 우연이 아니다. 둘 다 같은 랭크 1 구조에서 나왔다.
이번엔 문제 5의 상태가격. 같은 표를 세로로 한 번, 가로로 한 번 읽는다.
D = np.array([[1.0, 1, 1], [3, 2, 1], [1, 2, 4]])
pv = np.array([0.95, 1.95, 2.10])
q = np.linalg.solve(D, pv) # 세로 : 경우별 값
세웠다("5", D, q, pv)
같다("5", q, [0.30, 0.40, 0.25], "호황/보통/불황 1만원의 오늘 값")
print(f" det D = {np.linalg.det(D):.0f}, D^-1 이 정수인가 : "
f"{np.allclose(np.linalg.inv(D), np.round(np.linalg.inv(D)))}")
print(f" 상태가격의 합 {q.sum():.2f} = 무위험 할인계수, "
f"무위험이자율 {1/q.sum()-1:.4%}")
print(f" 전부 양수인가 : {bool((q > 0).all())} <- 아니면 차익거래가 있다")
theta = np.linalg.solve(D.T, [1, 0, 0]) # 가로 : 복제 방법
세웠다("5", D.T, theta, [1, 0, 0], "같은 표를 전치해 읽으면 복제 포트폴리오")
같다("5", theta, [-6, 2, 1], "국채 6주 팔고 우량주 2주, 성장주 1주")
같다("5", theta @ pv, q[0], "복제 비용이 상태가격과 같다")[5] 맞다 |Ax-b|_max = 0.00e+00
[5] 맞다 차이 1.11e-16 호황/보통/불황 1만원의 오늘 값
det D = -1, D^-1 이 정수인가 : True
상태가격의 합 0.95 = 무위험 할인계수, 무위험이자율 5.2632%
전부 양수인가 : True <- 아니면 차익거래가 있다
[5] 맞다 |Ax-b|_max = 0.00e+00 같은 표를 전치해 읽으면 복제 포트폴리오
[5] 맞다 차이 0.00e+00 국채 6주 팔고 우량주 2주, 성장주 1주
[5] 맞다 차이 8.33e-16 복제 비용이 상태가격과 같다
문제 6의 부트스트랩은 자료가 이미 삼각꼴로 도착한 경우다. 소거를 할 필요가 없다.
C = np.array([[100.0, 0, 0], [5, 105, 0], [6, 6, 106]])
P = np.array([95.0, 99.25, 101.20])
print("C 가 하삼각인가 :", np.allclose(C, np.tril(C)))
d = np.zeros(3) # 전진대입을 손으로
for i in range(3):
d[i] = (P[i] - C[i, :i] @ d[:i]) / C[i, i]
같다("6", d, [0.95, 0.90, 0.85], "할인계수")
같다("6", d, np.linalg.solve(C, P), "손으로 한 전진대입 == solve")
for k, dk in enumerate(d, 1):
print(f" z_{k} = {dk ** (-1 / k) - 1:.4%}")C 가 하삼각인가 : True
[6] 맞다 차이 1.11e-16 할인계수
[6] 맞다 차이 1.11e-16 손으로 한 전진대입 == solve
z_1 = 5.2632%
z_2 = 5.4093%
z_3 = 5.5667%
문제 7의 부품표. 순환이 없으면 이 멱영이고, 그래서 무한급수가 저절로 멈춘다.
N = np.array([[0, 1, 2, 0],
[0, 0, 0, 4],
[0, 0, 0, 2],
[0, 0, 0, 0]], dtype=float)
이름 = ["자전거", "프레임조립", "바퀴", "볼트"]
for k in (1, 2, 3):
print(f"N^{k} 이 0 인가 : {np.allclose(np.linalg.matrix_power(N, k), 0)}")
T = np.eye(4) + N + N @ N
같다("7", T[0], [1, 1, 2, 8], "자전거 한 대의 전체 부품표")
같다("7", T, np.linalg.inv(np.eye(4) - N), "I+N+N^2 == (I-N)^-1")
print(f" det(I - N) = {np.linalg.det(np.eye(4) - N):.0f} "
"<- 삼각이고 대각이 1 이라 언제나 1")
주문 = np.array([10.0, 0, 5, 0])
같다("7", 주문 @ T, [10, 10, 25, 90], "자전거 10대 + 바퀴 5개")N^1 이 0 인가 : False
N^2 이 0 인가 : False
N^3 이 0 인가 : True
[7] 맞다 차이 0.00e+00 자전거 한 대의 전체 부품표
[7] 맞다 차이 0.00e+00 I+N+N^2 == (I-N)^-1
det(I - N) = 1 <- 삼각이고 대각이 1 이라 언제나 1
[7] 맞다 차이 0.00e+00 자전거 10대 + 바퀴 5개
문제 8의 재고. 차분의 역이 누적이라는 것을 행렬이 그대로 말한다.
Dm = np.eye(4) - np.eye(4, k=-1)
S = np.linalg.inv(Dm)
같다("8", S, np.tril(np.ones((4, 4))), "D^-1 은 누적합 행렬")
수요 = np.array([20.0, 30, 25, 25])
I0 = 8.0
생산 = (수요.sum() - I0) / 4 # 기말 재고 0
같다("8", 생산, 23.0, "평준화 생산량")
재고 = S @ (생산 - 수요) + I0
같다("8", 재고, [11, 4, 2, 0], "주말 재고")
# 점화식으로 따라가도 같은가
따라 = []
남 = I0
for t in range(4):
남 = 남 + 생산 - 수요[t]
따라.append(남)
같다("8", 따라, 재고, "한 주씩 따라간 것과 행렬로 한 것이 같다")[8] 맞다 차이 0.00e+00 D^-1 은 누적합 행렬
[8] 맞다 차이 0.00e+00 평준화 생산량
[8] 맞다 차이 0.00e+00 주말 재고
[8] 맞다 차이 0.00e+00 한 주씩 따라간 것과 행렬로 한 것이 같다
문제 9는 상자 2다. 만으로는 안 되고 "합이 1"이라는 줄을 더해야 답이 하나로 정해진다.
G = np.array([[-2.0, 4, 0, 0],
[2, -6, 4, 0],
[0, 2, -6, 4],
[0, 0, 2, -4]])
print("열 합 :", G.sum(axis=0), " -> 전부 0 이라 G 는 특이하다")
print(f"rank G = {np.linalg.matrix_rank(G)} / 4 -> 영공간이 1차원")
# 합이 1 이라는 줄을 마지막 줄에 덮어쓴다
증강 = np.vstack([G[:-1], np.ones(4)])
우변 = np.array([0.0, 0, 0, 1])
p = np.linalg.solve(증강, 우변)
같다("9", p, np.array([8, 4, 2, 1]) / 15, "정상분포")
print(f" 거절률 = p_3 = {p[3]:.4f} = 1/15 = {1/15:.4f}")
print(f" 평균 재고 = {np.arange(4) @ p:.4f} = 11/15")
print()
print("선반을 늘리면 :")
for 칸 in range(3, 8):
비 = 0.5 ** np.arange(칸 + 1)
print(f" 선반 {칸}칸 -> 거절률 {비[-1]/비.sum():.4%}")열 합 : [0. 0. 0. 0.] -> 전부 0 이라 G 는 특이하다
rank G = 3 / 4 -> 영공간이 1차원
[9] 맞다 차이 1.39e-17 정상분포
거절률 = p_3 = 0.0667 = 1/15 = 0.0667
평균 재고 = 0.7333 = 11/15
선반을 늘리면 :
선반 3칸 -> 거절률 6.6667%
선반 4칸 -> 거절률 3.2258%
선반 5칸 -> 거절률 1.5873%
선반 6칸 -> 거절률 0.7874%
선반 7칸 -> 거절률 0.3922%
선반 하나마다 거절률이 반씩 준다. 공비 이 도착률과 처리율의 비이기 때문이다.
문제 10의 오분류 되돌리기. 여기서 진짜 물음은 "관측이 틀렸으면 답이 얼마나 흔들리나"였다.
M오분류 = np.array([[0.9, 0.1, 0], [0.1, 0.8, 0.1], [0, 0.1, 0.9]])
b = np.array([0.32, 0.45, 0.23])
x = np.linalg.solve(M오분류, b)
세웠다("10", M오분류, x, b)
같다("10", x, [0.30, 0.50, 0.20], "진짜 등급 비율")
w = np.linalg.eigvalsh(M오분류)
같다("10", np.sort(w), [0.7, 0.9, 1.0], "고윳값")
print(f" 증폭 배율 = 1/lambda_min = {1/w.min():.4f}")[10] 맞다 |Ax-b|_max = 0.00e+00
[10] 맞다 차이 5.55e-17 진짜 등급 비율
[10] 맞다 차이 0.00e+00 고윳값
증폭 배율 = 1/lambda_min = 1.4286
d = np.array([100.0, 50, 50])
참 = np.linalg.solve(np.eye(3) - A, d)
부분 = np.zeros(3)
for k in range(160):
부분 = 부분 + np.linalg.matrix_power(A, k) @ d
if k in (2, 5, 10, 20, 159):
남 = 1 - 부분.sum() / 참.sum()
print(f" {k+1:2d} 라운드 : 누적 {부분.sum():7.2f} / {참.sum():.0f}"
f" 남은 몫 {남:.6f} 0.8^{k+1} = {0.8**(k+1):.6f}")
같다("11", 부분, 참, "급수를 160 항까지 더하면 (I-A)^-1 d 와 같다")
w = np.linalg.eigvals(A)
print()
print(f"고윳값 {np.round(w, 4)}")
print(f"rho(A) = {max(abs(w)):.4f} <- 열 합 0.8 이 그대로 고윳값이다")
print(f"1^T A = {np.ones(3) @ A} (= 0.8 * 1^T)")
print(f"총승수 1/(1-0.8) = {1/(1-0.8):.1f}, 실제 {참.sum()/d.sum():.1f}") 3 라운드 : 누적 488.00 / 1000 남은 몫 0.512000 0.8^3 = 0.512000
6 라운드 : 누적 737.86 / 1000 남은 몫 0.262144 0.8^6 = 0.262144
11 라운드 : 누적 914.10 / 1000 남은 몫 0.085899 0.8^11 = 0.085899
21 라운드 : 누적 990.78 / 1000 남은 몫 0.009223 0.8^21 = 0.009223
160 라운드 : 누적 1000.00 / 1000 남은 몫 -0.000000 0.8^160 = 0.000000
[11] 맞다 차이 1.71e-13 급수를 160 항까지 더하면 (I-A)^-1 d 와 같다
고윳값 [0.8+0.j 0.3+0.1j 0.3-0.1j]
rho(A) = 0.8000 <- 열 합 0.8 이 그대로 고윳값이다
1^T A = [0.8 0.8 0.8] (= 0.8 * 1^T)
총승수 1/(1-0.8) = 5.0, 실제 5.0
문제 12의 기본행렬. 서술에서 강조한 것은 이것이 레온티예프 역행렬과 같은 꼴이라는 점이었다.
Q = np.array([[0.0, 0, 0], [1, 0, 1], [0, 0.2, 0]])
Nb = np.linalg.inv(np.eye(3) - Q)
print(f"det(I-Q) = {np.linalg.det(np.eye(3)-Q):.4f}")
같다("12", Nb[:, 0], [1, 1.25, 0.25], "가공에서 출발했을 때 거치는 평균 횟수")
같다("12", Nb[1, 0], 1 / (1 - 0.2), "검사 횟수 = 1/(1-0.2)")
R = np.array([[0.0, 0.7, 0], [0, 0.1, 0]])
같다("12", (R @ Nb)[:, 0], [0.875, 0.125], "양품 7/8, 폐기 1/8")
print(f" 하루 1000개 -> 검사 {1000*Nb[1,0]:.0f}회, "
f"재작업 {1000*Nb[2,0]:.0f}회, 양품 {1000*(R@Nb)[0,0]:.0f}개")
print()
print("같은 꼴인가 : (I-A)^-1 과 (I-Q)^-1 은 둘 다 I + M + M^2 + ... 이다")
급수Q = sum(np.linalg.matrix_power(Q, k) for k in range(60))
같다("12", 급수Q, Nb, "Q 의 급수도 역행렬로 수렴한다")det(I-Q) = 0.8000
[12] 맞다 차이 0.00e+00 가공에서 출발했을 때 거치는 평균 횟수
[12] 맞다 차이 0.00e+00 검사 횟수 = 1/(1-0.2)
[12] 맞다 차이 0.00e+00 양품 7/8, 폐기 1/8
하루 1000개 -> 검사 1250회, 재작업 250회, 양품 875개
같은 꼴인가 : (I-A)^-1 과 (I-Q)^-1 은 둘 다 I + M + M^2 + ... 이다
[12] 맞다 차이 1.11e-16 Q 의 급수도 역행렬로 수렴한다
문제 13의 펌프. "두 대니까 두 배"가 왜 틀렸는지를 숫자로 본다.
QT = np.array([[-2.0, 0], [2, -1]])
m = np.linalg.solve(-QT.T, [1, 1])
같다("13", m, [1.5, 1.0], "평균 잔여 시간 (2대 출발, 1대 출발)")
print(" 한 대짜리 1시간의 2배가 아니라 1.5배다")
print(f" 2대 구간 평균 {1/2:.2f}시간 + 1대 구간 평균 {1/1:.2f}시간 = {1/2+1:.2f}")
QT2 = np.array([[-2.0, 2], [2, -3]]) # 수리공 한 명
m2 = np.linalg.solve(-QT2.T, [1, 1])
같다("13", m2, [2.5, 2.0], "수리공을 붙이면")
print(f" 수리공 한 명이 버티는 시간을 {m2[0]/m[0]:.2f}배로 만든다")[13] 맞다 차이 0.00e+00 평균 잔여 시간 (2대 출발, 1대 출발)
한 대짜리 1시간의 2배가 아니라 1.5배다
2대 구간 평균 0.50시간 + 1대 구간 평균 1.00시간 = 1.50
[13] 맞다 차이 0.00e+00 수리공을 붙이면
수리공 한 명이 버티는 시간을 1.67배로 만든다
Ad = np.array([[0.8, 0.1, 1.0], [0.2, 0.7, 0.0], [0.0, 0.2, 0.0]])
print("열 합 :", Ad.sum(axis=0))
w, V = np.linalg.eig(Ad)
print(f"고윳값 {np.round(np.real(np.sort(w)), 4)} <- 이산이라 1 을 찾는다")
k = np.argmin(abs(w - 1))
pi = np.real(V[:, k]); pi = pi / pi.sum()
같다("14", pi, np.array([15, 10, 2]) / 27, "정상분포")
print(f" 고장으로 서는 비율 {pi[2]:.4f} = 2/27, 가동률 {1-pi[2]:.4f}")
print()
lam2 = sorted(abs(w))[-2]
print(f"lambda_2 = {lam2:.4f} -> 기억이 이 비율로 죽는다")
for e in (2, 4, 8, 16):
v = np.linalg.matrix_power(Ad, e)[:, 0]
print(f" A^{e:<2d} 첫 열 {np.round(v,4)} pi 와의 거리 "
f"{np.abs(v-pi).max():.6f} 0.4^{e} = {0.4**e:.6f}")열 합 : [1. 1. 1.]
고윳값 [0.1 0.4 1. ] <- 이산이라 1 을 찾는다
[14] 맞다 차이 2.22e-16 정상분포
고장으로 서는 비율 0.0741 = 2/27, 가동률 0.9259
lambda_2 = 0.4000 -> 기억이 이 비율로 죽는다
A^2 첫 열 [0.66 0.3 0.04] pi 와의 거리 0.104444 0.4^2 = 0.160000
A^4 첫 열 [0.5726 0.359 0.0684] pi 와의 거리 0.017044 0.4^4 = 0.025600
A^8 첫 열 [0.556 0.3701 0.0739] pi 와의 거리 0.000437 0.4^8 = 0.000655
A^16 첫 열 [0.5556 0.3704 0.0741] pi 와의 거리 0.000000 0.4^16 = 0.000000
from scipy.linalg import expm
G2 = np.array([[-2.0, 2, 0], [2, -3, 2], [0, 1, -2]])
print("열 합 :", G2.sum(axis=0))
wc = np.linalg.eigvals(G2)
print(f"고윳값 {np.round(np.sort(np.real(wc)), 4)} <- 연속이라 0 을 찾는다")
Vc = np.linalg.eig(G2)[1]
kc = np.argmin(abs(wc))
pic = np.real(Vc[:, kc]); pic = pic / pic.sum()
같다("15", pic, [0.4, 0.4, 0.2], "정상분포, 장기 가용도 80%")
def 가용도(t):
return 1 - (expm(G2 * t) @ np.array([1.0, 0, 0]))[2]
print()
print("두 대 다 멀쩡한 상태에서 출발하면 :")
for t in (0, 0.25, 0.5, 1, 2, 5):
공식 = 4/5 + (1/3)*np.exp(-2*t) - (2/15)*np.exp(-5*t)
print(f" t={t:4.2f}h A(t) = {가용도(t):.6f} 닫힌 꼴 {공식:.6f}")
print(" 30분 뒤는 91.2% 이지 80% 가 아니다")열 합 : [0. 0. 0.]
고윳값 [-5. -2. 0.] <- 연속이라 0 을 찾는다
[15] 맞다 차이 1.67e-16 정상분포, 장기 가용도 80%
두 대 다 멀쩡한 상태에서 출발하면 :
t=0.00h A(t) = 1.000000 닫힌 꼴 1.000000
t=0.25h A(t) = 0.963976 닫힌 꼴 0.963976
t=0.50h A(t) = 0.911682 닫힌 꼴 0.911682
t=1.00h A(t) = 0.844213 닫힌 꼴 0.844213
t=2.00h A(t) = 0.806099 닫힌 꼴 0.806099
t=5.00h A(t) = 0.800015 닫힌 꼴 0.800015
30분 뒤는 91.2% 이지 80% 가 아니다
T = np.array([[0.8, 0.12, 0.08],
[0.2, 0.72, 0.08],
[0.2, 0.12, 0.68]])
print("부서가 준 표 T 의 행 합 :", T.sum(axis=1), " <- 행이 1 이다")
print(" T 의 열 합 :", T.sum(axis=0), " <- 열은 1 이 아니다")
Mo = T.T # 옳게 : 전치한다
Mx = T # 틀리게 : 그대로 쓴다
def 정상(M):
w, V = np.linalg.eig(M)
k = np.argmin(abs(w - 1))
v = np.real(V[:, k])
return v / v.sum()
옳 = 정상(Mo)
틀 = 정상(Mx)
같다("16", 옳, [0.5, 0.3, 0.2], "전치하고 푼 정상 점유율")
print(f"[16] 전치를 빼먹으면 : {np.round(틀, 4)} <- 균등분포가 나온다")
print(f" A사 점유율이 {옳[0]:.1%} 가 아니라 {틀[0]:.1%} 로 나온다. "
f"{abs(옳[0]-틀[0])*100:.1f}%p 틀린다.")부서가 준 표 T 의 행 합 : [1. 1. 1.] <- 행이 1 이다
T 의 열 합 : [1.2 0.96 0.84] <- 열은 1 이 아니다
[16] 맞다 차이 2.22e-16 전치하고 푼 정상 점유율
[16] 전치를 빼먹으면 : [0.3333 0.3333 0.3333] <- 균등분포가 나온다
A사 점유율이 50.0% 가 아니라 33.3% 로 나온다. 16.7%p 틀린다.
전치를 빼먹으면 정상 점유율이 균등분포로 나온다. 행 합이 1인 행렬을 그대로 쓰면 이 오른쪽 고유벡터가 되기 때문이다.
이 실수가 무서운 이유는 답이 그럴듯해 보인다는 데 있다. 셋이 0.333씩이면 "아 점유율이 비슷해지는구나"라고 읽고 넘어가기 쉽다.
M16 = Mo
pi16 = np.array([0.5, 0.3, 0.2])
같다("16", M16, 0.4 * np.outer(pi16, np.ones(3)) + 0.6 * np.eye(3),
"M = 0.4 pi 1^T + 0.6 I 라는 랭크1 구조")
w16 = np.linalg.eigvals(M16)
print(f" 고윳값 {np.round(np.sort(np.real(w16))[::-1], 4)} <- 1, 0.6, 0.6")
x0 = np.array([1.0, 0, 0])
print()
print(f"{'t':>3} {'실제로 곱한 것':>28} {'닫힌 꼴':>28}")
x = x0.copy()
for t in range(1, 4):
x = M16 @ x
닫 = pi16 + 0.6 ** t * (x0 - pi16)
print(f"{t:>3} {str(np.round(x,4)):>28} {str(np.round(닫,4)):>28}")[16] 맞다 차이 1.11e-16 M = 0.4 pi 1^T + 0.6 I 라는 랭크1 구조
고윳값 [1. 0.6 0.6] <- 1, 0.6, 0.6
t 실제로 곱한 것 닫힌 꼴
1 [0.8 0.12 0.08] [0.8 0.12 0.08]
2 [0.68 0.192 0.128] [0.68 0.192 0.128]
3 [0.608 0.2352 0.1568] [0.608 0.2352 0.1568]
a = np.array([2.0, -1, 0, 1, -2])
b = np.array([5.0, -3, 1, 2, -5])
print(f"성분합 : a {a.sum():.0f}, b {b.sum():.0f} -> 이미 중심화돼 있다")
beta = a @ b / (a @ a)
같다("17", beta, 2.5, "베타")
e = b - beta * a
같다("17", e, [0, -0.5, 1, -0.5, 0], "잔차")
같다("17", a @ e, 0.0, "잔차가 시장과 직교")
같다("17", b @ b, beta**2 * (a @ a) + e @ e, "피타고라스")
print(f" R^2 = {beta**2*(a@a)/(b@b):.4f}")
# 자료가 평균이 0 이 아니면
a2, b2 = a + 3, b + 8
print()
print(f"[17] 평균을 안 뺀 채로 같은 공식을 쓰면 : "
f"beta = {a2 @ b2 / (a2 @ a2):.4f} <- 2.5 가 아니다")
ac, bc = a2 - a2.mean(), b2 - b2.mean()
print(f" 중심화하고 다시 하면 : beta = {ac @ bc / (ac @ ac):.4f} <- 돌아온다")성분합 : a 0, b 0 -> 이미 중심화돼 있다
[17] 맞다 차이 0.00e+00 베타
[17] 맞다 차이 0.00e+00 잔차
[17] 맞다 차이 0.00e+00 잔차가 시장과 직교
[17] 맞다 차이 0.00e+00 피타고라스
R^2 = 0.9766
[17] 평균을 안 뺀 채로 같은 공식을 쓰면 : beta = 2.6364 <- 2.5 가 아니다
중심화하고 다시 하면 : beta = 2.5000 <- 돌아온다
print("[3] 쿠르노 : c_3 를 올려 가며 q_3 를 본다")
for c3 in [30, 35, 40, 43, 43.34, 45, 50]:
q = np.linalg.solve(M쿠르노, 100 - np.array([10.0, 20, c3]))
표 = "" if q[2] > 0 else " <- 음수 생산. 뜻이 없다"
print(f" c_3 = {c3:6.2f} -> q = {np.round(q, 3)}{표}")
print(f" 임계값은 40 이 아니라 130/3 = {130/3:.4f}")[3] 쿠르노 : c_3 를 올려 가며 q_3 를 본다
c_3 = 30.00 -> q = [30. 20. 10.]
c_3 = 35.00 -> q = [31.25 21.25 6.25]
c_3 = 40.00 -> q = [32.5 22.5 2.5]
c_3 = 43.00 -> q = [33.25 23.25 0.25]
c_3 = 43.34 -> q = [33.335 23.335 -0.005] <- 음수 생산. 뜻이 없다
c_3 = 45.00 -> q = [33.75 23.75 -1.25] <- 음수 생산. 뜻이 없다
c_3 = 50.00 -> q = [35. 25. -5.] <- 음수 생산. 뜻이 없다
임계값은 40 이 아니라 130/3 = 43.3333
print("[11] 열 합을 1 에 붙여 가며 (I-A)^-1 이 어떻게 되는지")
print(f"{'열 합':>8} {'rho(A)':>10} {'총승수':>12} {'전부 양수':>10}")
for s in [0.5, 0.7, 0.8, 0.9, 0.95, 0.99, 1.0, 1.05]:
As = A * (s / 0.8)
r = max(abs(np.linalg.eigvals(As)))
try:
L = np.linalg.inv(np.eye(3) - As)
승수 = L.sum() / 3
양수 = bool((L > 0).all())
print(f"{s:>8.2f} {r:>10.4f} {승수:>12.2f} {str(양수):>10}")
except np.linalg.LinAlgError:
print(f"{s:>8.2f} {r:>10.4f} {'특이':>12} {'-':>10}")
print(" 열 합이 1 을 넘으면 역행렬에 음수가 섞인다. 음수 산출은 뜻이 없다.")[11] 열 합을 1 에 붙여 가며 (I-A)^-1 이 어떻게 되는지
열 합 rho(A) 총승수 전부 양수
0.50 0.5000 2.00 True
0.70 0.7000 3.33 True
0.80 0.8000 5.00 True
0.90 0.9000 10.00 True
0.95 0.9500 20.00 True
0.99 0.9900 100.00 True
1.00 1.0000 특이 -
1.05 1.0500 -20.00 False
열 합이 1 을 넘으면 역행렬에 음수가 섞인다. 음수 산출은 뜻이 없다.
print("[10] 검사원이 서툴러질수록 되돌리기가 얼마나 위험해지는가")
print(" 이 꼴의 고윳값은 1, 1-e, 1-3e 라 e = 1/3 에서 특이해진다")
print(f"{'오분류율 e':>12} {'1-3e':>10} {'lambda_min':>12} "
f"{'증폭 배율':>12} {'조건수':>10}")
for eps in [0.05, 0.1, 0.2, 0.3, 0.32, 0.333]:
Me = np.array([[1-eps, eps, 0],
[eps, 1-2*eps, eps],
[0, eps, 1-eps]])
w = np.linalg.eigvalsh(Me)
print(f"{eps:>12.3f} {1-3*eps:>10.4f} {w.min():>12.4f} "
f"{1/w.min():>12.2f} {np.linalg.cond(Me):>10.2f}")
print(" e -> 1/3 에서 lambda_min -> 0 이고 증폭 배율이 폭발한다.")
print()
print(" e 가 1/3 을 넘으면 :")
for eps in [0.4, 0.49]:
Me = np.array([[1-eps, eps, 0],
[eps, 1-2*eps, eps],
[0, eps, 1-eps]])
w = np.linalg.eigvalsh(Me)
print(f" e = {eps:.2f} -> lambda_min = {w.min():+.4f} "
f"부호가 뒤집혔다. 조건수 {np.linalg.cond(Me):.2f} 가 되돌아오지만")
print(f" 그 사이에 특이점을 지났으므로 안전해진 것이 아니다.")[10] 검사원이 서툴러질수록 되돌리기가 얼마나 위험해지는가
이 꼴의 고윳값은 1, 1-e, 1-3e 라 e = 1/3 에서 특이해진다
오분류율 e 1-3e lambda_min 증폭 배율 조건수
0.050 0.8500 0.8500 1.18 1.18
0.100 0.7000 0.7000 1.43 1.43
0.200 0.4000 0.4000 2.50 2.50
0.300 0.1000 0.1000 10.00 10.00
0.320 0.0400 0.0400 25.00 25.00
0.333 0.0010 0.0010 1000.00 1000.00
e -> 1/3 에서 lambda_min -> 0 이고 증폭 배율이 폭발한다.
e 가 1/3 을 넘으면 :
e = 0.40 -> lambda_min = -0.2000 부호가 뒤집혔다. 조건수 5.00 가 되돌아오지만
그 사이에 특이점을 지났으므로 안전해진 것이 아니다.
e = 0.49 -> lambda_min = -0.4700 부호가 뒤집혔다. 조건수 2.13 가 되돌아오지만
그 사이에 특이점을 지났으므로 안전해진 것이 아니다.
print("[18] 두 열이 가까워질수록 계수가 얼마나 흔들리는가")
x1 = np.array([-3.0, -1, 1, 3])
y18 = np.array([-6.0, -4, 2, 8])
print(f"{'x2 = a*x1 + 나머지':>22} {'det(X^TX)':>12} {'조건수':>10} {'beta':>22}")
for a_ in [0.0, 0.3, 0.5, 0.7, 0.9, 0.99, 1.0]:
x2 = a_ * x1 / 3 * 2 + (1 - a_) * np.array([-1.0, -1, 1, 1])
X = np.c_[x1, x2]
G_ = X.T @ X
if abs(np.linalg.det(G_)) < 1e-12:
print(f"{a_:>22.2f} {np.linalg.det(G_):>12.2e} {'inf':>10} "
f"{'갈라낼 수 없다':>22}")
continue
be = np.linalg.solve(G_, X.T @ y18)
print(f"{a_:>22.2f} {np.linalg.det(G_):>12.4f} "
f"{np.linalg.cond(G_):>10.1f} {str(np.round(be,3)):>22}")[18] 두 열이 가까워질수록 계수가 얼마나 흔들리는가
x2 = a*x1 + 나머지 det(X^TX) 조건수 beta
0.00 16.0000 34.0 [2. 1.]
0.30 7.8400 77.7 [1.714 1.429]
0.50 4.0000 165.6 [1.333 2. ]
0.70 1.4400 502.6 [0.444 3.333]
0.90 0.1600 4968.2 [-4. 10.]
0.99 0.0016 519046.2 [-64. 100.]
1.00 -3.55e-14 inf 갈라낼 수 없다
두 열이 비례에 가까워질수록 가 0으로 가고 계수가 폭주한다. 0이냐 아니냐가 아니라 얼마나 0에 가까우냐가 문제라는 것이 숫자로 보인다.
직교 설계는 왜 좋은가¶
문제 21의 23 설계는 였다. 인자를 하나 빼도 나머지 계수가 안 바뀐다는 것을 확인해 보자.
A_ = np.array([-1.0, 1, -1, 1, -1, 1, -1, 1])
B_ = np.array([-1.0, -1, 1, 1, -1, -1, 1, 1])
C_ = np.array([-1.0] * 4 + [1.0] * 4)
X = np.c_[np.ones(8), A_, B_, C_]
y8 = np.array([45.0, 57, 39, 51, 49, 61, 43, 55])
같다("21", X.T @ X, 8 * np.eye(4), "X^T X = 8I")
be = np.linalg.solve(X.T @ X, X.T @ y8)
같다("21", be, [50, 6, -3, 2], "평균, 온도, 압력, 냉각")
같다("21", X @ be, y8, "잔차가 0")
print()
print("냉각 열을 빼고 다시 풀면 :")
be2 = np.linalg.solve(X[:, :3].T @ X[:, :3], X[:, :3].T @ y8)
print(f" {np.round(be2, 4)} <- 앞 세 계수가 하나도 안 바뀐다")
print()
print("문제 18 처럼 직교가 아니면 :")
Xn = np.c_[np.ones(4), x1, np.array([-1.0, -1, 1, 1])]
be3 = np.linalg.solve(Xn.T @ Xn, Xn.T @ y18)
be4 = np.linalg.solve(Xn[:, :2].T @ Xn[:, :2], Xn[:, :2].T @ y18)
print(f" 세 열 : {np.round(be3, 4)}")
print(f" 둘만 : {np.round(be4, 4)} <- 계수가 바뀐다")[21] 맞다 차이 0.00e+00 X^T X = 8I
[21] 맞다 차이 0.00e+00 평균, 온도, 압력, 냉각
[21] 맞다 차이 0.00e+00 잔차가 0
냉각 열을 빼고 다시 풀면 :
[50. 6. -3.] <- 앞 세 계수가 하나도 안 바뀐다
문제 18 처럼 직교가 아니면 :
세 열 : [0. 2. 1.]
둘만 : [0. 2.4] <- 계수가 바뀐다
Ai = np.array([[-1.0, 1, 0], [0, -1, 1], [1, 0, -1]])
print("각 행의 합 :", Ai.sum(axis=1))
for 이름, 마지막 in [("오전", 1/1000), ("오후", 1/900)]:
y = np.array([np.log10(100), np.log10(10), np.log10(마지막)])
루프 = y.sum()
있음 = abs(루프) > 1e-9
print(f"\n{이름} : y = {np.round(y, 4)}")
print(f" 루프 합 = {루프:.6f} -> "
f"{'차익거래가 있다' if 있음 else '무차익'}")
print(f" 한 바퀴 돌면 {10**루프:.4f} 배")
if not 있음:
p = np.array([0.0, 2, 3])
같다("22", Ai @ p, y, "전위 p = (0,2,3) 이 y 를 만든다")
print()
print("통화가 늘면 확인할 관계가 몇 개인가 (오일러 공식 m-n+1)")
for n in [3, 5, 10, 20]:
m = n * (n - 1) // 2
print(f" 통화 {n:2d}개, 간선 {m:3d}개 -> 독립 루프 {m-n+1:3d}개 "
f"(모든 삼각형은 {n*(n-1)*(n-2)//6}개)")각 행의 합 : [0. 0. 0.]
오전 : y = [ 2. 1. -3.]
루프 합 = 0.000000 -> 무차익
한 바퀴 돌면 1.0000 배
[22] 맞다 차이 0.00e+00 전위 p = (0,2,3) 이 y 를 만든다
오후 : y = [ 2. 1. -2.9542]
루프 합 = 0.045757 -> 차익거래가 있다
한 바퀴 돌면 1.1111 배
통화가 늘면 확인할 관계가 몇 개인가 (오일러 공식 m-n+1)
통화 3개, 간선 3개 -> 독립 루프 1개 (모든 삼각형은 1개)
통화 5개, 간선 10개 -> 독립 루프 6개 (모든 삼각형은 10개)
통화 10개, 간선 45개 -> 독립 루프 36개 (모든 삼각형은 120개)
통화 20개, 간선 190개 -> 독립 루프 171개 (모든 삼각형은 1140개)
def 총비용(Q):
return 1.8e8 / Q + 500 * Q
Qs = np.array([200, 400, 600, 800, 1200], dtype=float)
print(f"{'Q':>8} {'주문비':>12} {'보관비':>12} {'총비용':>12}")
for Q in Qs:
print(f"{Q:>8.0f} {50000*3600/Q:>12.0f} {1000*Q/2:>12.0f} {총비용(Q):>12.0f}")
Q별 = np.linspace(100, 1500, 2000)
최적 = Q별[np.argmin(총비용(Q별))]
print(f"\n수치로 찾은 최적 Q = {최적:.1f}, 손으로 푼 값 = "
f"{np.sqrt(1.8e8/500):.1f}")
print(f"최적에서 주문비 {50000*3600/600:.0f} = 보관비 {1000*600/2:.0f}")
print()
print("선형인지 확인 : C(Q) 가 Q 에 1차인가?")
for Q in [300.0, 600.0, 900.0]:
print(f" C({Q:.0f}) = {총비용(Q):.0f}, "
f"C({2*Q:.0f}) = {총비용(2*Q):.0f}, "
f"2*C({Q:.0f}) = {2*총비용(Q):.0f} "
f"{'같다' if abs(총비용(2*Q)-2*총비용(Q))<1e-6 else '다르다 -> 1차가 아니다'}") Q 주문비 보관비 총비용
200 900000 100000 1000000
400 450000 200000 650000
600 300000 300000 600000
800 225000 400000 625000
1200 150000 600000 750000
수치로 찾은 최적 Q = 600.1, 손으로 푼 값 = 600.0
최적에서 주문비 300000 = 보관비 300000
선형인지 확인 : C(Q) 가 Q 에 1차인가?
C(300) = 750000, C(600) = 600000, 2*C(300) = 1500000 다르다 -> 1차가 아니다
C(600) = 600000, C(1200) = 750000, 2*C(600) = 1200000 다르다 -> 1차가 아니다
C(900) = 650000, C(1800) = 1000000, 2*C(900) = 1300000 다르다 -> 1차가 아니다
이므로 1차가 아니다. 사다리 셋째 칸이 비고, 상자 0이다.
6. 자기채점¶
스물넷의 상자 번호를 적어 넣고 채점해 보자. 틀린 것은 상자 번호가 아니라 어느 문장을 잘못 읽었는지를 찍는다.
정답 = {
1: 1, 2: 1, 3: 1, 4: 1, 5: 1, 6: 1, 7: 1, 8: 1, 9: 2, 10: 1,
11: 4, 12: 1, 13: 1, 14: 4, 15: 4, 16: 4, 17: 3, 18: 3, 19: 3,
20: 3, 21: 3, 22: 2, 23: 1, 24: 0,
}
힌트 = {
0: "미지수가 분모와 분자에 동시에 있는지 보라",
1: "줄 수와 미지수 수가 같고 독립인지 보라",
2: "줄이 모자라거나 겹쳐서 답이 하나로 안 정해지는 자리를 보라",
3: "줄이 미지수보다 많아 다 맞힐 수 없는 자리를 보라",
4: "같은 조작이 되풀이되는지, 장기 거동을 묻는지 보라",
5: "부피·크기·방향을 묻는지 보라",
}
# 여기에 자기 답을 적어 넣는다. 비워 두면 정답으로 채점한다.
내답 = dict(정답)
맞은수 = 0
for 번호 in sorted(정답):
참 = 정답[번호]
낸 = 내답.get(번호)
if 낸 == 참:
맞은수 += 1
else:
print(f" 문제 {번호:2d} : 상자 {낸} 라고 했는데 {참} 이다. {힌트[참]}")
print(f"\n{맞은수}/{len(정답)} 맞음")
from collections import Counter
print("상자별 분포 :", dict(sorted(Counter(정답.values()).items())))
24/24 맞음
상자별 분포 : {0: 1, 1: 12, 2: 2, 3: 5, 4: 4}
정리¶
| 확인한 것 | 결과 |
|---|---|
| 레온티예프가 총산출을 복원하는가 | 으로 정확히 |
| (문제 2 와 3 이 물려 있다) | 성립 |
| 같은 표를 와 로 읽기 | 값과 복제 방법이 둘 다 나온다 |
| 멱영이라 유한 항에서 끝난다 | |
| 차분의 역이 누적합 | 이 하삼각 전부 1 |
| 이라 수렴, 20라운드에 99% | |
| 기본행렬 도 같은 급수 | 성립 |
| 두 대면 두 배로 버티는가 | 아니다. 1.5배다 |
| 30분 뒤 가용도가 장기 평균인가 | 아니다. 91.2% 대 80% |
| 전치를 빼먹으면 | 정상 점유율이 균등분포로 나온다 |
| 중심화를 빼먹으면 | 베타가 2.5 가 아닌 값으로 나온다 |
| 열 합이 1 을 넘으면 | 역행렬에 음수가 섞인다 |
| 두 열이 비례에 가까워지면 | 계수가 폭주한다 |
| 직교 설계면 | 인자를 빼도 나머지 계수가 안 바뀐다 |
| EOQ 가 인가 | 아니다. |
가장 값진 줄은 "아니다"가 적힌 줄들이다. 세운 것이 맞았는지 재는 일보다 틀린 번역이 얼마나 그럴듯한 답을 내놓는지 보는 일이 더 무섭고 더 유익하다.
전치 한 번을 빼먹으면 점유율이 균등분포로 나온다. 그럴듯하다. 그래서 위험하다.