정규분포 관련 성질

정규분포 관련 성질 -> 더 자세한 내용은 [PDF] 다운로드.


표준정규분포의 pdf를 $\phi(x)$로 놓고 스케일링한 것을 $\phi_{\sigma}(x) = \phi(x/\sigma)/\sigma$라고 표기하면 편리하다. 이때 $N(\mu,\sigma^2)$의 pdf를 $\phi_\sigma(x-\mu)$로 쓸 수 있다.

표준정규분포 $\phi$의 도함수들은 $\phi$를 이용하여 표현할 수가 있다. 이 관계를 이용하면 미분 계산을 쉽게할 수 있다.
\begin{align}
\phi'(x) &= -x \phi(x)\\
\phi''(x) &= (x^2-1)\phi(x)
\end{align}

Convolution의 성질을 이용하면 복잡한 계산들을 매우 간단하게 할 수 있다. 우선 convolution은 다음과 같이 정의된다.
\[
(f*g)(x) = \int f(x-y) g(y) dy
\]
정규분포의 pdf에서 다음과 같은 성질을 확인할 수 있다.
\begin{align}
(\phi_{\sigma_1}*\phi_{\sigma_2})(x) = \phi_{\sqrt{\sigma_1^2+\sigma_2^2}}(x)
\end{align}
이 성질을 이용하여 $\int \phi^2(x) dx$, $\int \phi_{\sigma_1}(x)\phi_{\sigma_2}(x) dx$ 등을 쉽게 계산할 수 있다.

R 프로그래밍: R에서 C 함수 부르기

R 프로그래밍: R에서 C 함수 부르기

이 글에 사용된 코드와 명령은 MacOSX와 Linux에서 확인된 것이다. 특별히 설치되어야 하는 것은 없다. 기본적으로 R과 기본 개발도구가 설치되어 있으면 된다. MacOSX에서는 Xcode를 설치하면 된다. Fortran/C/C++ 무엇으로 짜야 하나 고민이 된다. 다음은 확인된 사실은 아니지만...

• R에서 vectorize할수 있는 데까지 다 했는데도 속도를 빠르게 하고 싶다면 Fortran/C/C++을 생각해야 한다.
• Fortran이 C보다 우월한가? 논쟁만 부르고 결론은 나지 않는 질문이다. 그럴 수도 있고 아닐 수도 있다.
• R에서는 공식적(?)으로 Fortran 77만 지원한다. Fortran 95로 작성한 코드를 R에서 부르는 방법이 없는 건 아니지만... 어떤 사람은 Fortran 77은 절대로 쓰지 말아야한다고도 한다. Fortran 95도 있고 최근 2010년에 표준화된 흔히 부르는 이름으로 "Fortran 2008"도 있다.

.Call

R 패키지를 개발할 때나 아니면 그냥 함수를 짤 때 C/C++/Fortran으로 짠 코드를 불러서 쓸 수 있다. R 함수 중 .C()는 컴파일된 C 코드를 호출할 때, .Fortran()는 컴파일된 Fortran 77 코드를 호출할 때 쓸 수 있다. 이것보다 좀더 유연하고 성능도 나은 것이 .Call()이다. 이것은 C/C++ 코드에 R 객체를 넘겨줄 때 사용된다. .External()도 마찬가지인데 .Call()은 고정된 수의 인자를 넘기고 .External()은 argument list로 넘긴다는 데에 차이가 있다. .Call()을 통해 사용할 C 코드를 작성할 때

• R.h, Rinternals.h 헤더를 반드시 include한다.
• R 객체를 받을 때 SEXP, REALSXP, INTSXP, STRSXP 타입을 이용한다.
• R 객체의 속성은 getAttrib()로 뽑아낼 수 있다.
• R 객체 안에 들어있는 데이터에 대한 포인터는 REAL(), INTEGER() 등을 이용하여 얻을 수 있다.
• R 객체를 만들 때는 allocVector(), allocMatrix()를 이용하여 SEXP 타입 객체로 만든다.
• 이때 항상 PROTECT()와 UNPROTECT()를 이용한다. 안전한 메모리 확보를 위한 장치이다.

예를 들어 다음과 같은 함수를 작성하고 myfun.c 파일에 저장하였다고 하자.

#include <R.h>
#include <Rinternals.h>

SEXP say(SEXP x) {
    return x
} 

컴파일은 다음과 같이 한다.

$ R CMD SHLIB myfun.c

