2012년 8월 4일 토요일

통계의 기본 개념을 잘 설명한 유용한 사이트

스마트폰을 뒤적거리다가 우연히 발견하게 된 사이트이다.

품질공학 강좌(Taguchi methods)

타구찌 법(Taguchi methods)라는 품질공학 기법을 소개하는 사이트인데, 제 6편에서 통계의 기초에 대한 내용을 매우 명쾌하게 설명해 주고 있어서 소개하고자 한다. 내가 종종 들추어보는 Biostatistics 책의 서문에서도 "이 책은 오로지 sample mean에 대한 내용을 다루고 있다"고 하지 않았던가? sample mean으로부터 모집단의 특성을 유추하고자 하는 것, 이것이 바로 통계를 공부하고 활용하는 중요 목적 중 하나이다. 이를 위해 기본 개념을 확실하게 다져 두는 것이 중요한 것이다.

2012년 8월 3일 금요일

가설 검정의 기본-모평균에 대한 가설 검정

[미리 밝히는 글] 이 글은 대학을 졸업한지 20년이 넘었지만 아직도 통계에 대한 기본 개념이 부족한 본인의 무지를 일깨우기 위해 다시 공부를 한다는 자세로 작성하고 있는 것이다. 따라서 오류나 잘못 설명된 부분이 있을 수 있다. 혹시 인터넷 서핑 중에 이 글을 발견한 독자가 있으시다면, 그저 수학적 지식이 부족한 전통적인 생물학자가 작성한 에세이 정도로 생각해 주시길.

보통의 기초 통계학 교과서에서는 점추정, 구간 추정에 대한 설명이 나온 뒤 바로 가설 검정에 관한 부분이 이어진다. 종이에 종 모양의 그림을 그리고 수직선을 몇 개 그은 뒤 여기가 기각역이고 유의수준이 어쩌고... 기계적으로 암기하고 단순한 문제를 풀 수는 있다. 그런데 진정한 이해가 뒷받침되지 않으니 내 의식은 20년이 넘게 무지와 몰이해의 늪을 헤매고 있는 듯하다. 정말 부끄러운 일이 아닐 수 없다.

가설 검정의 이론에 대한 설명은 여러 가지 방법이 가능하다. 내가 갖고 있는 기초 통계학 책만 해도 대여섯권은 될 터이니. 이 글에서는 UC Berkeley의 Steve Selvin이 지은 BIOSTATISTICS: How It Works (Pearson 2004)의 187쪽부터 나오는 부분을 공부하면서 느낀 바를 요약하고자 한다.

통계적 가설(statistical hypothesis)을 아주 단순하게 정의한다면, sampled variable의 확률 분포에 대해 짐작하는 것을 의미하다. 예를 들자면 이런 것이다. "특정한 관찰값을 나오게 한 모집단의 평균값은 얼마이다" 이런... 여기에서 그리스 알파벳을 어떻게 적어야 하는 거지.

