[data analyis] 확률과 주요통계량: 모멘트와 기대값 예제

예 1)

다음의 확률밀도함수(pdf)를 갖는 확률연속변수의 기대값은?

$$f(x)=\begin{cases}c(x^3+x^2+1)& 0\lt x \lt 10\\ 0&\text{otherwise} \end{cases}$$

위 함수에서 범위 (0, 10)에서의 적분값은 1이 되어야 합니다. 이 조건을 적용하여 상수 c가 계산될 수 있습니다.

적분 계산은 파이썬 라이브러리 sympy.integrate() 함수를 적용합니다. 또한 그 적분의 결과에서 표현되는 미지수 c는 sympy.solve() 함수를 사용하여 결정할 수 있습니다.

c, x=symbols("c, x")
f=c*(x**3+x**2+1)
F=f.integrate((x, 0, 10))
print(F)

8530*c/3

C=solve(Eq(F, 1), c);C

[3/8530]

위 결과를 c로 치환한 새로운 함수를 사용하여 기대값을 계산합니다.

f=f.subs(c, C[0])
print(f)

3*x**3/8530 + 3*x**2/8530 + 3/8530

위 결과인 확률함수에 대한 기대값의 계산은 식 1과 같습니다.

$$\tag{식 1}E(x) = \int^{10}_0 x\left(\frac{3}{8530}x^3+\frac{3}{8530}x+\frac{3}{8530} \right)$$

E=integrate(x*f, (x, 0, 10))
print(E)

6765/853

round(float(E), 3)

7.931

확률변수가 변수 x에 대한 함수 g(x)를 정의 할 수 있는 경우 기대값은 식 2와 같이 그 함수에 대한 확률 함수 f(x)와의 곱으로 계산할 수 있습니다.

$$\tag{식 2}E(g(x))=\begin{cases}\sum_{x\in S}g(x)f(x)& x:\text{이산변수},\; S: \text{표본공간}\\ \int^\infty_{-\infty}g(x)f(x)\, dx& x:\text{연속변수}\end{cases}$$

예 2)

주사위를 4번 던지는 경우 홀수가 나오는 횟수를 X하면 E(x) ?

이번 실험에서는 주사위가 홀수(1, 3, 5)이면 1로 설정하고, 주사위 숫자가 짝수(2, 4, 6)이면 0으로 설정합니다. 확률변수 x는 홀수가 나오는 횟수이므로 표본공간 S는 다음과 같습니다.

S = {0, 1, 2, 3, 4}

주사위를 한 번 던져 홀수와 짝수가 나오는 확률은 각각 0.5이며 각각의 시행은 독립입니다. 이 시행의 표본공간은 itertools.combinations_with_replacement(확률변수, 시행횟수) 함수로 나타낼 수 있습니다.

• itertools.combiations(x, r)
  • iterable 객체 x(확률변수)에 대해 지정한 수(r)만큼 선택한 조합의 모든 경우를 반환
  • iterable는 리스트, 튜플 형과 같이 각 요소의 순서에 맞게 호출할 수 있는 객체입니다.
  • x의 길이 = n → 0 ≤ r ≤ n
  • r > n 일 경우 itertools.combinations_with_replacement(x, r)을 사용

S=list(itertools.combinations_with_replacement([0,1], 4));S

[(0, 0, 0, 0), (0, 0, 0, 1), (0, 0, 1, 1), (0, 1, 1, 1), (1, 1, 1, 1)]

위 결과 각각의 확률함수(pmf)는 조합과 각각의 확률을 곱하는 것으로 결정할 수 있습니다. 예를 들어 홀수의 선택이 0개일 경우의 경우의 수는 5개 중에 0개를 선택하는 조합의 수로 ₅C₀이 됩니다. 이 경우의 수에서 홀수가 선택될 확률 $\left(\frac{1}{2}\right)^0$과 짝수가 선택될 확률 $\left(\frac{1}{2}\right)^5$이 곱으로 계산됩니다. 확률함수 f(x)는 식 3과 같이 나타낼 수 있습니다. 기대값은 각 확률변수에 대응하는 pmf와 곱들의 합으로 결정됩니다.

$$\begin{align}\tag{식 3} f(x)&={}_n\mathrm{C}_s \left(\frac{1}{2}\right)^s \left(\frac{1}{2}\right)^{n-s} \\ E(x)&=\sum xf(x)\end{align}$$

이 시행은 조합 계산을 위한 함수인 scipy.special.comb(n, k)를 사용하여 계산할 수 있습니다. 예를 들어 4번시행에서 홀수의 출현이 0인 경우는 다음과 같이 계산됩니다.

scipy.special.comb(4, 0)*(1/2)**0*(1/2)**4

0.0625

기대값은 위 함수를 사용하여 표본공간 S의 각 사건에 대한 확률함수를 계산하는 것으로 결정합니다.

S=np.array([0,1,2,3, 4])
p=np.array([scipy.special.comb(4, i)*(1/2)**i*(1/2)**(4-i) for i in S])
print(p)

[0.0625 0.25   0.375  0.25   0.0625]

E=np.sum(S*p); E

2.0

각 사건의 확률함수인 p를 시각화하면 그림 1와 같습니다.

그림 1. 4번의 주사위 시행에 대한 확률분포.

plt.figure(figsize=(4,3))
plt.bar([0, 1, 2, 3, 4], p)
plt.xlabel("# of odd")
plt.ylabel("prob")
plt.show()