컴파일하면 만들어지는 myfun.so를 R에서 불러 사용한다. R에서 다음과 같이 로딩한다.

> dyn.load("myfun.so")

로딩이 되었으면 다음처럼 실행할 수 있다.

> .Call("say", 3)
[1] 3

예제: 평균 계산하기

이제 C로 평균을 구하는 함수를 다음과 같이 작성하고 mystat.c 파일로 저장하였다고 하자.

// mystat.c
#include <R.h>
#include <Rinternals.h>

SEXP mean(SEXP x) {
  double *vec = REAL(x);
  int n = length(x);
  SEXP val = PROTECT(allocVector(REALSXP, 1));

  double sum = 0;
  for(int i=0; i < n; i++) sum += vec[i];
  REAL(val)[0] = sum / n;
  UNPROTECT(1);
  return val;
}

컴파일을 하면 mystat.so 파일이 생긴다.

$ R CMD SHLIB mystat.c

R에서 이것을 불러쓰기 위해 다음처럼 wrapper를 작성하자.

mymean <- function(x) {
    if(!is.loaded("mystat")) dyn.load("mystat.so")
    .Call("mean", as.numeric(x))
}

이제 R에서 mymean()이라는 함수를 사용하면 된다.

속도 비교

R에는 물론 평균을 계산해주는 mean() 함수가 있다. 만약 R에서 for으로 평균을 구하면 어떨까? iMac 3.06GHz Intel Core 2 Duo에서 1e7개의 난수를 더한 후 평균을 구하는 데에 다음처럼 오랜 시간이 걸린다.

> n <- 1e7
> x <- runif(n)
> system.time( {
+     total <- 0
+     for(i in 1:n) total <- total + x[i]
+     total / n
+ })
   user  system elapsed 
  8.436   0.076   9.165 

물론 R의 mean() 함수를 이용하면 빠르게 계산된다.

> system.time( mean(x) )
   user  system elapsed 
  0.079   0.001   0.086 

방금 C로 짠 함수의 경우에는:

> system.time( mymean(x) )
   user  system elapsed 
  0.020   0.000   0.023 

R의 내장 mean()보다 빠른 속도를 보여주고 있는데 데이터에 따라서 다를 수 있으니 항상 성능이 좋다고는 못한다. 어쨌든 R에서 for문 돌리는 것에 비하면 빠르다.

행렬

R에서 넘어온 행렬은 컬럼을 쌓은 벡터로 받아서 사용한다. $p\times q$ 행렬 $A$와 $q\times r$ 행렬 $B$의 행렬곱은 다음처럼 계산할 수 있다.

#include <R.h>
#include <Rinternals.h>

/**
 * @SmatA Sp by Sq matrix
 * @SmatB Sq by Sr matrix
 **/
SEXP prod(SEXP SmatA, SEXP SmatB, SEXP Sp, SEXP Sq, SEXP Sr) {

  int p = INTEGER(Sp)[0];
  int q = INTEGER(Sq)[0];
  int r = INTEGER(Sr)[0];

  // R matrix as a vector stacked column by column
  double* a = REAL(SmatA);
  double* b = REAL(SmatB);
  SEXP ans = PROTECT(allocVector(REALSXP, p*r));
  double sum;
  int i, j, k;
  for(i=0; i < p; i++) {
    for(k=0; k < r; k++) {
      sum = 0;
      for(j=0; j < q; j++) {
        sum += a[i + j*p] * b[j + k*q];
      }
      REAL(ans)[i + k*p] = sum;
    }
  }
  UNPROTECT(1);
  return ans;
}

R에서 parallel 패키지 이용하여 멀티코어 계산하기

R 3.0 기준으로 얘기하자면 parallel 패키지는 멀티코어 이용도 가능하게 해준다. multicore 패키지는 더이상 사용되지 않는다. 다음처럼 패키지를 로딩한다.

> require(parallel)

Linux, MacOSX에서 간단하게 mclappy와 같은 함수를 이용하는 예를 들면 다음과 같다. iMac 3.06GHz Intel Core 2 Duo에서 돌린 결과이다.

> system.time( out <- sapply(1:100000, function(.) mean(rnorm(100))))
   user  system elapsed
  2.916   0.039   2.959
> system.time( out <- mclapply(mc.cores=2, 1:100000, function(.) mean(rnorm(100))))
   user  system elapsed
  1.487   0.073   1.603

MS Windows에서는 mclapply가 제대로 작동하지 않는다. 대신 clusterApply를 쓸 수 있다. Intel Core i3 3.20GHz에 Windows 7이 설치된 PC에서 돌린 결과이다.

