Ensemble method: Bagging (bootstrap aggregating)

일반적인 ensemble 방법의 하나가 bootstrap aggregating인 Bagging 방법이다.

Leo Breiman에 의해서 1994년에 만들어 졌다.

bagging은 오리지널 트레이닝 셋으로 부터 bootstrap 샘플링 방법을 이용해서 여러개의 트레이닝 데이터를 생성하는 방법이다.

이렇게 생성된 데이터 셋은 하나의 모델을 생성하기 위해서 사용 된다.

상대적으로 Ensemble 기법들 중에서는 Bagging은 단순한 방법에 속하지만 꽤 잘 동작 한다.

특히 Decision Tree 알고리즘의 경우 약간의 입력 데이터의 변화에도 급격히 다른 결과값을 이끌어 내므로

Bagging으로 이것을 보완 할 수 있다.

ipred package는 bagged decision tree를 구현한 것이다.

resubstitution error는 한글로 재치환오류라고 하며 트레이닝 셋에서 훈련된 모델에서의 에러를 말한다.

다른 용어로는 In sample error라고도 한다.

당연히 수치는 긍정적으로 평가된다. 왜냐하면 같은 데이터로 훈련하고 테스트 하기 때문이다.

Bootstrap sampling

큰 데이터 셋을 추정하기 위해서 임의로 여러번 샘플링을 하는것을 말한다.

When this principle is applied to machine learning model performance, it implies the creation of several randomly-selected training and test datasets, which are then used to estimate performance statistics.

The results from the various random datasets are then averaged to obtain a final estimate of future performance.

어떤 점이 그럼 k-fold Cross Validation이랑 다른 것일까?

이 방법은 생성한 각각의 구획들이 단 한번만 나타나게 된다.

하지만 bootstrap의 방법은 복원 추출을 하기 때문에 여러번 나올 수 있다.

증명에 의하면 초기 데이터 셋의 수인 N이 엄청 크다면, bootstrap으로는 63.2%를 커버한다는 증명이 있다.

결국 36.8%는 트레이닝 셋이다.

original data의 size가 N 레코드라고 하자. 이때 bootstrap으로 샘플링을 수행하면 N번 수행하면 당연히 N개의 레코드가 있고

이러한 N개의 레코드는 원래 original data size인 N개 중 63.2%를 포함한다고 한다.

N번의 bootstrap sampling에서 단일 record가 선택 되어질 확률은 아래와 같다.

$$1-(1-\frac{1}{N})^{N}$$

이것은 asymptotically approaches $1-e^{-1}=0.632$에 근접한다.

이렇게 bootstrap sampling으로 포함되지 않은 것들은 test set으로 자동으로 포함된다.

이에반에 10-fold cv를 할경우 90%가 결국 트레이닝이 된다.

N이 충분히 크다면 bootstrap은 더 낮은 representative를 보인다.

당연히 63.2%의 트레이닝 셋을 가지고 테스트 하기 때문에 더 큰 트레이닝셋을 가지고 만든 모델보다 더 낮은 정확도를 가지게 된다.

최종적인 error rate의 계산식은 아래와 같다.

$$ error = 0.632 \times error_{test} + 0.368 \times error_{train}$$

cross validation 보다 더 좋은 bootstap의 장점은 더 작은 dataset에서 잘 동작하는 특성이 있다는 것이다.

또한 bootstrap sampling은 모델 정확도 향상의 효과가 있을 수도 있다

주어진 데이터로 부터 복원 표본을 구하는 작업을 여러 번 반복해 원하는 값을 추정 하는 것이다.

예를 들어 보면,

평균의 신뢰 구간을 구하는 경우를 생각해보자.

평균의 신뢰 구간은 평균이 속할 범위를 95% 신뢰도로 찾는 것이다.

부트스트래핑은 데이터로 부터 표본을 복원 추출로 구하고 이들의 평균을 기록하는 일을 반복한다.

이렇게 구한 평균들을 나열한 뒤 가장 작은 값으로 부터 2.5% 가장 큰 값으로 부터 97.5% 지점의 평균 두 개를 구한다.

그러면 두 평균 사이의 구간이 평균 95% 신뢰 구간이 된다.

실제 코드로 구현해 보자.

평균 30, 표준편차 3인 정규 분포로부터 1,000개의 난수를 만들고 이 난수로부터 역으로 평균을 추정해보자.

x <- rnorm(1000, mean=30, sd=3)
t.test(x)
	One Sample t-test
data:  x
t = 308, df = 999, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 29.72329 30.10447
sample estimates:
mean of x 
 29.91388 

t검정 결과 평균의 신뢰 구간이 (29.78891, 30.16317)로 구해졌다. 

library(foreach)
bootmean <- foreach(i=1:10000, .combine=c) %do% {
    return(mean(x[sample(1:NROW(x), replace=TRUE)]))
}

이제는 부트스트래핑으로 복원 추출을 통해 샘플을 구하고 평균을 구하는 작업을 10,000회 반복해보자.

평균의 신뢰구간은 이 평균들의 하위 2.5%, 상위 97.5%에 속하는 값의 구간이다.

bootmean <- sort(bootmean)
bootmean[c(10000*0.025, 10000*0.975)]
[1] 29.72578 30.10531

t검정과 유사한 값이 부트스트랩핑으로도 구해진것을 알 수 있다.

부트스트래핑을 사용한 모델의 정확도 평가도 마찬가지로 수행된다.

데이터의 일부를 트레이닝 데이터로, 일부를 테스트 데이터로 복원 추출을 통해서 구한다.

그리고 이를 수회 반복하여 정확도를 추정한다. 이것이 trainControl()에서 method="boot"가 하는 일이다.

.

주성분 분석은 데이터에 많은 변수가 있을 때 변수의 수를 줄이는 차원 감소(Dimensionality Reduction)기법 중 하나다.

PCA는 변수들을 주성분(Principal Component)이라 부르는 선형적인 상관관계가 없는 다른 변수들로 재표현한다.

무슨말인가 하면,

선형적인 관계가 있는 두 변수가 아래와 같이 있다고 가정해보자.

A = 2c+3

B= 3c

이다.

이때 Y = f(A,B)라는 어떤 모델이 있다고 하자.

이럴때 실제로 A,B 둘 사이에는 C라는 변수에 의한 선형적 관계가 있는데 둘 다 이용해서 Y를 예측하는 모델을 생성할 필요가 있을까?

이와 같을 때는 f(A,B)로 Y를 예측하기 보다는 그냥 f(C)로 Y를 예측하는것이 더 좋다.

Principal component들은 원 데이터의 분산(퍼짐 정도)을 최대한 보존하는 방법으로 구한다.

좀 더 자세한 내용은 아래의 이전 포스트를 참조하자.

practical machine learning (johns hopkins)

machine learning (standford)

R을 이용해서 실제로 적용해 보기 [R프로그래밍 실습서]

pincomp() 함수 또는 preComp()를 이용한다.

pincomp: 주성분 분석을 수행 한다.

아래와 같이 1:10을 저장한 x

약간의 noise를 추가한 y. 즉 x + noise = y 이다.

마지막으로 z는 x + y 에 noise를 추가한 것이다.

data frame으로 내용을 출력해보면 아래와 같다.

위 데이터들은 주성분들이 원 데이터의 분산 중 얼마만큼을 설명해 주는지를 알 수 있다.

Proportion of Variance 행을 보면 첫 번째 주성분은 데이터의 분산 중 92.86%를 설명해 주며,

두 번째 주성분은 데이터의 분산 중 7.1%를 설명함을 알 수 있다.

세 번째 주성분은 가장작은 0.03%의 분산을 설명 한다.

마지막 행의 Cumulative Proportion은 Proportion of Variance의 누적 값이다.

결국 원 데이터의 분산은 첫 번째와 두 번째 주성분에 의해 99.97%가 포함됨을 알 수 있다.

이들 두 주성분상의 좌표는 scores를 보고 구하면 된다.

결국 x,y,z의 데이터 분포는 그냥 2개의 차원으로 축소가 된것이다. 아래와 같은 좌표값을 가지게 된다.

알츠하이머(Alzheimer)를 이용한 예제 [Practical Machine Learning HW-2]

데이터를 불러온다.

위 데이터 셋에서 "IL"로 시작하는 트레이닝 셋 variables(변수)들에 대해서만 PCA를 수행한다.

preProcess()를 이용하며, PCA 수행결과 80%의 variance(분산)를 capture할 수 있는 변수들의 갯수를 결정해 보자.

먼저 grep과 정규표현으로 해당 열 데이터를 수집한다.

아래와 같이 caret package의 preProcess()를 이용해서 pca를 수행한다. 임계값은 80%이다.

실행 결과 총 12개의 IL 변수중 7개만으로 원래 데이터 분포의 80%수준의 분산을 설명 할 수 있음을 알 수 있다.

이전과 같이 princomp를 수행하면 아래와 같다.

summary(princomp(training[, IL_str]))

PCA가 정확도에 미치는 영향을 glm 모델을 만들어서 알아보자 [Practical Machine Learning HW-2]

PCA를 적용하지 않은 glm 모델

PCA를 적용해서 수행한 결과이다.

적용전 0.65

적용후 0.72

결국 dimentionality reduction이 정확도 향상에도 기여한다는 것을 알 수 있다.

TensorFlow 설치법과 PyDev 설치법은 이전에 다뤘다.

원할한 개발을 위해서 PyDev에 TensorFlow를 연결하는것을 다룬다.

우선 필자는 TensorFlow를 Virtualenv에다가 설치 했다. VirtualEnv를 설치하는 이유는 아래와 같다.

PyDev의 VirtualEnv 설정 연결

PyDev와 TensorFlow 그리고 virtualenv가 모두 설치되었어야 한다.

Eclipse 메뉴에서

Preferences -> PyDev > Interpreters > Python 메뉴에서

New를 눌러서 Browse를 통해서 설정한 virtualenvs를 등록해준다.

등록할때 모든 종속 라이브러릴 같이 등록해준다. 팝업이 나오면 다 선택된 상태로 등록하면 된다.

프로젝트 생성 및 "Hellow Tensor"실행

아래와 같이 프로젝트를 생성 할 때 interpreter를 생성한 virtualenv용 python 환경으로 잘 선택해 줘야 한다.

그다음 Hello tensorflow 코드를 실행하면 정상적으로 실행 되는것을 알 수 있다.