이 예제는 성공과 실패에 2가지 확률을 다룹니다. 이러한 시행에 대한 확률의 분포는 이항분포를 구성하며 이 분포의 확률함수의 전형적인 형태는 그림 1과 같습니다. 이항분포에서 각 확률변수에 대한 확률은 scipy.stats.binom.pmf(성공횟수, 총시행횟수, 시행당 확률) 함수를 사용하여 계산할 수 있으며 이 결과를 기반으로 기대값을 결정할 수 있습니다.

scipy.stats 모듈은 다양한 분포에 대해 확률함수, 기대값 등의 여러 통계량을 반환하는 함수들을 제공합니다.

p1=scipy.stats.binom.pmf(S, 4, 1/2);
print(p1)

[0.0625 0.25   0.375  0.25   0.0625]

E1=np.sum(S*p1); E1

2.0

또한 다음 코드와 같이 기대값을 직접적으로 확인할 수 있습니다.

scipy.stats.binom.stats(4, 1/2, moments='m')

2.0

예 3)

예제 2에서 확률변수를 X 대신 X²을 사용한다면 E(X²)?

S=np.array([0,1,2,3, 4])
E2=np.sum(SN**2*p1); E2

5.0

예 4)

다음은 연속확률변수의 PDF를 정의한 것입니다.

$$f(x)=\begin{cases}1& 0\lt x \lt 1\\ 0 & \text{otherwise} \end{cases}$$

연속확률변수 X를 기반으로 새로운 확률변수 Y = g(x) = e^x의 기대값과 e^x3의 기대값을 결정합니다.

기대값은 확률변수와 확률함수의 곱의 합으로 식 4와 같이 계산할 수 있습니다.

$$\begin{align}\tag{식 4} E(e^x)&=1\int^1_0 e^xf(x)\, dx\\&=\int^1_0 e^x\, dx\\&=e-1\\ E(e^{x^3})&=1\int^1_0 e^{x^3}f(x)\, dx\\&=\int^1_0 e^{x^3}\, dx\end{align}$$

식 4에서 두번째 함수 e^x3의 적분은 이 함수의 무한급수를 적용하여 실행할 수 있지만 간단하지 않습니다. 대신에 sympy의 integrate() 함수를 적용할 수 있으며 그 결과는 기호나 수로 구성된 식이 될 수 있습니다. 이러한 결과는 N()를 사용하여 수로 변환할 수 있습니다.

sympy 패키지는 수학의 다양한 식들의 계산을 위한 함수들을 제공합니다.

from sympy import *
x=symbols('x')
g1=exp(x)
E1=integrate(g1, (x, 0, 1));N(E1, 3)

1.72

g2=exp(x**3)
G=integrate(g2,(x, 0, 1))
simplify(G)

$- \frac{\left(-1\right)^{\frac{2}{3}} \gamma\left(\frac{1}{3}, e^{i \pi}\right)}{3} $

N(G, 3)

1.34

예 5)

통계학 강의에 두 개의 책이 사용된다고 합니다. 두 책의 구입여부는 독립적이라고 가정합니다. 다시말하면, 주교재의 구입여부가 부교재의 구입여부에 영향을 주지 않는다는 가정합니다. 그 가정하에 학생 당 책 구입 또는 구입하지 않을 확률은 동일하므로 확률변수로 고려할 수 있습니다. 확률 변수는 다음과 같이 구성됩니다.

• 두 권 모두 구입
• 주교재만 구입
• 부교재만 구입
• 두 권 모두 구입하지 않음

표 1은 과거의 학생들의 책 구입경향을 나타낸 것입니다.

표 1 교재 구입 현황

case	probiity
no both books	10%
main book	45%
sub-book	25%
both books	20%

주교재와 부교재의 가격이 각각 10000원과 7000원인 경우 강좌의 수강생 한명 당 평균적으로 지출할 책 구입비를 계산합니다.

책 구입을 위한 기대값을 계산하는 것으로 위에서 제시한 기존의 자료들을 사용하여 확률변수에 대해 다음과 같이 요약할 수 있습니다.

x=np.array([0,10000,7000,17000])
p=np.array([0.1,0.45,0.25,0.2])
da=pd.DataFrame([x, p], index=["가격","확률"], columns=[1,2,3,4]); da

	1	2	3	4
가격	0	10000	7000	17000
확률	0.1	0.45	0.25	0.2

da.product(axis=0).sum()

9650.0

예 6)

확률변수 X의 조건은 다음과 같습니다.

Rx =	{−3,	−2, −1, 0, 1, 2, 3}
f(x) =	1
	7
Rx: Range of random variable X

변수 x를 기반으로 하는 새로운 랜덤변수 y=2|x+1|의 범위와 기대값을 결정합니다.

함수 y에 의해 새로운 변수 범위를 결정합니다.

Rx=np.array([-3,-2,-1, 0, 1, 2, 3])
Ry=2*abs(Rx+1)
print(Ry)

[4 2 0 2 4 6 8]

변수 y의 범위 R_y로 부터 빈도와 확률을 계산합니다.

val, fre=np.unique(Ry, return_counts=True)
print(val)
print(fre)

[0 2 4 6 8]
[1 2 2 1 1]

prob=[Rational(i, 7) for i in fre]
print(prob)

[1/7, 2/7, 2/7, 1/7, 1/7]

data=pd.DataFrame([val, prob], index=['Y', "prob"]); data

	0	1	2	3	4
Y	0	2	4	6	8
prob	1/7	2/7	2/7	1/7	1/7

E=data.product(axis=0).sum()
print(E)

26/7