일반적으로 가설 검정은 귀무가설 H0과 대립가설 H1 중 하나를 고르는 과정이 된다. 물론 이 과정에는 오류가 있을 수 있다. H0가 옳은데도 기각하게 되는 것을 Type I error 혹은 false positive error라 한다. 반대로, H1이 옳은데도 이를 받아들이지 않는 오류를 Type II error 혹은 false negative error라 부른다. 그러면 무엇을 귀무가설로 설정할 것인가? 현재까지의 일반적인 믿음을 귀무가설로 하는 것이 정석이다. 예를 들어 고혈압 치료제의 효능을 검사하는 경우를 생각해 보자. 약을 투여한 사람에 대해서 혈압을 측정했더니 어떤 평균값이 얻어졌다고 하자. 이때의 귀무 가설은, 이들의 혈압은 투여하지 않은 상태의 평균치와 같다는 것이 귀무 가설이된다. 이 고혈압 환자의 수축기 평균 혈압은 140 mmHg이었다. 그런데 이들 중 20명을 취해서 약을 투여한 뒤 다시 혈압을 재어 평균을 내니 134가 나왔다고 하자. 분명 134는 140과 다른 수치이다. 그러나 모평균이 140인 집단에서 20명을 취하여 아무런 조치를 취하지 않은 상태에서 혈압을 재면 random variation에 의해 134라는 값이 나올 가능성은 얼마든지 있다. 따라서 이러한 차이는 약물에 의한 것이 아니라 정상적인 확률 분포에서 있을 수 있는 일이 벌어진 것이라고 결론을 내린다면, 결국 귀무 가설을 기각하지 못하는 것이다. 이번에는 약물 투여군의 평균 수축기 혈압이 131이라고 해 보자. 알려진 모평균보다는 좀 더 멀어졌다. 만약 이 표본 평균이 모평균에 비해 점점 멀어진다면, 이러한 차이는 확률 분포에 의해 자연스럽게 일어졌다고 보기는 점점 어려워진다. 차라리 약물투여군의 평균 혈압(모평균)이 아예 낮아졌다고 보는 것이 더 타당하게 된다. 이렇듯 알려진 모평균과 실제 측정한 표본평균 사이의 차이가 얼마나 멀어지면 귀무가설을 기각하는 것이 타당할까? 이에 대한 기준을 제시하는 것이 바로 유의수준(1%나 5%와 같은 확률로 주어짐)이고, 귀무가설을 기각할 수 있는 측정값의 경계가 되는 값이 바로 critical value(임계값)이다.

유의수준이란, 귀무가설이 참인데도 불구하고 대립가설을 채택하게 될 확률이다. 즉 유의수준이란 Type I error가 발생할 확률을 의미하며, 1 - (유의수준) = (신뢰도)이다. 다시 한 번 위의 예제를 들어 설명해 보자. 이 고혈압 환자의 평균 수축기 혈압은 140이고, 정규분포를 한다고 가정하면 이론적으로는 30이나 350과 같은 극단적인 수치가 나올 가능성이 있다. 그러나 이런 값이 나올 가능성은 매우 낮다. 100번에 1번(쉽게 말해 1%) 나올 정도로 모평균에서 먼 값이 나온다면, 차라리 모평균 자체가 다르다고 가정하는 것이 낫다. 여기서의 1%가 바로 유의수준이 되는 것이다.

여기서 또 부끄러운 고백을 아니할 수 없다. 가설 검정 문제에서 표본 평균이 항상 나오게 되므로, 나는 그동안 표본 평균으로부터 모평균이 포함될 신뢰구간을 구해 놓고서는 귀무가설에서의 모평균이 이 구간내에 있는지 혹은 바깥에 있는지를 가지고 가설을 검정하려고 했었다. 그런데 지금 와서 기초 통계 교과서를 살펴보면 이런 방식으로 접근하는 것이 보이지 않는다. 아마도 대부분의 교과서에서 신뢰구간 다음에 가설 검정이 나오기 때문에 이런 오해를 하게 된 것 같다. 이러한 방식의 풀이가 수학적으로 정말 맞는지를 증명할 실력은 되지 않는다. 일단은 잘못된(?) 믿음을 버리고, 일반적인 풀이의 방법으로 되돌아가도록 하자.

그러면 다시 Selvin의 책 190쪽에 나온 예제로 돌아가 보도록 하자. 어떤 집단의 혈압 모평균이 130 mm Hg이고 분산은 400으로 알려져 있다. 새로운 약을 투여했더니 16명의 무작위 표본으로부터 평균 120 mm Hg의 혈압 측정값을 얻을 수 있었다. 이 약은 정말로 혈압을 낮추는 효과가 있을까?

H0 : 혈압은 여전히 130이다(약물의 효과가 없다)
H1 : 혈압은 120이다(130보다 낮다).

표본평균이 130에 충분히 가깝다면 귀무가설을 채택해야 하고, 상당히 멀다면(130보다 상당히 작다면) 대립가설을 채택해야 한다. 그렇다면 가까운 정도를 정의해야 할 것이다. 첫번째 과정에서는 귀무가설을 기각시킬 수 있는 표본 평균의 값 c를 결정해야 한다. 이때 모든 c의 범위를 기각역이라고 한다. 본 예제에서는 표본평균 < c 이면 귀무가설을 기각하도록 하는 c를 구해야 한다. 임계치는 Type I error가 충분히 작아지도록 만드는 값으로부터 결정한다.