'''
Created on Nov 17, 2015

@author: root
'''
import tensorflow as tf
hello = tf.constant('Hello, TensorFlow!')
sess = tf.Session()
print sess.run(hello)

a = tf.constant(10)
b = tf.constant(32)
print sess.run(a+b)

IntelliJ 기반의 Pycharm

설치방법: 이전포스트

실행 모습

아래와 같이 Tensorflow가 있는 virtualenv를 interpreter로 설정을 해주어야 한다.

웹기반 Jupyter로 실행하기

설치방법: 이전포스트

2015.11.9일에 공개된 오픈프로젝트 TensorFlow는 Google에서 공개한 DeepLearning을 위한 2세대 system이다.

TensorFlow의 내용은 이전 Post를 참조 한다.

이 시스템의 설치와 실행을 다룬다.

Python VirtualEnv에 기반한 설치 방법

필자가 사용하는 Desktop PC의 환경은 아래와 같다.

설치환경

필자는 VirtualEnv-based Installation을 수행한다(VirtualEnv 설치).

VirtualEnv 환경을 python 자체가 지원하므로 고립된 수행 환경을 생성할 수 있는 이점이 있다. 문제 발생시에도 모든 사람이 같은 상황이므로 해결방법을 그대로 적용할 수 있는 장점이 있다.

#Next, set up a new virtualenv environment. To set it up in the directory ~/tensorflow, run:
$ virtualenv --system-site-packages ~/tensorflow
$ cd ~/tensorflow

tensorflow는 이제 가상환경이 저장되는 디렉터리가 된다.

--system-site-packages 옵션의 의미는 가상환경에게 global site-package의 접근 권한을 허용 하는 것이다.

아래의 경로에 해당 가상환경에만 적용되는 site-package들이 저장되어 있다.

아래의 경로에는 해당 가상환경에서 사용 할 수 있는 python의 종류이다.

virtualenv를 활성화 시킨다.

$ source (설치경로)/bin/activate  # If using bash
source ./tensorflow/bin/activate #나의 경우 

비활성화는 deactivate라고 치면된다.

만약 csh 쉘을 쓴다면 아래의 명령어이다.

$ source bin/activate.csh  # If using csh

최종적으로 아래처럼 프롬프트 창의 이름 앞에 (tensorflow)가 들어가는 것을 볼 수 있다.

(tensorflow)root@jemin-virtual-machine:

실행된 virtualenv의 환경에서 TensorFlow를 설치한다.

설치 종류는 CPU-only version과 GPU-enabled version이 있으므로 자신이 원하는것을 선택하자.

필자는 일단 CUDA sdk를 설치하지 않았으므로 CPU-only를 설치했다.

# For CPU-only version
$ pip install https://storage.googleapis.com/tensorflow/linux/cpu/tensorflow-0.5.0-cp27-none-linux_x86_64.whl

# For GPU-enabled version (only install this version if you have the CUDA sdk installed)
$ pip install https://storage.googleapis.com/tensorflow/linux/gpu/tensorflow-0.5.0-cp27-none-linux_x86_64.whl

이렇게까지 하면 tensorFlow는 설치가 된것이다.

이제 간단한 예제를 실행해보자.

Hello, TensorFlow

Hello, TensorFlow
Python 2.7.6 (default, Jun 22 2015, 17:58:13) 
[GCC 4.8.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import tensorflow as tf
>>> hello = tf.constant("Hello, TensorFlow!")
>>> sess = tf.Session()
I tensorflow/core/common_runtime/local_device.cc:25] Local device intra op parallelism threads: 8
I tensorflow/core/common_runtime/local_session.cc:45] Local session inter op parallelism threads: 8
>>> print sess.run(hello)
Hello, TensorFlow!
>>> a = tf.constant(10)
>>> b = tf.constant(32)
>>> print sess.run(a+b)
42
(tensorflow)root@jemin-virtual-machine:~/tensorflow# 

# Make 100 phony data points in NumPy.
x_data = np.float32(np.random.rand(2, 100)) # Random input
y_data = np.dot([0.100, 0.200], x_data) + 0.300

# Construct a linear model.
b = tf.Variable(tf.zeros([1]))
W = tf.Variable(tf.random_uniform([1, 2], -1.0, 1.0))
y = tf.matmul(W, x_data) + b

# Minimize the squared errors.
loss = tf.reduce_mean(tf.square(y - y_data))
optimizer = tf.train.GradientDescentOptimizer(0.5)
train = optimizer.minimize(loss)

# For initializing the variables.
init = tf.initialize_all_variables()

# Launch the graph
sess = tf.Session()
sess.run(init)

# Fit the plane.
for step in xrange(0, 201):
    sess.run(train)
    if step % 20 == 0:
        print step, sess.run(W), sess.run(b)

# Learns best fit is W: [[0.100  0.200]], b: [0.300]

I tensorflow/core/common_runtime/local_device.cc:25] Local device intra op parallelism threads: 8
I tensorflow/core/common_runtime/local_session.cc:45] Local session inter op parallelism threads: 8
Hello, TensorFlow!
42
0 [[ 0.3138563   0.52770293]] [-0.0221094]
20 [[ 0.17014106  0.30057704]] [ 0.20599112]
40 [[ 0.12187628  0.23032692]] [ 0.27126268]
60 [[ 0.10676718  0.20918694]] [ 0.29121917]
80 [[ 0.10208294  0.2027912 ]] [ 0.29731771]
100 [[ 0.10063919  0.20084961]] [ 0.29918078]
120 [[ 0.10019579  0.20025893]] [ 0.29974979]
140 [[ 0.1000599   0.20007896]] [ 0.2999236]
160 [[ 0.10001831  0.2000241 ]] [ 0.29997668]
180 [[ 0.1000056   0.20000736]] [ 0.29999286]
200 [[ 0.1000017   0.20000222]] [ 0.29999784]

일단 gitHub에서 TensorFlow 소스코드를 복제(clone)해와야 한다.

tensorflow)root@jemin-virtual-machine: git clone https://github.com/tensorflow/tensorflow.git

정상적으로 다운이 완료 되면, 예제를 실행 해보자.

이 예제는 "LeNet -5 - like convolutional MNIST model"이라고 한다.

(tensorflow)$ cd tensorflow/models/image/mnist
(tensorflow)$ python convolutional.py

실행결과

