이 포스팅은 Jupiter Notebook로 작성되었습니다.
그러니 직접 Python 코드를 실행시키고 싶으신 분은 아래 Link에서 다운로드 받으셔서 실행할 수 있습니다.
https://github.com/DevStarSJ/Study/tree/master/Blog/Python/DoingMathWithPython
파이썬으로 풀어보는 수학
- 원서명 : Doing Math with Python: Use Programming to Explore Algebra, Statistics, Calculus, and More! (ISBN 9781593276409)
- 지은이 : 아미트 사하(Amit Saha)
- 원서 및 관련자료 : https://www.nostarch.com/doingmathwithpython
- 번역서 : http://www.acornpub.co.kr/book/doing-math-with-python
5장 집합과 확률
1. 집합은 무엇인가 ?
집합(set)은 개별 객체의 모음(collection)입니다. 통상 객체를 원소(element)나 구성원(member)라고 부릅니다. 집합과 모음의 다른점으로는 집합은 동일한 2개의 원소를 가질 수 없습니다.
SymPy를 사용해 파이썬에서 집합으로 작업하는 방법에 대해서 살펴보겠습니다.2. 집합 생성
수학기호로 집합을 표현할 때는
{ 2, 4, 6 }와 같이 중괄호를 사용합니다. 파이썬에서는 FiniteSet 클래스를 사용하여 표현할 수 있습니다.
In [3]:
from sympy import FiniteSet
s = FiniteSet(2, 4, 6)
s
Out[3]:
In [4]:
type(s)
Out[4]:
동일한 집합에 정수, 분수, 부동소수점수를 같이 저장할 수도 있습니다.
In [6]:
from fractions import Fraction
FiniteSet(1, 1.5, Fraction(1,5))
Out[6]:
카디널리티(cardinality)란 집합의 구성원 수를 의미합니다. len()함수를 이용해 계산이 가능합니다.
In [7]:
len(s)
Out[7]:
대상 집합에 숫자가 존재하는지 여부는 in연산자를 사용해 알 수 있습니다.
In [8]:
3 in s
Out[8]:
In [9]:
4 in s
Out[9]:
공집합(empty set)을 생성하기 위해서는 인자가 없이 생성하면 됩니다.
In [10]:
FiniteSet()
Out[10]:
리스트나 튜플을 인자로 전달해서 집합을 생성할 수도 있습니다.
In [12]:
members = [1, 2, 3]
FiniteSet(*members)
Out[12]:
앞서 언급했듯이 집합 내에는 중복된 값을 허용하지 않습니다. 같은 값을 여러번 넣어도 한 번만 추가되고 나머지는 모두 무시됩니다.
In [16]:
s = FiniteSet(1, 3, 2, 3)
s
Out[16]:
그리고, 입력 순서와 저장되는 순서는 무관합니다. 집합 내의 구성원에 대한 순서는 별도로 저장하고 있지 않기 때문입니다.
In [18]:
for m in s:
print(m)
두 집합의 구성원이 서로 다른 순서로 저장되더라도, 모든 요소들이 같다면 두 집합은 같은 집합으로 취급합니다.
In [19]:
s = FiniteSet(3, 5, 7)
t = FiniteSet(5, 7, 3)
s == t
Out[19]:
3. 부분집합, 초집합, 파워집합
집합 s의 모든 구성원이 집합 t의 구성원일 경우 s는 t의
부분집합(subset)이라고 정의합니다. 파이썬에서는 is_subset() 함수를 사용해서 확인이 가능합니다.
In [20]:
s = FiniteSet(1)
t = FiniteSet(1, 2)
s.is_subset(t)
Out[20]:
In [21]:
t.is_subset(s)
Out[21]:
In [22]:
s.is_subset(s)
Out[22]:
공집합은 모든 집합의 부분집합이며, 모든 집합은 자기 자신이 부분집합입니다.
초집합(superset)은 부분집합의 반대 개념으로 집합t가 집합s의 모든 구성원을 포함할 경우 집합t는 집합s의 초집합이라고 부릅니다. 파이썬에서는 is_superset() 함수를 사용해서 확인이 가능합니다.
In [23]:
s.is_superset(t)
Out[23]:
In [24]:
t.is_superset(s)
Out[24]:
In [25]:
s.is_superset(s)
Out[25]:
파워집합(powerset)은 모든 가능한 부분집합입니다. 모든 집합은 2 ** cadinality 만큼의 부분집합을 가집니다. (공집합, 자기자신을 포함) 파이썬에서는 powerset() 함수를 사용해 찾아낼 수 있습니다.
In [26]:
s = FiniteSet(1, 2, 3)
ps = s.powerset()
ps
Out[26]:
In [27]:
len(ps)
Out[27]:
동일하지 않은 집합t와 집합s가 있을 경우, 집합s가 집합t의 부분 집합일 경우 집합t는 집합s의 초집합이라고 할 수 있습니다.
is_proper_subset() , is_proper_superset() 함수를 사용해서 부분집합, 초집합 관계를 확인 할 수 있습니다.
In [28]:
t = FiniteSet(1, 2, 3, 5)
s.is_proper_subset(t)
Out[28]:
In [29]:
t.is_proper_superset(s)
Out[29]:
In [30]:
s.is_proper_subset(s)
Out[30]:
In [31]:
s.is_proper_superset(s)
Out[31]:
4. 집합 연산
합집합
두 집합의 모든 구성원을 포함하는 집합입니다. (중복된 구성원은 하나만 포함됩니다.)
In [1]:
from sympy import FiniteSet
s = FiniteSet(1, 2, 3)
t = FiniteSet(2, 4, 6)
s.union(t)
Out[1]:
교집합
두 집합에 공통적으로 존재하는 구성원들로만 이루어진 집합입니다.
In [2]:
s.intersect(t)
Out[2]:
카르테지안 곱
두 집합의 구성원들을 택해 모든 가능한 쌍으로 구성된 집합입니다.
카르테지안 곱의 카디널리티는 개별 집합의 카디널리티의 곱입니다.
In [4]:
p = s*t
p
Out[4]:
In [5]:
for e in p:
print(e)
In [6]:
len(p) == len(s) * len(t)
Out[6]:
지수연산자를 이용해서 설정한 횟수만큼의 곱도 가능합니다.
In [7]:
u = FiniteSet(1, 2)
p = u**3
p
Out[7]:
In [8]:
for e in p:
print(e)
다중변수 집합에 공식 적용
다음 수식은 추의 길이별 주기값을 구하는 수식입니다.
- T : 추가 한번 왕복하는데 소요되는 시간
- L : 추의 길이
- 상수값 : pi, g(지역의 중력가속도 : 대략 9.8 m/s^2)
길이에 따라 추의 주기가 어떻게 변하는지 알고 싶다면 위 수식에
L값을 변경해가면서 입력해보면 됩니다.
In [10]:
from sympy import FiniteSet, pi
def time_period(length):
g = 9.8
T = 2 * pi * (length / g)**0.5
return T
L = FiniteSet(15, 18, 21, 22.5, 25)
for e in L:
t = time_period(e/100)
print('Length: {0} cm Time Period: {1:.3f} s'.format(float(e), float(t)))
입력으로 사용된 집합의 값이 cm 단위이기 때문에 함수에 전달시 100으로 나누어서 전달하였습니다.
서로 다른 3곳(중력이 다른 곳)에서 실험을 한다고 가정해 보겠습니다. (적도 : 9.78 , 북극 9.83, 호주 9.8)
In [14]:
def time_period(length, g):
return 2 * pi * (length / g)**0.5
L = FiniteSet(15, 18, 21, 22.5, 25)
G = FiniteSet(9.8, 9.78, 9.83)
print('{0:^15}{1:^15}{2:^15}'.format('Length(cm)','Gravity(m/m^2)','Time Period(s)'))
for e in L * G:
t = time_period(e[0]/100, e[1])
print('{0:^15}{1:^15}{2:^15.3f}'.format(float(e[0]),float(e[1]), float(t)))
5. 확률
- 실험(experiment) : 각각 가능한 확률에 대한 테스트, 실험을 한 번 실행하는 것을 시도(trial)이라고 함. 예를 들어 주사위 던지기, 카드 뽑기
- 표본공간(
S) : 모든 가능한 실험 결과들의 집합. 예를 들어 6면 주사위를 한 번 던진 경우 표본공간S = {1, 2, 3, 4, 5, 6} - 사건(
E) : 표본공간의 부분집합. 예를 들어 6면 주사위의 표본공간 중 숫자 3이 나올 확률
특정 사건이 일어날 확률(
P(E))는 해당 사건의 개수(n(E))를 전체 표본공간의 개수(n(S))로 나눈 값입니다. 주사위를 던져서 3이 나올 확률은 다음과 같습니다.S = { 1, 2, 3, 4, 5, 6 } E = { 3 } n(S) = 6 n(E) = 1 P(E) = 1 / 6
이를 함수로 작성해 보겠습니다.
In [23]:
def probability(space, event):
return len(event) / len(space)
위 함수를 이용해서 1 ~ 20 사이의 숫자 중 소수(Prime number)일 확률을 구해보도록 하겠습니다.
In [31]:
def check_prime(number):
if number != 1:
for factor in range(2,number):
if number % factor == 0:
return False
else:
return False
return True
def get_primes(number):
primes = [x for x in range(2, number + 1) if check_prime(x) ]
return FiniteSet(*primes)
space = FiniteSet(*range(1,21))
event = get_primes(20)
probability(space, event)
Out[31]:
In [32]:
event
Out[32]:
사건A나 사건B의 확률
A 사건과 B 사건의 집합을 합집합(
union)한 다음에 확률을 계산하면 됩니다.- S : 6면주사위를 1번 던짐 -> { 1, 2, 3, 4, 5, 6 }
- A : 소수 -> { 2, 3, 5 }
- B : 홀수 -> { 1, 3, 5 }
일 경우의 사건 A나 B일 확률은 다음과 같습니다.
In [19]:
S = FiniteSet(*range(1,7))
A = FiniteSet(2, 3, 5)
B = FiniteSet(1, 3, 5)
len(A.union(B)) / len(S)
Out[19]:
사건 A 이면서 사건 B일 확률
두 집합의 교집합(
intersect)의 확률을 계산하면 됩니다.
In [20]:
len(A.intersect(B)) / len(S)
Out[20]:
6. 랜덤 숫자 생성
파이썬에서 랜덤숫자를 생성하려면 먼저 표준 라이브러리
random을 포함시켜야 합니다. 주로 사용되는 랜덤함수는 다음의 2가지 정도만 알아도 됩니다.- randint(from, to) : from 에서 to 사이의 숫자를 리턴합니다. int값을 넣어줘야 합니다.
- random() : 0에서 1사이의 부동소수점 숫자를 생성합니다.
주사위를 굴려서 총합이 20될 때까지 몇 번을 던져야 하는지를 랜덤을 통해서 구현해 보겠습니다.
In [24]:
import random
target_score = 20
def roll():
return random.randint(1,6)
def play_game():
score = 0
num_rolls = 0
while score < target_score:
dice = roll()
num_rolls += 1
print('Rolled: {0}'.format(dice))
score += dice
print('Score of {0} reached in {1} rolls'.format(score, num_rolls))
play_game()
In [25]:
play_game()
목표점수 달성이 가능한가 ?
이번에는 목표로 한 점수가 최대던지기 횟수 내에 달성이 가능한지 그 여부 및 확률을 계산해주는 프로그램을 작성해 보겠습니다.
In [1]:
from sympy import FiniteSet
from random import randint
def find_prob(target_score, max_rolls):
dice = FiniteSet(*range(1,7))
space = dice**max_rolls
eventNum = 0
for e in space:
n = sum(e)
if (n >= target_score):
eventNum += 1
return eventNum / len(space)
find_prob(10,2)
Out[1]:
In [7]:
find_prob(20, 3)
Out[7]:
주사위를 3번 던져서 최고값은 18이므로 20이 나올 수 있는 확률은 0% 입니다.
위에서 사용한 랜덤은 균일 랜덤 숫자(uniform random number)였습니다. 비균일 랜덤을 생성하려면 어떻게 해야 할까요 ? 가장 간단하게 생각해 볼 수 있는 방법은 높은 확률에 더 넓은 영역을 지정하여 랜덤을 수행하는 것입니다.
동전을 던졌을 때 앞면(True)이 나올 확률이 2/3이고, 뒷면(False)가 나올 확률이 1/3일 경우
In [9]:
import random
def nonUniformToss():
if random.random() < 2/3:
return True
return False
trueNum = 0
for i in range(10000):
if nonUniformToss() == True:
trueNum += 1
trueNum
Out[9]:
위 함수를 조금 더 일반적으로 구현하여 재사용 가능하도록 해보겠습니다. 인자로 확률값을 넣은 리스트를 전달받아서 호출시 마다 해당 index를 리턴받도록 하겠습니다
In [19]:
import random
def get_index(probabilities):
max = sum(probabilities)
acc = 0
rand = random.random() * max
for idx, percent in enumerate(probabilities):
acc += percent
if rand < acc:
return idx
return len(probabilities)
이 함수를 이용해서 다음의 경우에 대해서 시뮬레이션 해보겠습니다.
다음 확률로 지폐를 배분하는 ATM기의 경우 만번의 지폐를 배분했을 때 각각 몇회 배분되었는지를 구해보겠습니다.
- $5 : 1/6
- $10 : 1/6
- $20 : 1/3
- $50 : 1/3
In [20]:
from collections import defaultdict
d = defaultdict(int)
keys = ['$5', '$10', '$20', '$50']
probabilities = [ 1, 1, 2, 2 ]
for i in range(10000):
d[keys[get_index(probabilities)]] += 1
d
Out[20]:
프로그래밍 연습
1. 벤다이어그램을 사용하여 집합 간의 관계를 가시화
matplotlib_venn 패키지를 이용하여 벤다이어그램을 그릴 수 있습니다.
예를 들어서 20이하의 소수와 홀수를 그리는 벤다이어그램을 구현한 코드를 살펴보겠습니다.
In [35]:
def get_odds(number):
odds = [x for x in range(1, number + 1) if x % 2 == 1]
return FiniteSet(*odds)
odds = get_odds(20)
odds
Out[35]:
In [36]:
primes = get_primes(20)
primes
Out[36]:
In [40]:
from matplotlib_venn import venn2
import matplotlib.pyplot as plt
def draw_venn(sets, lables):
venn2(subsets=sets, set_labels=lables)
plt.show()
draw_venn([odds, primes], ['odds', 'primes'])
(학번, 축구를 좋아하는지 여부, 다른 운동을 좋아하는지 여부)의 3개의 컬럼을 가진 csv 파일을 읽어서 그 결과를 보여주는 벤다이어그램을 출력하세요. 좋아한다는 것은 1, 안 좋아한다는 것은 0으로 표시하도록 하겠습니다. (csv를 읽어서 list를 리턴하는 함수는 Ch.03에서 작성해 놓은 함수를 재사용 하겠습니다.)
In [43]:
import csv
def CsvToList(filename, colHeaderLen, rowHeaderLen):
reader = csv.reader(open(filename))
columns = len(next(reader)) # pass column header
for i in range(1, colHeaderLen):
next(reader)
data = []
for i in range(rowHeaderLen, columns):
data.append([])
for row in reader:
for i in range(rowHeaderLen, columns):
data[i - rowHeaderLen].append(float(row[i]))
return data
In [69]:
lists = CsvToList('files/ch05.venn.data.csv', 1, 0)
lists[2]
footballs = FiniteSet(*[ int(id) for id,football in zip(lists[0], lists[1]) if football == 1 ])
others = FiniteSet(*[ int(id) for id,other in zip(lists[0], lists[2]) if other == 1 ])
draw_venn([footballs, others],['Football','Other'])
In [70]:
sum(lists[1])
Out[70]:
In [72]:
sum(lists[2])
Out[72]:
2. 대수의 법칙 (기대값 계산)
기대값은 모든 경우의 수의 값에 그 확률을 곱한 값을 의미합니다.
그럼 주사위를 던졌을 경우의 기대값은 다음과 같습니다.
In [73]:
e = 1*(1/6) + 2*(1/6) + 3*(1/6) + 4*(1/6) + 5*(1/6) + 6*(1/6)
e
Out[73]:
주사위를 던지는 회수를 늘려가면서 기대값대로 나오는지를 입증하는 프로그램을 작성하세요.
먼저, 재사용 가능한 형태로 기대값을 구하는 함수부터 구하는 만들어 보겠습니다.
In [75]:
def GetExpectedValue(events, probabilities):
elements = [x*p for x,p in zip(events, probabilities)]
return sum(elements)
위 함수가 제대로 동작하는지 앞에서 직접 계산해본 주사위를 던졌을 경우의 기대값을 함수를 통해서 구해보겠습니다.
In [76]:
events = [1,2,3,4,5,6]
probabilities = [1/6,1/6,1/6,1/6,1/6,1/6]
GetExpectedValue(events, probabilities)
Out[76]:
그럼 이제 주사위를 던지를 횟수를 늘려가면서 실제로 비슷하게 나오는지 살펴보겠습니다.
In [79]:
def GetExpectedValueInMultipleDrawing(num):
acc = 0
for i in range(num):
acc += events[get_index(probabilities)]
return acc / num
In [84]:
GetExpectedValueInMultipleDrawing(100)
Out[84]:
In [83]:
GetExpectedValueInMultipleDrawing(10000)
Out[83]:
In [82]:
GetExpectedValueInMultipleDrawing(100000)
Out[82]:
3. 돈이 떨어지기 전에 토스 시도 횟수는 ?
동전을 던져서 앞면이 나오면 1를얻고,뒷면이나오면 1.5를 잃는 게임이 있습니다. (누가봐도 불공평한 이딴 겜을 누가 할지는 모르겠습니다만... 어쨌든 하는 사람이 있다고 합시다.) 잔고가 0이 되면 게임이 끝납니다. 이 게임을 시뮬레이션하는 코드를 작성하세요.
In [94]:
events = [1, -1.5]
probabilities = [0.5, 0.5]
def PlayDrawGame(balances):
cnt = 0
while balances > 0:
index = get_index(probabilities)
balances += events[index]
cnt += 1
print('{0:3} {1} ! Current amount: {2:.1f}'.format(cnt, 'Heads' if index == 0 else 'Tails', balances))
print('Game Over ! Coin tosses {0:3} times'.format(cnt))
In [100]:
PlayDrawGame(10)
4. 카드뭉치 섞기
52장의 트럼프카드를 섞는(shuffling)하는 프로그램을 작성하세요.
먼저 파이썬 표준라이브러리의
random 모듈에 있는 shuffle()함수에 대해 살펴보겠습니다.
In [109]:
import random
x = [1, 2, 3, 4]
random.shuffle(x)
x
Out[109]:
트럼프에는 총 52장의 카드가 있으니 숫자 1 ~ 52를 이용해서 셔플을 한 후에 해당 숫자를 트럼프 카드와 매핑을 해주는 클래스를 이용하는 방식으로 구현해 보겠습니다. 먼저 1 ~ 52까지의 숫자를 트럼프 카드로 매핑하는 클래스를 만들어 보겠습니다.
In [114]:
class Trump:
suits = ['spades','diamonds','clubs','hearts']
ranks = ['ace','two','three','four','five','six','seven','eight','nine','ten','jack','queen','king']
def __init__(self, idx):
self.suit = self.suits[idx // 13]
self.rank = self.ranks[idx % 13]
def print(self):
print('{0} of {1}'.format(self.rank, self.suit))
이제 위 클래스를 이용해서 카드를 셔플해 보겠습니다.
먼저 0 ~ 52 사이의 숫자를 셔플하겠습니다.
In [112]:
import random
cards = list(range(52))
random.shuffle(cards)
print(cards)
이 결과를
Trump 클래스를 이용해서 출력해 보겠습니다.
In [118]:
trumps = [Trump(card) for card in cards]
for trump in trumps:
trump.print()
5. 원의 면적 추정
원의 면적은
pi * r^2이라는 공식이라는 것을 알고 있습니다만, 이것이 정말로 맞는지를 다른 방법으로 추정해 보고자 합니다. 한쪽면의 길이가 2r인 정사각형 안에 반지름 r인 원으로 구성된 다트보드가 있다고 생각해 봅시다. 여기에 다트를 던졌을 경우 원 안에 명중한 갯수를 N개라고 하고 나머지는 원 안이 아닌 정사각형 보드위에 맞췄다고 가정하겠습니다. (M) 정사각형 보드 밖으로 나가는 경우는 생각하지 않겠습니다. 이럴 경우 원안에 명중할 확률은 f = N / (N + M)이 됩니다. 정사각형의 면적이 A일 경우 A * f는 원의 면적과 같을 것입니다.
이것을 직접 프로그램으로 구현하여 던지는 횟수를 다르게 하여 실제 원의 면적과 얼마나 차이가 나는지를 살펴보세요.
먼저 실제 원의 면적을 구해보겠습니다. (반지름이 2이라 가정하겠습니다.)
In [122]:
from sympy import pi, N
r = 2
N(pi *r**2)
Out[122]:
이제 다트를 던졌을 때 원안에 맞은 확률을 구해야 합니다. 이 확률을 구하기 위해서 다트를 던지는 것의 시뮬레이션 결과 중 명중한 것을
1, 그렇지 않은 것을 0으로 하겠습니다. 그 결과에 정사각형의 면접과 곱한한 결과가 위에 계산결과와 맞는지를 비교해 보겠습니다.
원 안에 맞았다는 것을 어떻게 구현해야 할까요 ? 다트가 명중한 곳과 원의 중심간의 거리가 반지름(2)보다 작은 경우 명중했다고 생각하겠습니다. 계산상 편의를 위해 원의 중심을
(0, 0)으로 생각하고 좌표를 랜덤으로 생성하겠습니다. (-2 ~ 2)
In [131]:
import random
def Dart():
x = random.uniform(-2, 2)
y = random.uniform(-2, 2)
dist = (x*x + y*y)**0.5
return True if dist < 2 else False
def ProbabilityOfDart(number):
shootings = [1 if Dart() else 0 for x in range(number)]
return sum(shootings) / number
def GetExpectedOfDart(number, length):
return ProbabilityOfDart(number) * length * length
In [135]:
GetExpectedOfDart(100, 4)
Out[135]:
In [136]:
GetExpectedOfDart(1000, 4)
Out[136]:
In [137]:
GetExpectedOfDart(10000, 4)
Out[137]:
In [138]:
GetExpectedOfDart(100000, 4)
Out[138]:
댓글 없음:
댓글 쓰기