(흐유, 한컴오피스 한/글에서 수식 입력하느라 힘들었다) 위의 식에서 c의 값을 계산할 수 있다. 실제 계산을 해 보면 c = 121.8이다. 만일 표본평균이 이보다 작으면, 귀무가설은 기간한다. 단측검정 혹은 양측검정의 여부는 대립가설의 형태에 따라 정해진다. 

Type I error의 확률에 기반하여 기각역을 결정했다면, 다음으로는 Type II error를 생각해 보도록 하자. Type II error는 어떤 경우에 발생하는가? 모평균이 변해서 이에 따라 표본평균도 다르게 나왔는데, 이를 random variation에 의해 차이가 발생했다고 생각하고 귀무가설을 그대로 인정한다면 바로 Type II error를 범하게 되는 것이다. 

2012년 8월 2일 목요일

통계는 어렵다?

참 이상한 일이다. 대학과 대학원때 각각 3학점씩 하는 [확률과 통계] 과정을 두 번이나 수강하였고, 전부 A학점을 받은 것으로 기억하는데도, 왜 biodata의 분석에 써먹을 정도의 기본 개념이 남아있지 않은 것일까? 단순한 망각일까, 시험 문제 풀이는 겨우 했지만 제대로 된 개념을 결국은 잡지 못한 것일까, 아니면 실제 데이터를 가지고 충분한 연습을 하지 않았기 때문일까? KOBIC에서 열리는 제10차 차세대 생명정보학 교육 워크샵을 수강하면서, 그나마 실낱과 같이 남아있는 용어들의 정의에 대한 기억을 되살릴 수 있음을 감사하고 있다.

KRIBB에 와서도 부족한 개념을 보충하기 위해 KAIST 산업공학과 채경철 교수님의 통계학 수업을 잠시 청강한 적이 있다. 한 학기를 제대로 채우지 못했지만, 이 과정을 통해서 남긴 가장 중요한 '개념'은 "확률변수란 표본점에 실수를 대응시키는 함수"라는 것이다. 변수라는 이름이 붙어 있기에 오해를 하기 쉬운데, 실제 일어날 수 있는 모든 이벤트라는 정의구역에 대해 '실수'를 대응시키는 함수이다. 쉬운 예를 들자면, 20마리의 실험용 생쥐에 독성 물질을  주사하여 테스트하는 경우를 생각해 보자. 어떤 물질을 투여했을때 죽은 쥐의 수가 바로 확률 변수의 한 예에 해당한다.

확률변수가 함수라는 개념을 확실히 갖고 있지 않으면, P(X=x)라는 표기를 이해하기 어렵다. 그냥 p(x)라고 하면 되는 것 아닌가? 그렇지 않다. X는 확률 변수(함수)이고, 이 함수의 값이 x(즉 X = x)일 경우에 해당하는 확률을 P(X=x)로 표시하는 것이다. 만약 이 확률변수가 f(x)라는 분포함수를 갖고 있다면, P(X=x) = f(x)라고 표시할 수 있겠다. 오늘 UNIST 남덕우 교수의 강의에 따르면 확률에서 대문자는 변수요, 소문자는 상수라고 하였다. 예외가 하나 있는데, 그것은 N(모집단의 크기)과 n(표본의 크기)이다.

다음으로 이해하기 어려운 것은 표본 표준편차를 구할 때 표본 크기(n)가 아니라 n - 1로 나누는 것이다. 이것은 어느 기초 통계학 책에서나 충분히 설명을 하고 있으니 걱정할 것은 없다.

부끄러운 이야기지만 아직까지도 나를 혼동스럽게 하는 것은 가설 검정과 구간 추정 부분이다. 양측 검정과 단측 검정을 각각 어떤 경우에 해야 하는지도 아직 혼동할 정도이니... 정답을 먼저 이야기하자면 대립가설이 어떤 형태냐에 따라 달라지게 된다.