Initialized!
Epoch 0.00
Minibatch loss: 12.054, learning rate: 0.010000
Minibatch error: 90.6%
Validation error: 84.6%
Epoch 0.12
Minibatch loss: 3.285, learning rate: 0.010000
Minibatch error: 6.2%
Validation error: 7.0%
Epoch 0.23
Minibatch loss: 3.473, learning rate: 0.010000
Minibatch error: 10.9%
Validation error: 3.7%
Epoch 0.35
Minibatch loss: 3.221, learning rate: 0.010000
Minibatch error: 4.7%
Validation error: 3.2%
Epoch 0.47
Minibatch loss: 3.193, learning rate: 0.010000
Minibatch error: 4.7%
Validation error: 2.7%
Epoch 0.58
Minibatch loss: 3.301, learning rate: 0.010000
Minibatch error: 9.4%
Validation error: 2.5%
Epoch 0.70
Minibatch loss: 3.203, learning rate: 0.010000
Minibatch error: 6.2%
Validation error: 2.8%
Epoch 0.81
Minibatch loss: 3.022, learning rate: 0.010000
Minibatch error: 4.7%
Validation error: 2.6%
Epoch 0.93
Minibatch loss: 3.131, learning rate: 0.010000
Minibatch error: 6.2%
Validation error: 2.1%
Epoch 1.05
Minibatch loss: 2.954, learning rate: 0.009500
Minibatch error: 3.1%
Validation error: 1.6%
Epoch 1.16
Minibatch loss: 2.854, learning rate: 0.009500
Minibatch error: 0.0%
Validation error: 1.8%
Epoch 1.28
Minibatch loss: 2.825, learning rate: 0.009500
Minibatch error: 1.6%
Validation error: 1.4%
Epoch 1.40
Minibatch loss: 2.938, learning rate: 0.009500
Minibatch error: 7.8%
Validation error: 1.5%
Epoch 1.51
Minibatch loss: 2.767, learning rate: 0.009500
Minibatch error: 0.0%
Validation error: 1.9%
Epoch 1.63
Minibatch loss: 2.771, learning rate: 0.009500
Minibatch error: 3.1%
Validation error: 1.4%
Epoch 1.75
Minibatch loss: 2.844, learning rate: 0.009500
Minibatch error: 4.7%
Validation error: 1.2%
Epoch 1.86
Minibatch loss: 2.694, learning rate: 0.009500
Minibatch error: 0.0%
Validation error: 1.3%
Epoch 1.98
Minibatch loss: 2.650, learning rate: 0.009500
Minibatch error: 0.0%
Validation error: 1.5%
Epoch 2.09
Minibatch loss: 2.667, learning rate: 0.009025
Minibatch error: 1.6%
Validation error: 1.4%
Epoch 2.21
Minibatch loss: 2.658, learning rate: 0.009025
Minibatch error: 1.6%
Validation error: 1.2%
Epoch 2.33
Minibatch loss: 2.640, learning rate: 0.009025
Minibatch error: 3.1%
Validation error: 1.2%
Epoch 2.44
Minibatch loss: 2.579, learning rate: 0.009025
Minibatch error: 1.6%
Validation error: 1.1%
Epoch 2.56
Minibatch loss: 2.568, learning rate: 0.009025
Minibatch error: 0.0%
Validation error: 1.2%
Epoch 2.68
Minibatch loss: 2.554, learning rate: 0.009025
Minibatch error: 1.6%
Validation error: 1.1%
Epoch 2.79
Minibatch loss: 2.503, learning rate: 0.009025
Minibatch error: 0.0%
Validation error: 1.2%
Epoch 2.91
Minibatch loss: 2.487, learning rate: 0.009025
Minibatch error: 0.0%
Validation error: 1.2%
Epoch 3.03
Minibatch loss: 2.463, learning rate: 0.008574
Minibatch error: 1.6%
Validation error: 1.2%
Epoch 3.14
Minibatch loss: 2.458, learning rate: 0.008574
Minibatch error: 1.6%
Validation error: 1.1%
Epoch 3.26
Minibatch loss: 2.410, learning rate: 0.008574
Minibatch error: 0.0%
Validation error: 1.4%
Epoch 3.37
Minibatch loss: 2.496, learning rate: 0.008574
Minibatch error: 3.1%
Validation error: 1.3%
Epoch 3.49
Minibatch loss: 2.399, learning rate: 0.008574
Minibatch error: 1.6%
Validation error: 1.1%
Epoch 3.61
Minibatch loss: 2.377, learning rate: 0.008574
Minibatch error: 0.0%
Validation error: 1.1%
Epoch 3.72
Minibatch loss: 2.333, learning rate: 0.008574
Minibatch error: 0.0%
Validation error: 1.1%
Epoch 3.84
Minibatch loss: 2.312, learning rate: 0.008574
Minibatch error: 0.0%
Validation error: 1.2%
Epoch 3.96
Minibatch loss: 2.300, learning rate: 0.008574
Minibatch error: 1.6%
Validation error: 1.1%
Epoch 4.07
Minibatch loss: 2.276, learning rate: 0.008145
Minibatch error: 0.0%
Validation error: 1.1%
Epoch 4.19
Minibatch loss: 2.250, learning rate: 0.008145
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 4.31
Minibatch loss: 2.233, learning rate: 0.008145
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 4.42
Minibatch loss: 2.217, learning rate: 0.008145
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 4.54
Minibatch loss: 2.324, learning rate: 0.008145
Minibatch error: 3.1%
Validation error: 1.0%
Epoch 4.65
Minibatch loss: 2.212, learning rate: 0.008145
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 4.77
Minibatch loss: 2.174, learning rate: 0.008145
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 4.89
Minibatch loss: 2.211, learning rate: 0.008145
Minibatch error: 1.6%
Validation error: 1.0%
Epoch 5.00
Minibatch loss: 2.193, learning rate: 0.007738
Minibatch error: 1.6%
Validation error: 1.0%
Epoch 5.12
Minibatch loss: 2.148, learning rate: 0.007738
Minibatch error: 3.1%
Validation error: 1.0%
Epoch 5.24
Minibatch loss: 2.153, learning rate: 0.007738
Minibatch error: 3.1%
Validation error: 1.0%
Epoch 5.35
Minibatch loss: 2.111, learning rate: 0.007738
Minibatch error: 1.6%
Validation error: 0.9%
Epoch 5.47
Minibatch loss: 2.084, learning rate: 0.007738
Minibatch error: 1.6%
Validation error: 0.8%
Epoch 5.59
Minibatch loss: 2.054, learning rate: 0.007738
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 5.70
Minibatch loss: 2.043, learning rate: 0.007738
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 5.82
Minibatch loss: 2.134, learning rate: 0.007738
Minibatch error: 3.1%
Validation error: 1.0%
Epoch 5.93
Minibatch loss: 2.006, learning rate: 0.007738
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 6.05
Minibatch loss: 2.048, learning rate: 0.007351
Minibatch error: 3.1%
Validation error: 0.9%
Epoch 6.17
Minibatch loss: 1.988, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 1.1%
Epoch 6.28
Minibatch loss: 1.957, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 6.40
Minibatch loss: 1.971, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 6.52
Minibatch loss: 1.927, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 6.63
Minibatch loss: 1.912, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 6.75
Minibatch loss: 1.901, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 6.87
Minibatch loss: 1.886, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 6.98
Minibatch loss: 1.894, learning rate: 0.007351
Minibatch error: 1.6%
Validation error: 1.0%
Epoch 7.10
Minibatch loss: 1.859, learning rate: 0.006983
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 7.21
Minibatch loss: 1.844, learning rate: 0.006983
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 7.33
Minibatch loss: 1.836, learning rate: 0.006983
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 7.45
Minibatch loss: 1.887, learning rate: 0.006983
Minibatch error: 3.1%
Validation error: 0.9%
Epoch 7.56
Minibatch loss: 1.808, learning rate: 0.006983
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 7.68
Minibatch loss: 1.822, learning rate: 0.006983
Minibatch error: 1.6%
Validation error: 0.9%
Epoch 7.80
Minibatch loss: 1.782, learning rate: 0.006983
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 7.91
Minibatch loss: 1.772, learning rate: 0.006983
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 8.03
Minibatch loss: 1.761, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 8.15
Minibatch loss: 1.773, learning rate: 0.006634
Minibatch error: 1.6%
Validation error: 0.9%
Epoch 8.26
Minibatch loss: 1.742, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 8.38
Minibatch loss: 1.744, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 8.49
Minibatch loss: 1.719, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 8.61
Minibatch loss: 1.700, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 8.73
Minibatch loss: 1.700, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 8.84
Minibatch loss: 1.801, learning rate: 0.006634
Minibatch error: 1.6%
Validation error: 0.8%
Epoch 8.96
Minibatch loss: 1.666, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 9.08
Minibatch loss: 1.666, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 9.19
Minibatch loss: 1.649, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 9.31
Minibatch loss: 1.676, learning rate: 0.006302
Minibatch error: 1.6%
Validation error: 0.8%
Epoch 9.43
Minibatch loss: 1.626, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 9.54
Minibatch loss: 1.621, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 9.66
Minibatch loss: 1.606, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 9.77
Minibatch loss: 1.596, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 9.89
Minibatch loss: 1.602, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.9%
Test error: 0.8%

실행 시간은 아래와 같다.

- MNIST -

real 26m26.382s

user 77m58.714s

sys 19m21.539s

가상화 환경에서 CPU만 가지고 처리해서 오래걸린것 같다.

그래도 i7 4세대 프로세서, ram 24gb, ssd pro 256gb에서 실행한것인데, 대략 Android 5.0 컴파일 시간의 반 정도의 시간이 걸린것 같다.

이론적인 부분은 아래의 이전 포스트를 확인한다.

Confusion Matrix

가장 쉽게 생성할 수 있는 방법은 table()을 이용 하는 것이다.

vec1 <- factor(c("ham","ham","spam","ham"))
vec2 <- factor(c("spam","ham","spam","ham"))
> table(vec1,vec2)
      vec2
vec1   ham spam
  ham    2    1
  spam   0    1

좀더 상세한 결과를 원한다면, gmodels package의 CrossTable() function을 사용하자.

default로써 CrossTable() output으로써 proportion값을 포함하고 있다.

> class(sms_test_pred)
[1] "factor"
> class(sms_raw_test$type)
[1] "factor"
> CrossTable(sms_test_pred, sms_raw_test$type)

 
   Cell Contents
|-------------------------|
|                       N |
| Chi-square contribution |
|           N / Row Total |
|           N / Col Total |
|         N / Table Total |
|-------------------------|

 
Total Observations in Table:  1390 

 
              | sms_raw_test$type 
sms_test_pred |       ham |      spam | Row Total | 
--------------|-----------|-----------|-----------|
          ham |      1202 |        31 |      1233 | 
              |    16.109 |   106.250 |           | 
              |     0.975 |     0.025 |     0.887 | 
              |     0.996 |     0.169 |           | 
              |     0.865 |     0.022 |           | 
--------------|-----------|-----------|-----------|
         spam |         5 |       152 |       157 | 
              |   126.514 |   834.437 |           | 
              |     0.032 |     0.968 |     0.113 | 
              |     0.004 |     0.831 |           | 
              |     0.004 |     0.109 |           | 
--------------|-----------|-----------|-----------|
 Column Total |      1207 |       183 |      1390 | 
              |     0.868 |     0.132 |           | 
--------------|-----------|-----------|-----------|

위 데이터를 가지고 직접 precision, recall, accuracy, error rate 등을 계산 하면 된다.

Beyond accuracy, other measures of performance

자동 계산을 위해서는 또 다른 package인 caret을 사용 한다.

해당 package는 classification과 regression training에 대한 정보를 담고 있으며, Max Kuhn에 의해서 만들어 졌다.

preparing, training, evaluating, and visualizing machine learning models and data 등을 지원 한다.

Journal of statistical software에 "Building predictive models in R using the caret package."을 보면 상세한 내용을 알 수 있다.

아래와 같이 사용하며, 단지 positive label이 무엇인지만 정의해 주면 된다. 

> library(caret)
필요한 패키지를 로딩중입니다: lattice
필요한 패키지를 로딩중입니다: ggplot2
> confusionMatrix(sms_test_pred, sms_raw_test$type, positive = "spam")
Confusion Matrix and Statistics

          Reference
Prediction  ham spam
      ham  1202   31
      spam    5  152
                                          
               Accuracy : 0.9741          
                 95% CI : (0.9643, 0.9818)
    No Information Rate : 0.8683          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.8795          
 Mcnemar's Test P-Value : 3.091e-05       
                                          
            Sensitivity : 0.8306          
            Specificity : 0.9959          
         Pos Pred Value : 0.9682          
         Neg Pred Value : 0.9749          
             Prevalence : 0.1317          
         Detection Rate : 0.1094          
   Detection Prevalence : 0.1129          
      Balanced Accuracy : 0.9132          
                                          
       'Positive' Class : spam    

이제 일반적인 statistics를 알아보자.

■ The kappa statistic

■ Leave One Out Cross - Validation

cross validation 방법의 일종으로써

총 N개의 샘플을 생성 한다.

이때 N개는 해당 데이터 셋의 샘플들의 총 수를 말한다

즉, 100개의 데이터 샘플이 있다면 N은 100이다.

이에대해서 이제 N개의 모델을 생성 한다. 그리고 각각의 모델을 생성할때

하나의 샘플씩을 제외한다. 즉 N개의 모델은 N-1개의 샘플 수로 만들어지는 것이다.

그리고 각각의 모델에 대해서 그 제외한 샘플 1개를 이용해서 test를 수행한다.

그렇게 나온 test 결과 N개에대한 평균 값을 구하게 된다.

MSE가 측정 방법이라고 한다면 아래의 식과 같다.

$MSE_{ i }=(y_{ 2 }-\hat { y_{ 2 } } )^{ 2 }$

구하려고 하는 성능의 척도는 결국 아래와 같다.

$CV_{(n)}=\frac{1}{n}\sum_{i=1}^{n}{MSE_{i}}$

LOOCV의 장점은 결국 모든 샘플에 대해서 다 한번씩은 test하는 것이 되므로

어떠한 랜덤성도 존재하지 않게 된다. 따라서 굉장히 stable한 결과를 얻을 수 있다.