> system.time( out <- sapply(1:100000, function(.) mean(rnorm(100))) )
   user  system elapsed 
   2.66    0.00    2.67 
> cl <- makeCluster(4)
> system.time( out <- clusterApply(cl, 1:100000, function(.) mean(rnorm(100))) )
   user  system elapsed 
  12.85    8.23   21.94 
> stopCluster(cl)

이해가 가지 않는 결과이다. 코어 하나만 쓰는 것보다 코어 4개를 모두 쓰면 더 느려진다. 작업관리자에서 CPU 사용현황을 모니터링하면 멀티코어로 실행되는 것은 분명히 맞다. 이것은 core에 작업을 할당하는 문제 때문이다. 위 실험의 경우 100000개의 작업을 4개의 코어에 분배하여야한다. 그리고 하나의 작업은 매우 연산량이 매우 작다. 그 분배 처리 작업 때문에 하나의 코어에서 10만개 작업을 순차적으로 할당하는 것보다 느리게 되는 것이다. 코어가 4개이므로 10만번 시뮬레이션을 25000번의 시뮬레이션 4개로 나누어 할당하면 속도가 빨라진다.

> simul <- function(n) sapply(1:n, mean(rnorm(100)))
> system.time( out <- sapply(rep(25000, 4), simul) )
   user  system elapsed 
   2.60    0.01    2.68 
> cl <- makeCluster(4)
> system.time( out <-clusterApply(cl, rep(25000, 4), simul) )
   user  system elapsed 
   0.00    0.02    1.04 

조건부 확률의 조건부 확률. posterior predictive distribution.

조건부 확률

조건부 확률을 보통 $A$, $B$ 기호를 가지고 이야기하지만 여기서는 좀 헤깔릴 수가 있으니 $Y$, $A$로 표시해 보자. $Y$는 관심 사건이고 $A$는 주어진 조건에 해당하는 사건이다. 조건부 확률의 기본적인 정의에 따라 쓰면 다음과 같다.

\[
P(Y|A) = \frac{P(Y, A)}{P(A)}
\]


그리고 전확률 법칙(Law of Total Probability)라는 것도 있다. 표본 공간이 $B_1,B_2,\cdots,B_n$으로 분할(partition)되어있을 때
\[
P(Y) = \sum_{i=1}^n P(Y|B_i)P(B_i)
\]
이 성립한다는 법칙이다. 이것은 $P(Y)$를 구하기가 힘들 때 $P(Y|B_i)$들은 구하기 쉽고 각 $P(B_i)$들을 이미 알고 있는 경우에 유용하게 쓸 수 있는 법칙이다.

그런데, 이런 경우를 생각해보자. $P(Y|A)$를 구해야하는데 이것을 구하기가 힘들다. 그런데, $B_i$로 이 문제를 분할하여 구하는 것은 쉽다고 하자. 조건부 확률의 조건부 확률을 생각해보는 것이다. 다음처럼 구할 수 있다.
\[
P(Y|A) = \sum_{i=1}^n P(Y|B_i,A)P(B_i|A)
\]
이 식을 바로 생각하려고 하면 헤깔릴 때가 있다. $Y$를 $A$라는 조건 아래서 생각하고 있는 데 거기에 $B_i$라는 조건을 추가로 고려하면 어떻게 될까, 하는 것인데. 두번째 항을 $P(A|B_i)$인 것으로 착각을 할때가 있다. 식을 따져보면 당연한 것이기는 하다.
\[
\sum_{i=1}^n P(Y|B_i,A)P(B_i|A)
= \sum_{i=1}^n \frac{P(Y,B_i,A)}{P(B_i,A)}\frac{P(B_i,A)}{P(A)}
= \sum_{i=1}^n \frac{P(Y,B_i,A)}{P(A)}
= \frac{P(Y,A)}{P(A)}
= P(Y|A)
\]

Posterior Predictive Distribution

미지의 평균과 분산이 $\theta =(\mu, \sigma^2)$인 대상에 대해 $y=(y_1,\cdots,y_n)$을 관측했다. $\tilde y$는 다음 번 관측 대상의 값이다. 이때 $\tilde y$의 분포를 posterior predictive distribution이라고 한다. 다음처럼 구한다.

