Probability and Statistics – Understanding Exp Data

(** 인터넷 익스플로러에서 보면 파이썬 코드 줄이 제대로 적용되어 보이지 않을 수 있습니다.**)

실험에서 얻은 데이터를 파이썬 프로그램을 작성하여 분석할 수 있다. 스프링 운동 실험에서 스프링에 알려진 무게의 추를 달고 늘어난 길이를 측정한 데이터를 가정하자.

Distnace (m) Mass (kg)
0.0865 0.1
0.1015 0.15
0.11 0.203873598
0.125 0.244648318
0.19 0.30479103
0.27 0.326197757
0.28 0.387359837
0.24 0.448521916
0.345 0.508664628
0.321 0.540265036
0.371 0.612640163
0.417 0.673802243
0.451 0.713659531
0.435 0.739041794
0.435 0.775739042
0.442 0.838939857
0.4416 0.9
0.4304 0.95
0.437 1.0

첫 번째 줄에 헤더가 있고, 그 다음 줄부터 스프링 실험에서 측정한 거리와 질량에 대한 숫자가 공백으로 구분되어 있다. 예를 들어, 헤더 다음의 첫 번째 줄의 숫자 0.0865와 0.1은 0.1kg 무게의 추를 달면 스프링이 0.0865m 늘어나는 실험 결과를 보여준다.

스프링에 관한 물리 공식으로 훅의 법칙 F = -k x이 있다. k는 스프링 특성 상수, x는 늘어난 길이, F는 힘이다. 마이너스 기호는 힘의 방향이 스프링이 늘어나는 방향과 반대라는 의미이다.

사용한 스프링의 특성 상수 k를 알고자 위 실험을 수행하였다. 실험 데이터에 있는 추의 무게와 이 추를 달았을 때 늘어난 스프링의 길이를 알면 스프링 특성 상수를 훅의 법칙으로 계산할 수 있다. 첫 번째 실험 데이터 수치로 k를 다음과 같이 계산할 수 있다.

k = - F / x 
  = - M g / x                   (g: 중력 가속도, M: 추의 무게, x: 스프링이 늘어난 길이)
  = - 1kg x 9.81m/s^2 / 0.1m 
  = -98.1 (N/m)                 (N: 뉴튼, kg x m/s^2)

스프링 실험 데이터는 무게와 길이가 있는데 이를 힘과 길이로 약간 가공한 다음 x축은 힘, y축은 길이로 두어 각 데이터의 점을 찍어 그래프를 그려보자.

import pylab

# 파일 SpringData.txt에 스프링 실험 데이터를 저장하고
# 이 파일이 저정된 경로에 맞추어 다음 문자열을 작성한다.
path = 'your path to SpringData.txt'  

def getData(fileName):
    dataFile = open(fileName, 'r')
    distances = []
    masses = []
    discardHeader = dataFile.readline()
    for line in dataFile:
        d, m = line.split(' ')
        distances.append(float(d))
        masses.append(float(m))
    dataFile.close()
    return (masses, distances)

def plotData(inputFile):
    masses, distances = getData(inputFile)
    masses = pylab.array(masses)
    distances = pylab.array(distances)
    forces = masses * 9.81 ## Vector operation
    pylab.plot(forces, distances, 'bo',
               label = 'Measured displacements')
    pylab.title('Measured Displacement of Spring')
    pylab.xlabel('|Force| (Newtons)')
    pylab.ylabel('|Distance (meters)')

위 파이썬 프로그램을 실행해서 그래프를 살펴보자.

>>> plotData(path)
>>> pylab.show()

참고로, 이 파이썬 프로그램에서 파일 열고 닫기와 파일에서 텍스트 읽기를 사용하였다. 텍스트 파일에서 한 줄씩 처리할 때 for문으로 간결하게 처리할 수 있다. 문자열의 split 함수는 인자로 주어진 구분자(delimeter)인 공백 ' '으로 여러 주어진 문자열(line)을 여러 문자열들(d, m)로 나눈다.

dataFile = open(fileName, 'r')
...
for line in dataFile:    // 텍스트 파일의 한 줄을 읽어 for문 몸체를 실행
    d, m = line.split(' ')
...
dataFile.close()

>>> dataFile = open(datafile, 'r')
>>> print(dataFile)
<_io.TextIOWrapper name='C:\\Work\\pythonstudy\\SpringData.txt' mode='r' encoding='cp949'>
>>> s = dataFile.readline()
>>> print(s)
Distnace (m) Mass (kg)

>>> s = dataFile.readline()
>>> print(s)
0.0865 0.1

>>> x,y = s.split(' ')
>>> print(x)
0.0865
>>> print(y)
0.1

>>> x,y = "abc,123".split(',')
>>> print(x)
abc
>>> print(y)
123