단점으로는 N개의 모델을 만들어서 N번 테스트 하게 되므로 computing time이 굉장히 오래 걸리게 된다.

또한 K-fold Cross validation에 비해서 bias의 비중이 높게 된다. 증 다양성을 포함하기 어렵게 된다.

아래는 작업을 수행한 코드이다.

# LOOCV 처리방법
library(ISLR)
glm.fit=glm(mpg~horsepower, data=Auto) #mpg를 horsepower에 대해서 linear regression을 힙니다.
library(boot) #cv.glm을 사용하기 위한 라이브러리입니다.
cv.err=cv.glm(Auto, glm.fit) #cv.glm함수를 통해 LOOCV를 시행합니다.
cv.err$delta #delta는 cross-validation 결과들을 담고 있습니다.

#linear model, polynomial model에 한해서 LOOCV에 간단한 형태의 공식을 사용할 수 있습니다.
loocv=function(fit){
    h=lm.influence(fit)$h
    mean((residuals(fit)/(1-h))^2)
}

loocv(glm.fit)

cv.error=rep(0,5)
degree=1:5 #1차부터 5차함수까지 fitting해 봅시다.
for(d in degree){
    glm.fit=glm(mpg~poly(horsepower,d), data=Auto) #poly함수를 이용해 d차 함수를 fit할 수 있습니다.
    cv.error[d]=loocv(glm.fit) #LOOCV 결과를 각 벡터에 저장합니다.
}
plot(degree,cv.error,type="b") #LOOCV 결과를 plot합니다.

5차원에서 가장 에러가 적은것을 알 수 있다.

자료출처: https://smlee729.github.io/r/machine%20learning/2015/03/19/1-loocv.html

Bayesian rule에 대해서는 이전 theorem post를 참조 한다.

Bayes algorithm을 이용해서 mobile phone spam을 filtering하는 방법

Short Message Service (SMS)를 구분 하는 작업을 수행 한다.

Step 1 - collecting data

naive bayes classifer를 개발하기 위해서 SMS Spamcollection으로 부터 알고리즘을 적용 한다.

http://www.dt.fee.unicamp.br/~tiago/smsspamcollection/

이 스팸 데이터의 포맷은 아래와 같다.

spam 라벨은 junk message를 의미하고, ham 라벨은 legitimate message를 의미한다.

ham massage의 예:

• Better. Made up for Friday and stuffed myself like a pig yesterday. Now I feel bleh. But at least its not writhing pain kind of bleh.
• If he started searching he will get job in few days. He have great potential and talent.
• I got another job! The one at the hospital doing data analysis or something, starts on monday! Not sure when my thesis will got finished

spam message의 예:

• Congratulations ur awarded 500 of CD vouchers or 125gift guaranteed & Free entry 2 100 wkly draw txt MUSIC to 87066
• December only! Had your mobile 11mths+? You are entitled to update to the latest colour camera mobile for Free! Call The Mobile Update Co FREE on 08002986906
• Valentines Day Special! Win over £1000 in our quiz and take your partner on thetrip of a lifetime! Send GO to 83600 now. 150p/msg rcvd.

#데이터 읽기
sms_raw <- read.csv("sms_spam.csv",stringsAsFactors =  FALSE)

#구조 분석 파악
> str(sms_raw)
'data.frame':	5559 obs. of  2 variables:
 $ type: chr  "ham" "ham" "ham" "spam" ...
 $ text: chr  "Hope you are having a good week. Just checking in" "K..give back my thanks." "Am also doing in cbe only. But have to pay." "complimentary 4 STAR Ibiza Holiday or 짙10,000 cash needs your URGENT collection. 09066364349 NOW from Landline not to lose out"| __truncated__ ...

#데이터 타입 분석
> class(sms_raw)
[1] "data.frame"

#변경을 수행 한다.
> sms_raw$type <- factor(sms_raw$type)

#변경 결과를 확인하면 정상적으로 Factor가 된것을 알 수 있다.
> str(sms_raw)
'data.frame':	5559 obs. of  2 variables:
 $ type: Factor w/ 3 levels "All done, all handed in. Don't know if mega shop in asda counts as celebration but thats what i'm doing!",..: 2 2 2 3 3 2 2 2 3 2 ...
 $ text: chr  "Hope you are having a good week. Just checking in" "K..give back my thanks." "Am also doing in cbe only. But have to pay." "complimentary 4 STAR Ibiza Holiday or 짙10,000 cash needs your URGENT collection. 09066364349 NOW from Landline not to lose out"| __truncated__ ...

이상하게도 기본적으로 제공되는 csv 파일을 가지고 작업을 수행할 경우 문제가 발생한다.
따라서 아래의 코드 작업을 수행 해야 한다.

# 부족한 데이터를 삽입한다.
sms_raw[1072,"type"] <- factor("ham")
sms_raw[1072,"text"] <- c("All done, all handed in. Don't know if mega shop in asda counts as celebration but thats what i'm doing!")

부족한 데이터 삽입전에는 아래와 같이 요인이 3개가 된다. 하나가 요인 값이랑 text 값이 잘못 되어 있기 때문이다.

> table(sms_raw$type)
All done, all handed in. Don't know if mega shop in asda counts as celebration but thats what i'm doing! 
                                                                                                       1 
                                                                                                     ham 
                                                                                                    4811 
                                                                                                    spam 
                                                                                                     747

잘못된 데이터를 수정 하면 아래와 같이 바뀐다.

> table(sms_raw$type)
 ham spam 
4812  747 

■ Data preparation - processing text data for analysis

Document Classification은 주어진 문서를 하나 이상의 분류로 구분하는 문제이다.

지금 다루는 스팸 인지 아닌지를 구분하는 것이 문서 분류의 가장 흔한 예라고 할 수 있다.

또 다른 예로는 제품 리뷰 글을 보고 해당 리뷰가 제품에 대한 긍정적인 리뷰인지 부정적인 리뷰인지를 구분하는 

감성 분석 (Sentiment Analysis)이 있다.

text 처리를 위해서 특화된 package인 Text Mining (tm)의 사용에 대해서 다루 겠다.

TM은 Ingo Feinerer에 의해서 Vienna University에서 박사학위 논문의 일환으로 만들어진 Package 이다.

전공은 Economics and Business이다. 해당 package에 대해서 좀 더 깊이 있는 학습을 원할경우 아래의 사이트를 참조 하자.

http://tm.r-forge.r-project.org/

install.packages("tm")

library(tm)

우선 corpus function을 이용해서 text 데이터를 오브젝트로 생성 한다.

포멧팅 오류를 막기 위해서, corpus 실행전에 UTF-8 en-coding을 수행 하자.

그다음 corpus 함수를 수행 한다.

아래와 같이 corpus 함수는 text를 정의하고 있는 파라메터를 입력 받는다. 이를 위해서 text를 일련의 vector로 만들었었다. 이것을 넘겨줄때 VectorSource()를 이용해서 vector임을 corpus에게 알려주게 된다.

최종적으로 corpus는 sms_corpus라고 이름 지어진 Object를 반환하게 된다. 이 Object를 이용해서 이제 Text Data 처리를 쉽게 할 수 있다.

# try to encode again to UTF-8
sms_raw$text <- iconv(enc2utf8(sms_raw$text),sub="byte")
sms_corpus <- Corpus(VectorSource(sms_raw$text))
# examine the sms corpus
print(sms_corpus)
inspect(sms_corpus[1:3])
Metadata:  corpus specific: 0, document level (indexed): 0
Content:  documents: 3
[[1]]
Metadata:  7
Content:  chars: 49
[[2]]
Metadata:  7
Content:  chars: 23
[[3]]
Metadata:  7
Content:  chars: 43

Corpus는 정말로 유연한 문서 읽기 라이브러리이다. 

PDF나 Microsoft Word 문서 같은 것들을 읽을 수 있다. 좀 더 조사하고 싶다면, Data Input section에 대해서

tm package의 vignette를 이용해서 관련 문서를 확인하다.

command: print(vignette("tm"))

본격적으로 text를 분할해서 단어를 분석하기 전에 기본적으로 

데이터를 클린하는 것이 필요하다.

즉, punctuation과 의미 없는 단어들을 말한다.

각각의 단계는 아래와 같다.

# clean up the corpus using tm_map()
# 소문자로 모두 변경 한다.
corpus_clean <- tm_map(sms_corpus, content_transformer(tolower))
# 숫자를 모두 제거 한다.
corpus_clean <- tm_map(corpus_clean, removeNumbers)
# stop words로 알려진, "to", "and", "but", "or"를 모두 제거 한다.
# 이것을 위해서 stop word를 새로 정의하지 않고, 알려진 리스트인 stopword()를 이용 한다.
corpus_clean <- tm_map(corpus_clean, removeWords, stopwords())
# punctuation을 제거 한다.
corpus_clean <- tm_map(corpus_clean, removePunctuation)
# single space로 모두 만든다. 즉, whitespace(여백)을 모두 제거 한다.
corpus_clean <- tm_map(corpus_clean, stripWhitespace)

각각 변경된 내부 구조를 살펴 보면 아래와 같다.

> lapply(sms_corpus[1:3],as.character)
$`1`
[1] "Hope you are having a good week. Just checking in"
$`2`
[1] "K..give back my thanks."
$`3`
[1] "Am also doing in cbe only. But have to pay."
> lapply(corpus_clean[1:3],as.character)
$`1`
[1] "hope good week just checking "
$`2`
[1] "kgive back thanks"
$`3`
[1] " also cbe pay"

data pre-processing의 마지막 단계는 tokenization 단계이다.

해당 단계에서는 각각의 word를 하나의 컴포넌트로 구분하게 된다.

DocumentTermMatrix() 함수를 이용해서 sparse matrix를 구성 한다.

기본적으로 그냥 생성하면 SMS message를 나타내는 것은 5,559개 이며, term을 나타내는 것은 7,000개 이다.

> sms_dtm <- DocumentTermMatrix(corpus_clean)
> sms_dtm
<>
Non-/sparse entries: 42635/44407129
Sparsity           : 100%
Maximal term length: 40
Weighting          : term frequency (tf)

이제 sparse matrix을 생성 했으므로 word frequency를 수행 할 수 있다.

■ Data preparation - creating training and test datasets 

트레이닝과 테스트로 데이터를 나눠서 나중에 평가를 위해서도 사용 할 수 있도록 한다.