\begin{align}
p(\tilde y|y) &= \int p(\tilde y, \theta | y) d\theta\\
&= \int p(\tilde y|\theta, y)p(\theta|y) d\theta\\
&= \int p(\tilde y|\theta)p(\theta | y) d\theta.\\
\end{align}
마지막 행이 성립하는 것은 $\tilde y$와 $y$가 given $\theta$에서 조건부 독립이기 때문이다.

예제

여자아이가 태어날 확률을 $\theta$라고 하자. 서로 독립적인 $n$명의 신생아를 관찰하였을 때 여자아이의 수는 이항분포 $Binomial(n, \theta)$를 따른다. 사전분포로는 균일분포 $\theta \sim Unif(0,1)$을 사용하자. 1000명의 신생아를 관찰하였더니 485명이 여자아이였다. 1001번째 신생아가 여자아이일 확률은?

$n$명 중 여자아이의 수 $y$라고 하자. likelihood는
\[
p(y|\theta) = \binom{n}{y} \theta^y(1-\theta)^{n-y}
\]
piror distribution는
\[
\pi(\theta) = 1, \qquad 0 < \theta < 1
\]
posterior distribution은
\begin{align}
p(\theta|y) &\propto p(y|\theta)\pi(\theta)\\
&\propto \theta^y (1-\theta)^{n-y}
\end{align}
형태이므로 $Beta(y+1, n-y +1)$ 분포임을 알 수 있다.

이제 새로 태어날 아이에 대한 관측을 $\tilde y$라고 하면 여자아이일 확률은
\begin{align}
p(\tilde y = 1| y) &= \int_0^1 p(\tilde y = 1 | \theta,y)p(\theta|y) d\theta\\
&= \int_0^1 p(\tilde y = 1 | \theta)p(\theta|y) d\theta \\
&= \int_0^1 \theta p(\theta|y) d\theta\\
&= E(\theta|y) & (\theta|y \sim Beta(y+1, n-y+1))\\
&= \frac{y+1}{n+2} & (\because X\sim Beta(\alpha, \beta),\; EX = \alpha/(\alpha+ \beta))\\
\end{align}

따라서, 새로 태어날 아이가 여자아이일 확률을 486/1002.

베이즈통계 책

• Berger, J. O. (1985). Statistical Decision Theory and Bayesian Analysis. Springer- Verlag. 
• Ghosh, J.K. and Ramamoorthi (2002). Bayesian Nonparametrics. Springer. 
• Lange, K. (1998). Numerical Analysis for Statisticians. Springer. 
• Robert, C. P. and Casella, G. (1999), Monte Carlo Statistical Methods. Springer. 
• Norris, J.R. (1997). Markov Chains. Cambridge. 
• Chen. M.-H., Shao, Q.-M. and Ibrahim, J.G. (2000). Monte Carlo Methods in Bayesian Computation. Springer.  
• Gelman, A., Carlin, J.B., Stern, H.S. and Rubin, D.B. (2003). Bayesian Data Analysis. Chapman & Hall/CRC. 
• Ghosh, J.K., Delampady, M. and Samanta, T. (2006). An Introduction to Bayesian Analysis. Springer. 
• Congdon, P. (2003). Applied Bayesian Modelling. Wiley. 
• Congdon, P. (2005). Bayesian Models for Categorical Data. Wiley. 
• Nyhoff, L. R. and Leestma, S. C. (1996). Introduction to Fortran 90 for Engineers and Scientists. Prentince Hall.

선형모형 책

• [석사] Alvin C. Rencher and G. Bruce Schaaalje. Linear Models in Statistics (2nd ed.). 2008. Wiley. google books
• [대학원] Youngjo Lee, John A. Nelder, Yudi Pawitan. Generalized Linear Models with Random Effects: Unified Analysis via h-likelihood. Chapman & Hall/CRC. 2006. [google books]

확률론, 확률과정 교재

• [학부 확률과정] Rick Durrett. Essentials of Stochastic Processes. Springer. 2010. [google books]
• George F. Lawler. Introduction to Stochastic Processes (2nd ed.). Chapman & Hall/CRC. 2006. [google books] 
• Geoffrey Grimmett and David Stirzaker. Probability and Random Process (3rd ed.). Oxford Unisersity Press. 2001. [google books]
• Olav Kallenberg. Foundations of Modern Probability. Springer. 2002. [google books]
• Sidney I. Resnick. Adventures in Stochastic Processes. Birkhauser. 1992. [google books]
• Zdislaw Brzezniak and Tomasz Zastawniak. Basic Stochastic Processes. Springer. 2000. [google books]