파이썬 공부에서 다시 확률 통계 공부로 돌아와서 스프링 실험에 대해서 생각해보자. 스프링 실험을 한 목적은 스프링 특성 상수를 측정하는 것이다. 앞에서 그린 그래프가 y = a x 방정식 형태를 갖춘다면 k를 쉽게 알아낼 수 있다. 따라서 이제부터 할 일은 실험 데이터 그래프와 일치하는(fit) 방정식을 찾는 것이다. 스프링의 경우 F = -k x 식을 알고 있기 때문에 1차 방정식을 찾으면 된다. 하지만 만일 이러한 식을 사전에 알지 못한다면 실험 데이터 그래프가

ax + b = 0
ax^2 + bx + c = 0
ax^3 + bx^2 + cx + d = 0
...

과 같이 1차, 2차, 3차 그리고 그 이상의 차수의 방정식과 일치할 수도 있으니 차수를 다르게 하며 일치하는 방정식을 찾아볼 필요가 있다.

파이썬 함수 pylab.polyfit는 데이터와 차수를 인자로 주면 가장 일치하는 방정식을 계산한다.

### fitData 버전 1:
def fitData(inputFile):
    masses, distances = getData(inputFile)

    masses = pylab.array(masses)
    distances = pylab.array(distances)
    
    forces = masses * 9.81   ## 벡터 연산
    pylab.plot(forces, distances, 'bo',
               label = 'Measured displacements')
    pylab.title('Measured Displacement of Spring')
    pylab.xlabel('|Force| (Newtons)')
    pylab.ylabel('|Distance (meters)')

    #find linear fit
    a,b = pylab.polyfit(forces, distances, 1)  ## forces, distances 데이터에 일치하는 1차 방정식 ax+b=0
    predictedDistances = a * pylab.array(forces) + b   ## 벡터 연산
    k = 1.0/a
    pylab.plot(forces, predictedDistances,
               label = 'Displacements predicted by\nlinear fit, k = '
               + str(round(k, 5)))
    pylab.legend()

위 파이썬 프로그램을 실행하는 방법은 유사하다.

>>> fitData(path)
pylab.show()

앞에서 그렸던 스프링 실험 데이터를 힘과 거리로 가공한 데이터를 점 그래프로 보여주고 이 데이터와 일치하는 1차 방정식을 선으로 그린 그래프를 함께 보여준다. 먼저 1차 방정식 (ax+b=0의 a와 b)를 polyfit 함수로 계산하고, forces 리스트를 x축 좌표로 두어 y축 좌표 predictedDistances 리스트를 만든다. 이때 스프링 상수 k는 1/a이다. (k = - F/x이므로)

이번에는 스프링 실험 데이터와 일치하는 3차 방정식을 구해보자.

### fitData 버전 2:
def fitData(inputFile):
    masses, distances = getData(inputFile)

    masses = pylab.array(masses)
    distances = pylab.array(distances)
    
    forces = masses * 9.81  ## 벡터 연산
    pylab.plot(forces, distances, 'bo',
               label = 'Measured displacements')
    pylab.title('Measured Displacement of Spring')
    pylab.xlabel('|Force| (Newtons)')
    pylab.ylabel('|Distance (meters)')

    #find linear fit
    a,b = pylab.polyfit(forces, distances, 1)  ## forces, distances 데이터에 일치하는 1차 방정식 ax+b=0
    predictedDistances = a * pylab.array(forces) + b   ## 벡터 연산
    k = 1.0/a
    pylab.plot(forces, predictedDistances,
               label = 'Displacements predicted by\nlinear fit, k = '
               + str(round(k, 5)))

    #find cubic fit
    a,b,c,d = pylab.polyfit(forces, distances, 3) ## forces, distances 데이터에 일치하는 3차 방정식 ax^3+bx^2+cx+d=0
    predictedDistances = a * forces ** 3 + b * forces ** 2 + c * forces + d  ## 벡터 연산
    pylab.plot(forces, predictedDistances, 'b:', label = 'cubic fit')

    pylab.legend()

1차 방정식 그래프와 3차 방정식 그래프를 보면 3차 방정식이 스프링 데이터의 점 그래프와 더 일치하는 것을 알 수 있다. 하지만 3차 방정식 그래프를 보면 맨 오른쪽에서 x축의 힘이 증가하는데도 y축의 스프링이 늘어난 거리가 줄어들고 있다. 이것은 우리가 가지고 있는 스프링에 대한 상식과 맞지 않다.

사실 19개의 스프링 데이터 중에서 마지막 6개 데이터는 측정시 오류가 있었던 잘못된 데이터였다. 이 6개의 데이터에 맞추느라 3차 방정식이 1차 방정식보다 더 일치하는 것처럼 보인 것 뿐이다. 실험 데이터를 해석할 때 이러한 상황이 자주 발생하는데 이를 과도한 일치(overfitting)라고 부른다.