75 vs 25로 분할한다. 그리고 이미 임의로 섞여 있기 때문에 그냥 나누면 된다.

raw data / document-term matrix / corpus 이렇게 세가지 종류의 데이터 셋을 모두 분할 한다.

# creating training and test datasets
sms_raw_train <- sms_raw[1:4169, ]
sms_raw_test  <- sms_raw[4170:5559, ]

sms_dtm_train <- sms_dtm[1:4169, ]
sms_dtm_test  <- sms_dtm[4170:5559, ]

sms_corpus_train <- corpus_clean[1:4169]
sms_corpus_test  <- corpus_clean[4170:5559]

비율이 정상적으로 분할 되었는지는 아래의 명령어를 통해서 간단히 알 수 있다.

> # check that the proportion of spam is similar
> prop.table(table(sms_raw_train$type))
      ham      spam 
0.8647158 0.1352842 
> prop.table(table(sms_raw_test$type))
      ham      spam 
0.8683453 0.1316547 

■ Visualizing text data - word clouds

word cloud는 각 단어의 빈도를 text data로 묘사하게 된다.

여기서 보여지는 단어들의 배치는 랜덤이며 그 글자 크기가 빈도를 나타낸다.

이것을 통해서 경향을 파악할 수 있기 때문에 word cloud는 많이 사용 된다.

spam과 ham을 구분하는 예제에서는 word cloud를 그려 봄으로서 생성하려는 classifier가 성공적으로 동작 할지 안할지를

가늠 할 수 있게 된다.

install.packages("wordcloud")
library(wordcloud)
wordcloud(sms_corpus_train, min.freq = 30, random.order = FALSE)

# subset the training data into spam and ham groups
spam <- subset(sms_raw_train, type == "spam")
ham  <- subset(sms_raw_train, type == "ham")
wordcloud(spam$text, max.words = 40, scale = c(3, 0.5))
wordcloud(ham$text, max.words = 40, scale = c(3, 0.5))

findFreqTerms(sms_dtm_train, 5)

# indicator features for frequent words
sms_dict<- findFreqTerms(sms_dtm_train, 5)
sms_train <- DocumentTermMatrix(sms_corpus_train, list(dictionary = sms_dict))
sms_test  <- DocumentTermMatrix(sms_corpus_test, list(dictionary = sms_dict))

# convert counts to a factor
convert_counts <- function(x) {
  x <- ifelse(x > 0, 1, 0)
  x <- factor(x, levels = c(0, 1), labels = c("No", "Yes"))
}
# apply() convert_counts() to columns of train/test data
sms_train <- apply(sms_train, MARGIN = 2, convert_counts)
sms_test  <- apply(sms_test, MARGIN = 2, convert_counts)

## Step 3: Training a model on the data ----
library(e1071)
sms_classifier <- naiveBayes(sms_train, sms_raw_train$type)
sms_classifier
$tables$yesterday
                  yesterday
sms_raw_train$type          No         Yes
              ham  0.995284327 0.004715673
              spam 0.996453901 0.003546099

## Step 4: Evaluating model performance ----
sms_test_pred <- predict(sms_classifier, sms_test)

> library(gmodels)
> CrossTable(sms_test_pred, sms_raw_test$type,
+            prop.chisq = FALSE, prop.t = FALSE, prop.r = FALSE,
+            dnn = c('predicted', 'actual'))

 
   Cell Contents
|-------------------------|
|                       N |
|           N / Col Total |
|-------------------------|

 
Total Observations in Table:  1390 

 
             | actual 
   predicted |       ham |      spam | Row Total | 
-------------|-----------|-----------|-----------|
         ham |      1202 |        31 |      1233 | 
             |     0.996 |     0.169 |           | 
-------------|-----------|-----------|-----------|
        spam |         5 |       152 |       157 | 
             |     0.004 |     0.831 |           | 
-------------|-----------|-----------|-----------|
Column Total |      1207 |       183 |      1390 | 
             |     0.868 |     0.132 |           | 
-------------|-----------|-----------|-----------|

## Step 5: Improving model performance ----
sms_classifier2 <- naiveBayes(sms_train, sms_raw_train$type, laplace = 1)
sms_test_pred2 <- predict(sms_classifier2, sms_test)
CrossTable(sms_test_pred2, sms_raw_test$type,
           prop.chisq = FALSE, prop.t = FALSE, prop.r = FALSE,
           dnn = c('predicted', 'actual'))
   Cell Contents
|-------------------------|
|                       N |
|           N / Col Total |
|-------------------------|

 
Total Observations in Table:  1390 

 
             | actual 
   predicted |       ham |      spam | Row Total | 
-------------|-----------|-----------|-----------|
         ham |      1203 |        31 |      1234 | 
             |     0.997 |     0.169 |           | 
-------------|-----------|-----------|-----------|
        spam |         4 |       152 |       156 | 
             |     0.003 |     0.831 |           | 
-------------|-----------|-----------|-----------|
Column Total |      1207 |       183 |      1390 | 
             |     0.868 |     0.132 |           | 
-------------|-----------|-----------|-----------|

Rounding of Numbers

정확한 계산을 필요로 하는 경우도 있지만, 때로는 반올림을 할 필요도 있다. 이 경우를 위해서 R은 아래와 같은 

몇가지 반올림 관련 함수들을 제공 한다.

ceiling(x)

floor(x)

trunc(x, ...)

round(x, digits =)

signif(x, digits =)

digits는 반올림할 decimal point의 위치를 결정 한다. 

상수면 소수점(decimal point)를 나타내며, 음수면 10,100의 자리를 나타낸다 (multiples of 10, 100).

round(123.456, digits=2)
[1] 123.46
round(-123.456, digits=-2)
[1] 100

만약 유효숫자(significant digits)를 숫자의 크기에 상관없이 명시하고 싶다면,

signif()를 사용할 수 있다.

> signif(-123.456,digits=4)
[1] -123.5
> signif(-123.456,digits=5)
[1] -123.46
> signif(-123.456,digits=1)
[1] -100
> signif(-123.456,digits=1)
[1] -100

round와 signif는 5를 기준으로 작은지 큰지에 따라 rounding을 수행 한다. 이와 다르게 한쪽 방향으로만 rounding을 수행하는

함수들도 존재 한다.

floor(x) 가장 작은 정수로 x가 결정 된다.

> floor(123.45)
[1] 123
> floor(-123.45)
[1] -124

ceiling(x) 가장 큰 정수로 x가 결정 된다.

> ceiling (123.45)
[1] 124
> ceiling (-123.45)
[1] -123

trunc(x) 가장 가까운 정수로 반올림되는데 방향이 0이랑 가장 근접한 쪽이다.

> trunc(123.45)
[1] 123
> trunc(-123.45)
[1] -123

doBy 패키지에는 summaryBy(), orderBy(), sampleBy()와 같은 특정 값에 따라 데이터를 처리하는 유용한 함수들이 있다.

다음은 이 함수들의 특징을 요약한 것이다.

▣ summaryBy()

저것을 사용하기 위해서 formula 문법에 대해서 알아야 한다.

summaryBy는 특정 조건에 따라 컬럼의 값을 요약해서 보여주게 된다.

아래의 경우 Sepal.Length와 Sepal.With를 Species에 따라 출력해서 보는 방식이다.

> summaryBy(Sepal.Width + Sepal.Length ~ Species, iris)
     Species Sepal.Width.mean Sepal.Length.mean
1     setosa            3.428             5.006
2 versicolor            2.770             5.936
3  virginica            2.974             6.588

▣ orderBy

doBy::orderBy: 포뮬러에 따라 데이터를 정렬한다.

formula, # 정렬할 기준을 지정한 포뮬러

# ~의 좌측은 무시하며, ~ 우측에 나열한 이름에 따라 데이터가 정렬된다.

data, # 정렬할 데이터

기본은 ascending 이기 때문에 descending으로 하고 싶을 경우 - (마이너)를 붙인다.

orderBy(~-Mean, dataFrame) 이런 식이다.

# 다음의 코드는 Species, Sepal.Width 순으로 정렬한 예이다.
> orderBy(~Species + Sepal.Width, iris)
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
42           4.5         2.3          1.3         0.3     setosa
9            4.4         2.9          1.4         0.2     setosa
2            4.9         3.0          1.4         0.2     setosa
13           4.8         3.0          1.4         0.1     setosa
14           4.3         3.0          1.1         0.1     setosa

▣ sampleBy

샘플 함수는 정보를 추출하는 함수이다. 아래와 같다.

10개 중 5개를 랜덤으로 추출한다. replace=TRUE를 설정하면 중복을 허용하는 것이다.

sampleBy과 대응 되므로 같이 알아두면 좋다.

> sample(1:10,5)
[1]  1 10  7  3  8
> sample(1:10,5, replace=TRUE)
[1] 8 8 2 3 6

이러한 샘플링은 주어진 데이터를 Training Data와 Test Data로 분리하는데 유용하게 사용 할 수 있다.

절적한 모델 테스팅을 위해선 균일한 분푸로 Test Data를 생성하는 것이 좋다. 이럴떄 sampleBy가 유용하게 사용 된다.

다음은 Iris 데이터에서 각 Species 별로 10%의 데이터를 추출한 예이다. sampleBy()에서 Species를 지정했으므로 각 종별로 5개씩 데이터가 정확히 샘플로 추출된 것을 볼 수 있다.

> sampleBy(~Species, fra=0.1, data=iris)
              Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
setosa.13              4.8         3.0          1.4         0.1     setosa
setosa.24              5.1         3.3          1.7         0.5     setosa
setosa.27              5.0         3.4          1.6         0.4     setosa
setosa.35              4.9         3.1          1.5         0.2     setosa
setosa.40              5.1         3.4          1.5         0.2     setosa
versicolor.59          6.6         2.9          4.6         1.3 versicolor
versicolor.64          6.1         2.9          4.7         1.4 versicolor
versicolor.77          6.8         2.8          4.8         1.4 versicolor
versicolor.84          6.0         2.7          5.1         1.6 versicolor
versicolor.89          5.6         3.0          4.1         1.3 versicolor
virginica.102          5.8         2.7          5.1         1.9  virginica
virginica.104          6.3         2.9          5.6         1.8  virginica
virginica.121          6.9         3.2          5.7         2.3  virginica
virginica.132          7.9         3.8          6.4         2.0  virginica
virginica.144          6.8         3.2          5.9         2.3  virginica