우선 중요한 숫자 몇 가지만 기억하도록 하자. 정규분포에서 평균을 중심으로 하여 +/-1 표준편차 이내의 구간에는 전체 데이터의 68.3%가 존재한다. +/- 2 표준편차에는 95.4%,  +/- 3 표준편차 내에는 99.7%가 존재한다. 이를 68-95-99.7 법칙이라 한다.

자, 그러면 표본평균에서 추정한 모평균의 신뢰구간을 생각해 보자. 표본의 수가 충분히 크면, 모집단이 어떤 분포를 따르든 상관없이 표본평균은 정규분포를 따르게 된다. 이론적으로  정규분포는 극단에 치우친 어떤 값도 가질 수 있지만, 평균에서 지나치게 멀리 떨어진 값이 나올 확률은 매우 적어진다. 따라서 이러한 극단의 값이 나오게 된다면(기각역에 포함될 경우) 차라리 평균의 위치를 옮기는 것이 더 타당하게 된다.  그러면 도대체 평균에서 얼마나 멀어야 극단으로 정의할 것인가? 이를 정의하는 것이 바로 신뢰도와 오차율(alpha)이다.

신뢰구간은 매우 정확하게 이해해야 한다. 모평균은 이미 결정되어 있는 숫자이다. 신뢰구간은 한 번의 표본추출(표본 크기가 1이라는 뜻이 아니다!!)에서 얻은 표본 평균에서 계산으로 얻어진다. 표본 평균은 매번 다르게 나오므로, 신뢰구간의 모습도 매번 달라질 수 밖에 없다. 표본 추출을 100번해서 모평균의 구간을 추정했더니, 그중 95번은 신뢰구간 안에 실제 모평균이 들어가게 된다 - 이것이 바로 95% 신뢰구간의 정확한 의미이다.

2012년 8월 1일 수요일