이 6개의 데이터를 제외하고 다시 1차 방정식과 3차 방정식을 다시 계산해서 그려보자.

### fitData 버전 3 (최종):
def fitData(inputFile):
    masses, distances = getData(inputFile)

    masses = pylab.array(masses[:-6])    ## 파일에서 읽은 masses 데이터 중에서 마지막 6개는 제외
    distances = pylab.array(distances[:-6])  ## 파일에서 읽은 distances 데이터 중에서 마지막 6개는 제외
    
    forces = masses * 9.81   ## 벡터 연산
    pylab.plot(forces, distances, 'bo',
               label = 'Measured displacements')
    pylab.title('Measured Displacement of Spring')
    pylab.xlabel('|Force| (Newtons)')
    pylab.ylabel('|Distance (meters)')

    #find linear fit
    a,b = pylab.polyfit(forces, distances, 1)  ## forces, distances 데이터에 일치하는 1차 방정식 ax+b=0
    predictedDistances = a * pylab.array(forces) + b   ## 벡터 연산
    k = 1.0/a
    pylab.plot(forces, predictedDistances,
               label = 'Displacements predicted by\nlinear fit, k = '
               + str(round(k, 5)))

    #find cubic fit
    a,b,c,d = pylab.polyfit(forces, distances, 3) ## forces, distances 데이터에 일치하는 3차 방정식 ax^3+bx^2+cx+d=0
    predictedDistances = a * forces ** 3 + b * forces ** 2 + c * forces + d   ## 벡터 연산
    pylab.plot(forces, predictedDistances, 'b:', label = 'cubic fit')

    pylab.legend()

최종 버전의 fitData 함수를 실행하면 3차 방정식이 1차 방정식과 거의 동일하게 일치함을 알 수 있다. 물론 실험 데이터에 대한 법칙을 알고 있고, 마지막 6개의 데이터를 제외하는 일종의 휴리스틱을 적용했기 때문에 이러한 결과를 얻을 수 있었다. 하지만 실험 데이터를 제대로 해석하기 위해서는 많은 경험이 필요하다.

마지막으로 파이썬 프로그램에서 사용한 pylab 벡터 연산에 대해서 설명하고 이 강의를 마친다. 파이썬 모듈 pylab은 pylab.array 클래스를 제공하고, 이 클래스로 표현된 객체는 일종의 벡터로 다룰 수 있다. 예를 들어 pylab.array 객체 [1 2 3]은 파이썬 리스트 [1,2,3]과 다르다. 각각을 출력했을 때 콤마가 있는지 여부를 따져서 두 유형을 구분할 수 있다.

pylab.array 클래스 형의 객체들을 벡터 연산으로 다룰 수 있다. pylab.array 클래스 형으로 표현된 객체 [1 2 3]과 [4 5 6]을 더하면 첫번째 숫자 1과 4를 더해 5, 2와 5를 더해 7, 3과 6을 더해 9를 계산하고 다시 이 클래스 형으로 표현된 객체 [5 7 9]를 만든다. 이를 벡터 덧셈 연산이라 부른다.

참고로 벡터 연산에 반대로 스칼라 연산이라는 용어를 사용한다. 예를 들어, 1 + 4 = 5는 단일 값(스칼라) 1과 또 다른 스칼라 4를 더해 새로운 스칼라 5를 계산한다.

앞에서 소개한 파이썬 프로그램에서 pylab.array의 쓰임새를 이해할 수 있는 예제를 다음과 같이 정리하였다. pylab.array( list )는 주어진 리스트 list를 pylab.array 클래스형 객체로 변환하는 함수이다.

>>> masses,distances = getData(datafile)
>>> print(masses)
[0.1, 0.15, 0.203873598, 0.244648318, 0.30479103, 0.326197757, 0.387359837, 0.448521916, 0.508664628, 0.540265036, 0.612640163, 0.673802243, 0.713659531, 0.739041794, 0.775739042, 0.838939857, 0.9, 0.95, 1.0]
>>> type(masses)
<class 'list'>

>>> massesArr = pylab.array(masses)
>>> distancesArr = pylab.array(distances)
>>> print(massesArr)
[0.1        0.15       0.2038736  0.24464832 0.30479103 0.32619776
 0.38735984 0.44852192 0.50866463 0.54026504 0.61264016 0.67380224
 0.71365953 0.73904179 0.77573904 0.83893986 0.9        0.95
 1.        ]
>>> type(massesArr)
<class 'numpy.ndarray'>
>>> type(distancesArr)
<class 'numpy.ndarray'>

>>> forcesArr = massesArr * 9.81
>>> type(forcesArr)
<class 'numpy.ndarray'>
>>> print(forcesArr)
[0.981  1.4715 2.     2.4    2.99   3.2    3.8    4.4    4.99   5.3
 6.01   6.61   7.001  7.25   7.61   8.23   8.829  9.3195 9.81  ]
>>>