디버깅 방법은 아래와 같이 크게 2가지로 나눠 진다.

• print(), sprint(), cat()을 사용해서 메시지나 객체의 내용을 출력해보는 방법
• browser()를 사용한 코드 디버깅 방법

print()

print()문으로 중간 중간 데이터를 찍어 보면서 수행 하도록 한다.

paste() 함수와 같이 쓰는면 유용하다.

인자간의 공백을 사용하지 않는 함수는 paste0() 함수 이다.

fibo <- function (n) {
    if (n == 1 || n == 2) {
        print (" base case ")
        return (1)
    }
    print ( paste0 (" fibo (", n - 1, ") + fibo (", n - 2, ")"))
    return ( fibo (n - 1) + fibo (n - 2))
}
> fibo(3)
[1] " fibo (2) + fibo (1)"
[1] " base case "
[1] " base case "
[1] 2

sprintf()

print()와 유사하지만 주어진 인자들을 특정한 규직에 맞게 문자열로 변환해 출력 한다.

sprintf("포멧팅 문자열",인자1,인자2,..인자n)
# %d / %f / %s (C언어의 문법과 유사하다)

cat()

행 바꿈을 수행하지 않고 출력하는 특징이 있다.

줄바꿈 하고 싶으면 "\n"을 수행 하면 된다.

이러한 특징 때문에 cat()을 사용해서 테스트 방식의 프로그레스 상태 메시지를 보여줄 수 있다.

> sum_to_ten <- function () {
+     sum <- 0
+     cat(" Adding ...")
+     for (i in 1:10) {
+         sum <- sum + i
+         cat(i, "...")
+     }
+     cat(" Done !", "\n")
+     return (sum)
+ }
> sum_to_ten()
 Adding ...1 ...2 ...3 ...4 ...5 ...6 ...7 ...8 ...9 ...10 ... Done ! 
[1] 55

Cat을 이용한 또다른 디버깅 코드이다.

#for detecting error, analyze each Corpus

sapply(1:length(vTitle), function(x){

    cat(x,'...')

    sms_corpus_clean <- tm_map(sms_corpus[x], content_transformer(tolower))

    if(x%%10 == 0){

        cat('\n')

    }

})

browser()

해당 함수가 호출 되면 명령어의 수행이 중지되고, 해당 시점부터 디버깅 모드가 시작된다.

원하는 부분에다가 직접 browser() 코드를 삽입 해야한다.

그럼 해당 코드 위치에 도달 했을때

Browse[1]> 이라는 명련 표시줄로 변경이 된다.

명령어는 gdb와 대부부 같다.

n,s,c

종료는 Q 이다.

browser()만의 장점은 다음과 같다.

• 어떠한 editor에서도 사용하다.
• 추가적으로 conditional breakpoint의 기능을 구현할수도 있다.

for (i in 1:1024) {
  start_work()
  if (i == 512)
    browser()
  finish_work()
}

R-Studio를 이용한 디버깅 방법

아래와 같이 편집기에서 마우스로 클릭하면 breakpoint를 설정 할 수 있다.

해당 라인에서 Shift + F9를 해도 된다.

breakpoint를 설정하면 아래와 같이 걸려있는 함수를 환경변수 패널에서 확인 할 수 있다.

만약 아래와 같은 경고가 나온다면 이유는 두가지 이다.

첫 번째는 source()를 통해서 해당 파일을 읽어 들이지 않아서 그래서 object가 생성되지 않았을 때이다.

두 번째는 editor 상의 file과 source()로 읽어들인 file간의 차이가 존재 했을 때이다.

위 두상황 모두 source()를 해주면 해결 된다.

정상적으로 breakpoint를 걸었다면,

이제 해당 코드를 실행하면 원하는 위치에서 멈추는거을 알 수 있다.

R Studio는 디버깅 명령어들을 모두 위즈윅 모드로 제공하므로 마우스를 써서 디버깅 하면 된다.

디버깅 실행 화면

만약 source 코드가 없는 상태에서 디버깅 하고 싶다면?

install.packages("devtools")

참고문헌: Debugging with RStudio

R Studio의 특징

대화식 R

R 그래픽 작업

코드 구조화, 다수의 프로젝트 관리

분석 내용 재현

R에 설치된 패키지 관리

분석 보고서 작성 및 공유

본인의 코드 공유 및 다른 사용자와의 협업

다양한 적용 가능 플랫폼

윈도우

리눅스

맥 OS X

웹

Comprehensive R Archive Network

R언어를 일단 잘알아야 RStudio를 통해서 통계분석이 가능하다.

위즈윅(WYSIWYG)에 익숙한 사람들, 즉 마우스로 클릭하고 데이터를 이동하는것에 익숙 했던 사람들은 여간 사용하기 어려운것이 아니다.

하지만, 분석의 재현을 위해서는 반드시 각 분석 단계를 소스코드 형태로 저장해야한다.

그래야 분석 과정을 재현하기 쉽다.

프로그래밍 방식은 또한 

단순히 사전식으로 해당 기능에대해서만 동작하는 위즈윅 방식의 툴보다 훨씬 유연하며 강력한 기능을 조합을통해서 달성 할 수 있다.

항상 최신 버전을 사용하라

View(news()) 를통해서 최신 정보를 알 수 있다.

update.packages()를 통해서 페키지를 업데이트 하라

설치방법

1. R을 설치해야한다.

Rstudio는 기본적으로 2.11 이상의 R을 설치해야한다. 

공식 사이트: http://www.r-project.org/

이곳에서 가장 가까운 서버인

http://cran.nexr.com/

위 사이트에서 Window용 R을 다운 받았다.

Tip: 항상 최신 버전의 알을 수시로 확인해서 다운받아서 설치하자

1.1 업데이트 방법

기본적으로 자동 업데이트는 제공하고 있지 않으므로, 이전 버전을 제거한다음 새로 설치하면 된다.

이전 package 설치본들을 복사한 다음 붙여 넣기 한다.

그 다음 현재 package들을 업데이트 한다.

update.packages(checkBuilt=TRUE, ask=FALSE) 

개인 라이브러리의 경우 아래의 경로에 저장 되어 있다.

현재 R에 설치된 라이브러리 리스트를 얻어서 그대로 재설치하는 방법도 있다.

library()

현재 R의 버전을 알아내는 방법

아래의 명령어를 통해서 설치된 package와 설치된 version을 알 수 있다.

installed.packages()[,c("Package","Version")]
> installed.packages()[,c("Package","Version")]
             Package        Version   
BH           "BH"           "1.58.0-1"
bitops       "bitops"       "1.0-6"   
brew         "brew"         "1.0-6"   
caTools      "caTools"      "1.17.1"  
colorspace   "colorspace"   "1.2-6"   
curl         "curl"         "0.9.3"   

2. R studio 설치

공식 사이트: http://www.rstudio.com/products/rstudio/download/

Tip: RStudio에서

Help > Check for Updates를 클릭하여 신규 업데이트 정보를 확인하자.

업데이트가 필요하면 위 공식 사이트에가서 신규버전으로 다시 설치하자.

아직까지는 자동 업데이트 기능은 없는것 같다.

3. 필수 package 설치 및 로드

R studio의 하단에 보면 package tab이 있다. 이곳에서 이미 설치된 package들의 리스트를 확인 할 수 있다. 

update 버튼을 통해서 기존 package들을 업데이트 할 수도 있다. 

그리고 install을 통해서 새로운 package들을 쉽게 설치 할 수 있다. 

추가로 설치한 리스트

설치된 package를 로드하기위해서는

library 명령어를 타이핑 해도되고

package 창에서 마우스로 check 해도 된다.

단축키 (shortcuts)

Keyboard Shortcuts: R support

GIt을 이용한 버전 관리

git을 설치한다.

자신의 운영체제에 맞춰서 설치파일을 다운로드 한다.

https://git-scm.com/downloads

Next > Next 만을 선택해서 설치한다.

설치확인 방법은 아래와 같이 적당한 디렉터리 아무곳에서나 오른쪽 마우스 클릭을해서

Git bash here 와 같은 명령어들이 모이면 정상적으로 설치된 것이다.

RStudio와 연동하기

메뉴 바에서

Toos -> Global Options -> Git/SVN 에서 Git을 설치한 디렉터리에서 git 실행파일을 연결해 준다.

이렇게 하고 재부팅을 해준다.

아래와 같이 정상적으로 git tab이 생성되는것을 알 수 있다.

자료형

각각의 자료형들간의 변환은 아래의 함수들을 이용 한다.

R에서는 데이터 분석을 위한 유용한 데이터 구조를 기본으로 제공해 준다.

Factor (범주형 자료)

생성 방법

factor(

x, #팩터로 표현하고자 하는 값(주로 문자열 벡터로 지정)

levels, #값의 레벨

ordered #TRUE면 순서형, FALSE면 명목형 데이터를 뜻한다. 기본값은 FALSE다.

)

팩터에서 레벨의 개수를 반환함.

nlevels(

x # 팩터 값

)

중간에 값 넣기

factor(append(as.character(proximity),"TRUE"))

factor level rename

http://www.cookbook-r.com/Manipulating_data/Renaming_levels_of_a_factor/

Vectors

백터는 동질적이다.

같은 자료형을 가진다.

같은 모드를 가지고 있어야 한다.

스칼라와 백터의 관계는 통상 엄밀히 구분하지만 R에서는 딱히 구분하지 않는다.

하나의 값으로 구성된 백터가 스칼라고

여러개 있으면 백터라고 생각하면 된다.

List

이질적이다.

다양한 자료형을 포함 할수 있다.

생성방법

list <- list(x,y,z) #이질적인 데이터도 가능함.

list <- list(mid=0.5, right=0.841, far=0.9)

접근방법