[Perl programming #2] substr() 함수 알뜰하게 사용하기

스트링 형태로 저장된 염기 혹은 단백질 서열을 60 문자마다 깨끗하게 줄바꿈을 하여 FASTA format으로 출력할 일이 아주 많다. BioPerl을 쓰면 되지만, 일부러 sequence object를 선언하여 쓰기에는 너무 번거롭다. 그동안 내가 사용한 코드는 스트링을 낱자로 하나하나 분해하여(split() 사용) 60개가 맞는지 센 다음에 출력하는, 그야말로 원시적인 것이었다. 작동하는데는 문제가 없으나, 여기에 공개하기에는 너무나 부끄러운 코드이다. 어제 Python 교육을 받으면서 머리를 스쳐가는 것이 있어서 다음과 같이 응용해 보았다.

$orig = 'AGCTAGGC...';
while ($orig) {
    print substr $orig, 0, $width;
    print "\n";
    $orig = substr $orig, $width;
}

대단히 간결하고 명확하지 않은가? while 구문을 이렇게 활용할 수 있을 것으로는 생각하지 못했다. 짧고, 읽기 쉽고, 능률적인 코딩 습관은 항상 추구해야 할 바라고 생각한다. 물론 한 가지의 목적을 달성하기 위하여 여러가지 해법이 있을 수 있다는 것인 Perl의 철학이지만, 조금만 생각을 더 하고 궁리를 해 보면 훨씬 능률적인 작업을 할 수 있음을 알아야 한다.

substr() 함수의 다양한 사용법은 여기에서...

http://perldoc.perl.org/functions/substr.html

아직 다 끝난 것이 아니다. substr()는 lvalue로 쓸 수 있다는 것. 이는 등호(=)의 왼쪽에 올 수 있음을 의미하는 것이다. 일반적인 함수라면 반환된 값을 내어 놓고, 이를 등호 왼쪽에 오는 변수에 저장하게 된다. substr()도 마찬가지이다. 인수로 주어진 스트링 내부의 지정된 부분을 끊어서 왼쪽의 변수에 자정하는 것이다.

$var = my_func(..)

substr()를 등호 왼쪽에 쓰게 되면, 얘기가 달라진다. substr()을 예로 든다면, 인수로 지정된 변수 내의 특정 위치를 오른쪽 변수에 저장된 값으로 바꿔친다는 것이다. 위에서 보인 샘플 코드를 더욱 단순하게 만들 여지가 있다는 것이다.

그런데 아무리 궁리를 해도 답이 안나온다^^ 위의 코드가 가장 간단한 코드일지도 모르겠다.

2012년 7월 31일 화요일

Python 공부하기

KOBIC에서 주최하는 제10회 차세대 생명정보학 교육 워크샵에서 이틀째 날을 보내고 있다. 어제의 리눅스 기초에 이어서 오늘은 파이썬 프로그래밍이다. Perl로 굳어진 사고 체계에 파이썬을 끼워 넣는 것은 쉽지 않다. 내장 함수를 써야 하나, 메쏘드를 써야 하나? 변수명 앞에 $ 혹은 @가 붙지 않으니 이것이 예약어인지 타당한 변수명인지 알 수가 없다.

Mutable sequence를 다루고 조작하는 방법은 매우 인상적이다. 원래 파이썬에 관심을 갖게 된 것은 남들이 만들어 둔 스크립트를 이용할 때 설정 등의 문제로 에러가 난 것을 다루기 위한 가벼운 마음으로 시작한 것인데, Perl로는 하기 어려운(사실 불가능한 것이 있으리라고는 생각하지 않지만) 작업을 능률적으로 하는데 도움을 줄 것으로 보인다. 어쩌면 파이썬과 Perl을 비교하면서 Perl에 대한 사랑이 더 깊어질지도^^

인실리코젠의 김형용씨가 강의와 실습을 진행하느라 하루 종일 수고를 하고 있다.

2012년 7월 12일 목요일

[Perl programming #1] my()를 이용한 변수 선언

[들어가는 글]

펄(Perl)을 업무에 활용한지도 꼭 12년째가 되어 간다. 내가 펄을 처음으로 공부하면서 접한 교재는 인터넷에서 찾아서 인쇄한 Russell Quong의 "Perl in 20 pages" version 2000c와 "Robert's Perl Tutorial Version 4.1.1"이다. 이 문서들이 아직도 인터넷 공간에 살아 있는지는 모르겠다. Univ. of Washington의 Phred/Phrap/Consed에 포함되어있는 스크립트가 조금씩 펄 스크립트 작성 실력을 키우면서 정말 큰 도움이 되었던 재료들이었다. 이제는 [펄 쿡북]을 참조하면서 필요한 부분은 그때마다 참조하고 있지만, 알고리즘이나 데이터 구조에 대한 깊이 있는 지식 없이 그저 필요한 기능만 구현하는데 집중하다 보니 항상 아쉬움이 많았다. 이러한 아쉬움을 달래고자, 조금 더 진지하게 펄을 공부해 보기로 하였다. 지금부터 내가 추구하는 것은 깔끔하고, 문법적으로 정확하며, 읽기 쉽고, 재사용이 가능한 펄 코드를 짜는 것이다.

[변수 초기화]

나에게는 아주 나쁜 버릇이 있다. 변수를 선언하는데 매우 불성실하다는 것이다. use strict나 my()라고 몇 자 더 치는 것을 게을리하여 오류가 난 코드를 고치는데 무척 애를 먹기도 한다. 변수명을 잘못 기입하면 오류의 원인을 찾기가 매우 어렵다.

#!/usr/bin/perl

use strict;
use warnings;

이렇게 스크립트를 시작하는 것과,

#!/usr/bin/perl -w

use strict;

이렇게 쓰는 것은 본질적으로 같(을 것이)다.

use strict 디렉티브를 사용하면, 모든 변수는 my()를 써서 선언해야 한다. 변수의 선언과 초기화는 다르다. 괄호를 언제 써야 하는가? 다음의 사례를 보라.

my $foo;  # 이것은 선언만 하는 것
my $foo = 123; # 선언과 초기화를 동시에
my ($foo, $bar); # 여러 변수를 한번에 선언

초기화하지 않은 변수가 갖는 값은 undef이다. 그러나 이는 $var eq undef 처럼 비교할 수 있는 것이 아니다. 그리고 undef는 엄연히 말하자면 null character 혹은 숫자로서의 0과는 다르다.

변수의 선언과 초기화가 다르다는 것을 다음의 스크립트에서 확인해 보자.


#!/usr/bin/perl

use strict;
use warnings;

my $foo;

if (defined $foo) {
    print 'Var $foo is defined!', "\n";
} else {
    print 'Var $foo is not defined!', "\n";
}

상세한 설명은 Perl 공식 문서의 Private Variables via my () 섹션을 참조하자. 참고로 이 마지막 문장은 이 글을 쓴지 10년이 지난 2022년 3월 13일에 추가하였다.

2012년 7월 9일 월요일

Consed 23.0와 DISPLAY 문제

3개월만의 포스팅이다. 주로 네이버 블로그(http://jeong_0449.blog.me/)에 글을 올리고 있고, 최근에는 genoglobe.com 및 genoglobe.net이라는 도메인을 구입하여 여기에도 글을 쓰고 있으니 당연히 구글 블로거는 소홀해질 수 밖에...

올해에는 consed version 22.0과 23.0이 근소한 시간차를 두고 발표되었다. 대부분의 성능 개선은 version 22.0에서 이루어졌고, 사소한 bug의 수정과 BamView의 성능 개선이 있었다고 한다. 22.0부터는 installConsed.perl이라는 친절한 설치 스크립트까지 포함되었다. mktrace, phd2fasta 및 454 데이터 관련 바이너리를 알아서 컴파일하여 설치해 주는 것으로 보인다. Consed의 home은 여전히 /usr/local/genome이다.

만약 consed_linux64bit를 실행시키는데 libstdc++.so.5가 없다는 에러 메시지가 나오면, compat-libstdc++-33 패키지를 설치하면 된다.

consed 22.0 or later와는 직접 관계가 있는 일인지 혹은 아닌지 잘 모르겠지만, CentOS 6.2가 설치된 컴퓨터에 consed를 깔고 이를 다른 컴퓨터에서 접속하여 실행하려고 하면 꼭 다음과 같은 에러가 난다. 다른 리눅스 서버에서 ssh -X로 접속하든, 혹은 윈도우 컴퓨터에서 Xmanager로 접속하든 같은 문제가 발생한다.

 CentOS 5.x이나 우분투에서는 절대로 아무런 문제 없이 실행이 된다.


[hyjeong@proton edit_dir]$ consed
no ~/consedrc file so no user resources will be used--that's ok
no consedrc file so no project-specific resources--that's ok
couldn't open readOrder.txt--that's ok
_X11TransSocketINETConnect() can't get address for localhost:6015: Name or service not known
Error: Can't open display: localhost:15.0

처음에는 이것이 Xmanager의 문제인줄 알고 넷사랑 컴퓨터에 질문을 하기도 했었다. 오늘에서야 부분적인 답을 알게 되었다. 

export DISPLAY=127.0.0.1:15.0

이렇게 한 다음에 consed를 실행하면 된다. consed가 설치된 컴퓨터에서 localhost의 위치를 찾지 못해서 생기는 문제라고 하는데, 왜 이런 문제가 발생하는지는 도무지 알 길이 없다. 게다가 이래가지고는 consed를 한번 실행해서 접속시마다 바뀌는 스크린 번호를 에러 메시지와 함께 알기 전에는 DISPLAY 변수를 설정하지 못한다. 

/etc/hosts내에서 localhost가 잘못 정의된 것일까? CentOS 6.2를 설치한 직후의 상태는 다음과 같고, 다른 host를 추가한 것 이외에는 건드린 것이 없다.

127.0.0.1 localhost localhost.localdomain localhost4 localhost4.localdomain4
::1 localhost localhost.localdomain localhost6 localhost6.localdomain6

어쩌면 consed 프로그램이 /etc/hosts 파일을 제대로 파악하지 못하는지도 모른다는 생각이 들었다. /etc/hosts를 다음과 같이 심플하게 변경해 보았다.

127.0.0.1   localhost

그랬더니... 잘 되네? 서버 컴퓨터의 OS를 CentOS 6.2로 업그레이드한 뒤 몇달간 골머리를 앓게 했던 문제를 이렇게 허무하게 해결하다니...