리스트 변수[[인덱스] # 값을 접근함

리스트 변수[인덱스] # (키,값) 의 형태로 서브 리스트를 반환 한다.

제거방법

x[-2] # without 2nd element

x[-c(2,3)] # without 2nd and 3rd

x["arrival_time"] <- NULL

빈리스트를 생성하고 인덱스를 이름으로 지정하는 방법

vecNames <- list.files("./ESLABCollector/")
filenames <- sapply(vecNames,list)  

Foreach를 사용할 때 list 인덱스의 이름을 유지하는 방법

jsonDfList <- foreach (x = rawDataList, .final = function(x) setNames(x, names(filenames))) %do% {

    jsonTodf(x)

}

.final 을 이용한다.

결국 setNames란 아래의 것을 축약한 기능이다.

setNames( 1:3, c("foo", "bar", "baz") )

# this is just a short form of

tmp <- 1:3

names(tmp) <-  c("foo", "bar", "baz")

tmp

List 중간에 값을 삽입 하기

append(testList, list(x=42), 3)

$`1`

[1] 1

$`2`

[1] 2

$`3`

[1] 3

$x

[1] 42

$`4`

[1] 4

append(lst, list(f=1:5), after=0)# after - position to append

Helper 함수를 생성하는 방법

lappend <- function (lst, ...){
lst <- c(lst, list(...))
  return(lst)
}

> a <- list()
> a
list()

> lappend(a,c(1,2,3))
[[1]]
[1] 1 2 3

> lappend(a, c(4,5,6), c(7,8,9))
[[1]]
[1] 4 5 6

[[2]]
[1] 7 8 9

Matrix

백터에 차원을 정해주면 바로 Matrix로 변경 된다.

백터의 dim 속성이 처음엔 NULL 이지만 이것을 주게 되면 matrix가 된다.

dim(2,3) # 2 by 3 행렬

list도 dim 변경을해서 행렬화 할 수 있다.

하지만 list의 특성상 이질적인 행렬을 만들어 낼 수 있으므로 그러한 것들은 조심해야 한다.

Array (배열)

3차원 matrix 부터는 배열이라는 이름을 쓴다.

Factor(요인)

백터처럼 생겼지만 특별한 속성이 있다.

R은 백터에 있는 고유한 값(unique value)의 정보를 얻어 내는데, 이 고유 값들을 요인의 '수준(level)'이라고 일컫는다.

요인들은 간단하고 효율적인 형태로 데이터 프레임에 저장된다.

다른 프로그래밍 언어에서는 요인을 '열거형 값(enumerated value)들로 이루어진 백터로 표현한다.

요인의 사용처는 아래와 같다.

범주형 변수>

요인 하나는 범주형 변수 하나를 나타낼 수 있다. 

범주형 변수들은 분할표, 선형 회귀, 변량 분석, 로지스틱 회귀 분석 등 많은 분야에서 사용 된다.

집단 분류

이것은 데이터 항목에다가 집다에 따른 라벨을 붙이거나 태깅을 할 때 쓰는 방법이다.

데이터 프레임 (Data Frame)

가장 중요한 데이터 타입이다.

표 형태의 데이터를 표현하기 위함이다.

특징

행은 모두 같은 길이이다. 열은 하나의 엘리먼트를 나타낸다.

행렬과 다른점은 각각의 행들이 서로다른 데이터 타입들을 가질 수 있다는 것이다.

즉, 각각의 열들은 같은 타입으로 이뤄지지만, 모든 열이 모두 같을 필요는 없다.

데이터 프레임의 속성은 row.names()를 통해서 알 수 있다.

통상 read.table() 또는 read.csv()를 통해서 생성 한다.

행렬을 data frame으로 변경 할 수 있다. data.matrix() 호출을 통해서

▣ 데이터 프레임 생성하기

data.frame()을 사용

g <- c("A","B")

x <- 1:3

dat <- data.frame(g,x)

▣ 리스트를 프레임 데이터로 변경

frame.data(list)

▣ 프레임의 구조를 알려준다.

str(data)

> str(ac)

'data.frame': 356 obs. of  2 variables:

 $ Time : num  0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 ...

 $ Power: num  14.8 14.1 19.9 14.6 14.1 ...

▣ 프레임에서 열들만 따로 출력하는 방법이다.

data$time

▣ 데이터 프레임의 Index를 알아내는 방법

> a <- data.frame(1:6,11:16)

> a

  X1.6 X11.16

1    1     11

2    2     12

3    3     13

4    4     14

5    5     15

6    6     16

> grep(15,a[,2])

[1] 5

▣ 해당 컬럼 이름이 몇번째 인지를 알아내는 방법

> grep("X11.16",colnames(a))

[1] 2

▣ 프레임에서 행과 열만 따로 출력하는 방법이다.

data$time[0]

▣ 위치로 데이터 프레임의 열 선택하기

스타일에 따라 열 또는 데이터프레임 두가지 자료형중 하나로 반환 된다.

데이터 프레임 스타일의 선택 방법

행렬 방식의 첨자 지정을 이용해서 선택 하는 방법

▣ 프레임 열의 이름을 변경하는 방법이다.

▣ 특정 열을 찾아서 이름을 변경하는 방법

> names(ac)[names(ac)=="Time"] <- c("newTime")

> names(ac)

[1] "newTime" "Power"  

> names(ac)[names(ac)=="newTime"] <- c("Time")

> names(ac)

[1] "Time"  "Power"

▣ 데이터 프레임에 열 추가하기

새로운 열의 이름을 설정하고 "NA"로 채운다.

data$newcol <- NA

새로운 열의 이름을 설정하고 백터를 삽입한다.

data$newcol <- vec

▣ 데이터 프레임에 열을 제거하기

해당열에 NULL을 대입한다.

data$badcol <- NULL

subset(data, select = -badcol)

subset(data, select = c(-badcol, -othercol)

▣ 데이터 프레임 열 순서 변경하기

#숫자로 해당열을 선택해서 변경한다.

dat <- dat[c(1,3,2)]

#이름으로 선택해서 변경 한다.

data <- dat[c("col1","col3","col2")]

▣ 데이터 프레임의 부분 집합 취하기

subset() 함수를 이용해서 부분 집합을 취할 수 있다.

일련의 조건을 기입하고 그 조건에 맞는 행등을 추출 할 수 있다.

subset(climate, Source == "Berkeley", select = c(Year, Anomaly10y))

# 데이터 셋 climate에서

# Source 열이 "Berkeley인 것에 대해서

# Year 와 Anomaly10y의 열에 대해서만 그 값을 추출 

논리 연산자를 조합하면 일정 범위 내의 데이터를 추출 할 수도 있다.

1900년에서 2000년 사이의 source가 "Berkeley"인 행들만 선택 한다.

subset(climate, Source == "Berkeley" & Year >= 1900 & Year <= 2000, select = c(Year, Anomaly10y))

▣ 여러가지 방법을 이용한 데이터 프레임 합치기

데이터 프레임 합치기

cbind # 두개의 데이터를 열로 쌓는다.

rbind # 두개의 데이터를 행으로 쌓는다.

연속적인 리스트라면,

do.call(rbind,mylist)를 하면된다.

주의: 위방법들은 모두 행열이 같은 크기일때만 가능함

크기가 서로 다르면 반복해서 아랫부분을 채워준다고 하지만 이상하게 동작하지 않는다.

서로다른 데이터 프레임 합치기

merge # SQL에서 join에 해당하는 명령어이다

default로 동작은 그냥 다순히 합치는 것이다.

일치하지 않는것이 있다면 데이터의양이 NxM의 양으로 늘어 난다.

일치하는 부분만 합치고 싶으면,

merge(A,B, by="열 이름") 

으로 설정한다.

merge의 다른 option들

## S3 method for class 'data.frame'
merge (x, y, by = intersect(names(x), names(y)),
      by.x = by, by.y = by, all = FALSE, all.x = all, all.y = all,
      sort = TRUE, suffixes = c(".x",".y"),
      incomparables = NULL, ...)

서로다른 크기의 데이터 프레임 합치기

일치하지 않는 컬럼을 단순히 <NA>로 채우고 싶다면,

plyr package의 rbind.fill을 사용하면 된다.

향상된 데이터 프레임 합치기

reshape2 package

본인 정리

Reshape and aggregate data with the R package reshape2

plyr package

본인정리

데이터 프레임 행과 열 크기 알아내기

nrow()

ncol()

help(package = "reshape2")

원본 참고 사이트: http://seananderson.ca/2013/10/19/reshape.html

reshape 2 website: http://had.co.nz/reshape/

서로 다른 데이터 타입들간의 데이터 상호 변환

R에서는 주로 벡터 기반으로 데이터를 처리 한다.

벡터기반 처리는 개별 요소를 for 루프 등으로 하나씩 처리하는 방식보다 빠르게 수행될 뿐만 아니라 손쉽게 병렬화가 가능하다.

따라서 대부분의 Machine Learning이나 통계 함수들은 입력값으로 벡터를 받는다.

따라서 데이터를 벡터로 변환하고 다시 벡터데이터를 가독성 높은 데이터 프레임으로 변경하는 것을 자유 자재로 할 수 있어야한다.

앞으로 보게될 많은 R 코드들은 모두 이러한 벡터 연산들을 이용하기 때문에 벡터를 다루는 함수들을 잘 알고 있어야 한다.

벡터화된 연산과 apply함수를 이용해서 일괄 실행함으로써 계산을 간소화하고 고속화를 도모 한다.

다양한 데이터를 임의의 함수에 적용하 결과를 얻는 방법으로 apply  함수들을 이용하는 방법이 있다.

이들 함수는 벡터, 행렬, 리스트, 데이터 프레임에 적용할 함수를 명시하는 형태로서, 함수형 언어 스타일에 가깝다고 볼 수 있다.

apply 함수들의 요약

▣ Apply 사용법

apply: 배열 또는 행렬에 함수 f를 d방향으로 적용하여 결과를 벡터, 배열 또는 리스트로 반환 한다.

apply( X, d, f)

X #배열 또는 행렬

d # 함수를 적용하는 방향, 1은 행방향, 2는 열 방향

  # c(1,2)는 행과 열 방향 모두를 의미

FUN # 적용할 함수

예제

#합 구하는 함수
> sum(1:10)
[1] 55
# 메트릭스 생성
> d <- matrix(1:9, ncol=3)
> d
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
# 각각의 행을 더해서 벡터로 반환 한다.
> apply(d,1,sum)
[1] 12 15 18
# 각각의 열을 더해서 벡터로 반환 하다.
> apply(d,2,sum) 
[1]  6 15 24

> head(iris[c(1,2,3,4)])
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2
3          4.7         3.2          1.3         0.2
4          4.6         3.1          1.5         0.2
5          5.0         3.6          1.4         0.2
6          5.4         3.9          1.7         0.4
# 벡터 스타일의 선택
> apply(iris[c(1,2,3,4)],2,sum)
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
       876.5        458.6        563.7        179.9 
# 메트릭스 스타일의 선택
> apply(iris[,1:4],2,sum)
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
       876.5        458.6        563.7        179.9

▣ lapply 사용법

리스트로 반환하는 특징이 있는 apply 계열 함수이다.

lapply (X, FUN, ...#추가인자) 

반환 값은 X와 같은 길이의 리스트이다.

X는 벡터, 리스트, 표현식, 데이터 프레임을 받을 수 있다.

리스트를 다른 타입으로 변환 하는 방법 

unlist( x, recursive=FALSE, use.name=TRUE )

반환 값은 벡터이다.

# 벡터 입력을 받아서 리스트로 반환하는 코드이다.
> result <- lapply(1:3, function(x){x*2})
> result
[[1]]
[1] 2

[[2]]
[1] 4

[[3]]
[1] 6

> result[[1]]
> unlist(result) # list 구조를  벡터로 변경 한다.
[1] 2 4 6

입력을 list로 받을 수도 있다.

> x <- list(a=1:3, b=4:6)
> x
$a
[1] 1 2 3

$b
[1] 4 5 6

> lapply(x,mean)
$a
[1] 2

$b
[1] 5

# iris의 데이터 구조를 보여줌
> head(iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

# list로 반환됨 각각의 열들의 평균을 구한다음
> lapply(iris[,1:4],mean)
$Sepal.Length
[1] 5.843333

$Sepal.Width
[1] 3.057333

$Petal.Length
[1] 3.758

$Petal.Width
[1] 1.199333

#각 열에 대한 평균을 구하게 된다.
> colMeans(iris[,1:4])
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
    5.843333     3.057333     3.758000     1.199333

데이터 프레임을 처리한 결과를 리스트로 얻은 뒤 해당 리스트를 다시 데이터 프레임으로 변활할 필요가 있다. 

이 변환은 몇 단계를 거쳐서 처리해야 한다.

• unlist()를 통해 리스트를 벡터로 변환한다.
• matrix()를 사용해 벡터를 행렬로 변환한다.
• as.data.frame()을 사용해 행렬을 데이터 프레임으로 변환한다.
• names()를 사용해 리스트로부터 변수명을 얻어와 데이터 프레임의 각 컬럼에 이름을 부여한다.

iris 데이터의 각 열의 평균을 구한다음

이 결과는 list로 반환 되므로 이것을 다시 데이터 프레임으로 변경해서 출력하는 예제이다.

># 데이터 프레임으로 변겨하는 코드이다.
# matrix()로 먼저 변환 하는 것은 vector를 바로 data frame으로 변경하면 단일 열로 밖에 변경을 못하기 때문이다.
 d <- as.data.frame(matrix(unlist(lapply(iris[,1:4],mean)),ncol=4,byrow=TRUE))
> d
        V1       V2    V3       V4
1 5.843333 3.057333 3.758 1.199333
> d <- as.data.frame(matrix(unlist(lapply(iris[,1:4],mean)),ncol=4,byrow=TRUE))
> d
        V1       V2    V3       V4
1 5.843333 3.057333 3.758 1.199333
#위와 같이 이름이 자동으로 설정되므로 가독성을 위해서 이름을 변경해 준다.
> names(d) <- names(iris[,1:4])
> d
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1     5.843333    3.057333        3.758    1.199333

unlist 방식의 경우 벡터로 바꾸기 때문에 데이터 타입이 서로 다를 경우 엉뚱한 값으로 바꿔버리게 된다.

> x <- list(data.frame(name="foo",value=1),data.frame(name="bar",value=2))
> x
[[1]]
  name value
1  foo     1

[[2]]
  name value
1  bar     2

> unlist(x)
 name value  name value 
    1     1     1     2 

이렇게 서로 다른 데이터로 구성된 list를 데이터 프레임으로 변경 해야한다면,

do.call() 함수를 이용 해야 한다.

> do.call(rbind,x)
  name value
1  foo     1
2  bar     2

하지만, 해당 방식은 속도면에서 문제가 발생 한다.

그래서 rbindlist()를 사용하게 된다.

▣ sapply 사용법

> sapply(iris[,1:4],mean)
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
    5.843333     3.057333     3.758000     1.199333 
> class(sapply(iris[,1:4],mean))
[1] "numeric"

sapply로 연산하고 벡터를 data frame으로 변경하는 방법이다.

> sapply(iris[,1:4],mean)
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
    5.843333     3.057333     3.758000     1.199333 
> class(sapply(iris[,1:4],mean))
[1] "numeric"
> x <- sapply(iris[,1:4],mean)
> as.data.frame(x)
                    x
Sepal.Length 5.843333
Sepal.Width  3.057333
Petal.Length 3.758000
Petal.Width  1.199333
> as.data.frame(t(x))
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1     5.843333    3.057333        3.758    1.199333

sapply를 이용해서 각 컬럼의 데이터 타입을 구할 수도 있다.

> sapply(iris,class)
Sepal.Length  Sepal.Width Petal.Length  Petal.Width      Species 
   "numeric"    "numeric"    "numeric"    "numeric"     "factor"

다음은 조건문을 이용해서 iris 데이터 집합중에서 3보다 큰지 여부를 판단하는 코드이다.

반환값은 행렬이 된다.

> y <- sapply(iris[,1:4], function(x) {x>3})
> class(y)
[1] "matrix"
> head(y)
     Sepal.Length Sepal.Width Petal.Length Petal.Width
[1,]         TRUE        TRUE        FALSE       FALSE
[2,]         TRUE       FALSE        FALSE       FALSE
[3,]         TRUE        TRUE        FALSE       FALSE
[4,]         TRUE        TRUE        FALSE       FALSE
[5,]         TRUE        TRUE        FALSE       FALSE
[6,]         TRUE        TRUE        FALSE       FALSE

주의해야할 사항은 벡터를 반환 하기 떄문에 sapply를 사용할 때는 데이터 타입을 혼합하여 처리하면 안된다.

데이터 타입을 혼합해야 할 경우에는 lapply를 사용하거나 plyr 페키지에 있는 ddply를 사용해야 한다.

▣ tapply() 사용법

tapply{

X 벡터

Index 데이터를 그룹으로 묶을 색인, 팩터를 지정해야 하며 팩터가 아닌 타입이 지정되면 팩터로 형 변환된다.

FUN 각 그룹마다 적용 할 함수

...., 추가 인자. 이 인자들은 FUN에 전달된다.

}

반환 값은 배열이다.

아래의 코드는 간단한 하나의 그룹 1에 대한 합을 나타낸다.

# rep는 1을 10번 반복하라는 뜻이다. 이렇게해서 1:10 까지의 벡터에 factor 데이터 1을 모두 맵핑해서 하나의 그룹으로 만들어 준다.
> rep(1,10)
 [1] 1 1 1 1 1 1 1 1 1 1
> tapply(1:10, rep(1,10),sum)
 1 
55

# 홀수 짝수를 나누어서 연산하는 것이다.
tapply(1:10,1:10%%2 ==1,sum)
FALSE  TRUE 
   30    25 

조금더 복잡한 연산으로 iris 데이터들중에서 length의 평균을 구하는데 그것을 specie에 따라서 구하는 것이다.

> head(iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa
> tapply(iris$Sepal.Length, iris$Species, mean)
    setosa versicolor  virginica 
     5.006      5.936      6.588 

이제 더 복잡한 예제를 살펴보자. 아래의 예제는 추후에 클러스터링 알고리즘을 적용한 다음에 (k-mean clustering)

같은 cluster에 속한 데이터들의 x와 y 좌표의 각가그이 평균을 구하는 용도로 사용 될 수 있다.

#데이터를 아래와 같이 생성 한다.
> m <- matrix(1:8, ncol=2, dimnames=list(c("spring","summer","fall","winter"),c("male","female")))
> m
       male female
spring    1      5
summer    2      6
fall      3      7
winter    4      8

위와 같이 데이터 구조를 생성 했을 때

이때 분기별 과 성별 셀의 합을 구해보자.

상반기 봄, 여름에 대한 남성의 셀 합은 1+2로 구할 수 있다. 그리고 여성의 합은 5+6이다.

하반기 가을, 겨울의 남성의 셀 합은 3+4이고 여성의 합은 7+8이다.

Index를 설정하기 위해서 (n,m)에서 n을 먼저 나열한 뒤 m 값을 나열한다. 

즉, 그룹 (n1, m1), (n2, m2)는 list( c(n1,n2), c(m1,m2))로 표현된다.

> tapply(m, list(c(1,1,2,2,1,1,2,2),c(1,1,1,1,2,2,2,2)),sum)
  1  2
1 3 11
2 7 15

▣ mapply() 사용법

sapply와 유사하지만 다수의 인자를 함수에 넘긴다는 점에서 차이가 존재 한다.

sapply 처럼 행렬과 벡터를 반환 하게 된다.

1) 여러개의 인자를 보내야 하는 경우에 사용 한다.

2) 어떤 함수에 여러개의 벡터를 입력으로 주고 싶은 경우

a,b를 받아서 처리하는 함수가 있을때

c(....), c(....) 이런식으로 2개의 벡터를 입력으로 보내서 연속적인 처리를 할 수 있다.

mapply: 함수에 리스트 또는 벡터로 주어진 인자를 적용한 결과를 반환 한다.

Fun #실행할 함수

... # 적용할 인자

rnorm의 경우 n개의 데이터를 mean에 대해서 sd인 정규 분포를 생성 한다.

그냥 rnorm을 세번 호출하면 되지만 편의를 위해서 mapply를 이용했다.

세개인자를 rnorm함수로 전달해서 연살을 수행하게 된다.

> mapply(rnorm, c(1,2,3),c(0,10,100),c(1,1,1))
[[1]]
[1] 0.5746612

[[2]]
[1] 10.026820  9.574186

[[3]]
[1]  99.53101  98.92185 101.26264

#iris 데이터의 경우 아래와 같이 수행 된다.
> mapply(mean,iris[,1:4])
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
    5.843333     3.057333     3.758000     1.199333 

> sapply(iris[,1:4],mean)
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
    5.843333     3.057333     3.758000     1.199333 

> lapply(iris[,1:4],mean)
$Sepal.Length
[1] 5.843333

$Sepal.Width
[1] 3.057333

$Petal.Length
[1] 3.758

$Petal.Width
[1] 1.199333