Numpy math

Numpy 배열을 사용하기 위한 자료는 google finance에서 호출한 kospi 일일 주가이다.

웹에서 자료를 호출하고 호출된 자료를 조절하기 위해 다음의 모듈이 필요하다.

import numpy as np
import pandas as pd
from datetime import datetime
from pandas_datareader import data, wb #웹에서 데이터 호출
import pandas_datareader.data as web
import matplotlib.pyplot as plt #그래프 작성을 위해 필요
%matplotlib inline
from sklearn.linear_model import LinearRegression

np.set_printoptions(precision=3)

startD=datetime(2016, 1,1)
endD=datetime(2017, 3, 31)
d=web.DataReader('KRX:kospi', 'google', startD, endD)
d.tail(3)

	Open	High	Low	Close	Volume
Date
2017-03-29	2172.31	2172.31	2162.04	2166.98	626792000
2017-03-30	2170.18	2174.16	2159.16	2164.64	643563000
2017-03-31	2166.62	2166.93	2159.80	2160.23	441640000

d1=np.array(d)
d1[:3]
array([[  1.954e+03,   1.955e+03,   1.919e+03,   1.919e+03,   3.590e+08],
      [  1.912e+03,   1.938e+03,   1.912e+03,   1.931e+03,   4.465e+08],
      [  1.934e+03,   1.934e+03,   1.912e+03,   1.925e+03,   5.946e+08]])

수학함수

NumPy는 sin, cos, exp, 제곱근와 같은 익숙한 수학 함수들로의 변환을 위해

np.sin(), np.cos(), np.exp(), np.sqrt()등의 함수를 제공한다.

이 함수들은 배열에 적용시 원소 별로 연산된다.

np.sin(d1[:5])
array([[ 0.389,  0.434,  0.684,  0.684,  0.775],
      [ 0.964,  0.712,  0.964,  1.   , -0.174],
      [-0.825, -0.825,  0.999,  0.358,  0.136],
      [-0.614, -0.576, -0.545,  0.501,  0.781],
      [-0.969,  0.953, -0.907,  0.949,  0.762]])

np.sum(배열, axis), np.prod(배열,axis): 배열의 합, 곱을 반환.

axis 설정하지 않으면 모든 원소들의 합, 곱을 반환

np.cumsum(배열, axis), np.cumprod(배열axis): 배열의 누적 합, 누적 곱을 반환.

axis 설정하지 않으면 모든 원소들의 누적합, 누적곱을 반환

np.argesort(배열, axis): 배열을 오름차순으로 정렬할 경우 각 값의 인덱스를 반환 np.sort(배열, axis): 오름차순으로 정렬된 새로운 배열을 반환(행의 배열이 디폴트)

np.sum(d1, axis=0) #객체의 각 열의 합
array([  6.188e+05,   6.209e+05,   6.161e+05,   6.187e+05,   1.171e+11])

np.prod(d1[:5], axis=0)
array([  2.616e+16,   2.707e+16,   2.512e+16,   2.605e+16,   1.612e+43])

np.cumsum(d1[:5], axis=0)
array([[  1.954e+03,   1.955e+03,   1.919e+03,   1.919e+03,   3.590e+08],
      [  3.866e+03,   3.892e+03,   3.831e+03,   3.849e+03,   8.056e+08],
      [  5.801e+03,   5.826e+03,   5.742e+03,   5.775e+03,   1.400e+09],
      [  7.716e+03,   7.753e+03,   7.644e+03,   7.679e+03,   1.793e+09],
      [  9.606e+03,   9.671e+03,   9.527e+03,   9.597e+03,   2.223e+09]])

#다음 데이터는 20개의 수를 랜덤하게 추출한 것이다.
x=np.floor(np.random.randn(20)*10)
np.argsort(x)
array([11,  7, 14, 17,  0, 13, 19,  2,  4,  6, 18,  1, 16,  9, 12,  3, 15,
       8, 10,  5], dtype=int64)

np.sort(x)
array([-20., -14., -12., -11.,  -8.,  -7.,  -4.,  -1.,   2.,   2.,   3.,
        8.,   9.,  10.,  12.,  12.,  13.,  14.,  16.,  18.])

np.all(배열객체) : 배열의 원소중에 0이 없으면 True를 반환,

NaN은 0이 아니므로 True np.apply_along_axis(함수, axis, 함수의 인자들) : 배열에 함수를 적용시킬 때 사용한다. axis의 값을 즉, 0 (열), 1(행)을 전달할 경우 정수만을 입력해야 한다.

(언어 R에서 apply()함수와 동일)

np.all(d1)
True

#객체의 최대값과 최소값의 합의 평균
def myfun(a):
   return (np.max(a)+np.min(a))/2
#위 kospi 일일 주가 자료의 각 열의 최대값과 최소값의 평균을 계산한다.
np.apply_along_axis(myfun, 0, d)
array([  2.011e+03,   2.020e+03,   1.994e+03,   2.007e+03,   5.137e+08])

np.zeros(차원): 0으로 구성된 행렬을 생성하며 그 차원은 수로 전달된다.

np.eye(N, M): 항등행렬을 생성하는 함수로 N은 대각원소(1)의 숫자,

N은 첫원소 1의 열의 위치 np.nonzero(객체): 객체에서 행단위로 0이 아닌 원소의 인덱스를 반환

np.empty(차원): np.zeros()와 같은 결과를 보이지만 각 원소의 값은 초기화되지 않음 np.trace(객체): 객체 중 대각원소의 합(행렬에서 trace) 반환

np.transpose(객체) : 행과 열을 전환 (전치행렬)=numpy배열.T

np.where(조건): 조건에 만족하는 값의 인덱스를 반환

np.fromfunction(함수, shape) : shape를 함수의 인수로 배열을 생성.

np.zeros((2,2))
array([[ 0.,  0.],
      [ 0.,  0.]])

x=np.eye(3)
x
array([[ 1.,  0.,  0.],
      [ 0.,  1.,  0.],
      [ 0.,  0.,  1.]])

np.nonzero(x)
array([0, 1, 2], dtype=int64), array([0, 1, 2], dtype=int64))

y=np.empty((2,3))#각 원소는 의미없는 값이다.
y
array([[ 0.,  0.,  0.],
      [ 0.,  0.,  0.]])

m=range(1, 7)
n=0
for i in [0,1]:
   for j in [0,1,2]:
       y[i,j]=m[n]
       n +=1
y
array([[ 1.,  2.,  3.],
      [ 4.,  5.,  6.]])

np.trace(y)
6.0

np.transpose(y)
array([[ 1.,  4.],
      [ 2.,  5.],
      [ 3.,  6.]])

#주가 자료인 d1의 0번째 열이 2160 이상인 값의 인덱스를 반환
np.where(d1[:,0]>=2160)
(array([298, 299, 301, 302, 303, 304, 305, 306, 307], dtype=int64),)

d1[np.where((d1[:,0]>=2150) & (d1[:,0]<=2170)), :]
#위의 경우 처럼 두개 이상의 조건을 만들 때는 하나의 조건마다 '()'으로 경계를 주어야 한다.
array([[[  2.155e+03,   2.157e+03,   2.143e+03,   2.150e+03,   4.024e+08],
       [  2.150e+03,   2.165e+03,   2.149e+03,   2.165e+03,   3.028e+08],
       [  2.160e+03,   2.161e+03,   2.152e+03,   2.157e+03,   2.649e+08],
       [  2.163e+03,   2.182e+03,   2.160e+03,   2.178e+03,   3.039e+08],
       [  2.160e+03,   2.170e+03,   2.159e+03,   2.168e+03,   3.944e+08],
       [  2.161e+03,   2.165e+03,   2.153e+03,   2.156e+03,   3.959e+08],
       [  2.167e+03,   2.169e+03,   2.158e+03,   2.163e+03,   3.940e+08],
       [  2.167e+03,   2.167e+03,   2.160e+03,   2.160e+03,   4.416e+08]]])

x=np.array([0,1,2,3,4,5,6,7])
y=np.array([3,4,5,6,7,8,9,10])
np.where(((x>=3) &(x<=5)))
(array([3, 4, 5], dtype=int64),)

다음의 예는 파이썬의 lambda를 이용하여 두 인수 i, j 값을 입력받아 그 합을 반환하는 함수이다.

또한 이렇게 생성된 배열은 행과 열이 3, 3인 배열 이다. 그 배열의 각 자리의 위치가 i,j의 인수가 된다.

아래의 예에서 그 인수는 (0,0), (0,1),(0,2), (1,0), (1,1), (1,2), (2,0) (2,1)(2,2)가 된다.

그러므로 람다함수에 전달되는 인수의 개수와 shape로 전달되는 튜플의 개수는 같아야 한다.

np.fromfunction(lambda i, j:i+j, (3,3))
array([[ 0.,  1.,  2.],
      [ 1.,  2.,  3.],
      [ 2.,  3.,  4.]])

행렬과 벡터 연산

np.dot(arr1, rr2): 행렬인 두 객체의 행렬곱

이 함수는 np.matmul(arr1, arr2)와 같다.

np.vdot(arr1, arr2) : 각 배열(행렬)의 산술 곱, 즉, 동일한 인덱스의 원소들의 곱 np.inner(arr1, arr2): 각 1차원 두 배열(벡터)들의 내적(1차원 배열의 행렬곱) np.outer(arr1, arr2): 두 벡터의 외적으로 계산.

외적의 계산은 다음과 같다.

위의 외적은 np.kron(arr1, arr2) 함수로도 계산할 수 있다.

그러나 이 함수는 행벡터로 반환된다.

np.diag(arr1): 배열의 대각원소들을 반환 np.linalg.matrix_power(arr, n): 배열(행렬)을 n번 행렬곱을 시행한 결과 반환

이 함수는 (from numpy import linalg as LA 인스톨 후에

LA.matrix_power() 실행) 예)

LA.matrix_power(x, 2)=

a=np.arange(6).reshape((2,3))
b=np.arange(6).reshape((3,2))
print(a)
print(b)

[[0 1 2]
[3 4 5]]
[[0 1]
[2 3]
[4 5]]

np.dot(a,b) #행렬곱
array([[10, 13],
      [28, 40]])

np.matmul(a, b)
array([[10, 13],
      [28, 40]])

a=np.arange(4).reshape((2,2))
b=np.arange(4).reshape((2,2))
print(a)
print(b)

[[0 1]
[2 3]]
[[0 1]
[2 3]]

np.vdot(a, b) #0*0+1*1+2*2+3*3=14, 행렬들의 산술곱연산
14

np.inner([3, 2,7],[5, 3, 6])
63

#위의 내적 연산은 다음과 같이 행렬 연산과 같다.
a=np.array([3,2,7])
b=np.transpose([5,3,6])
np.dot(a, b)
63

a=np.array([1,2,3])
b=np.array([4, 5])
print(a)
print(b)

[1 2 3]
[4 5]

np.outer(a, b) #두벡터의 외적을 계산한다.
array([[ 4,  5],
      [ 8, 10],
      [12, 15]])

#위의 외적의 연산 방식은 다음과 같이 이루어진다.
re=np.zeros((3,2))
for i in range(len(a)):
   for j in range(len(b)):
       re[i, j]=a[i]*b[j]
re        
array([[  4.,   5.],
      [  8.,  10.],
      [ 12.,  15.]])

np.kron(a, b)
array([ 4,  5,  8, 10, 12, 15])

a=np.arange(4).reshape((2,2))
a
array([[0, 1],
      [2, 3]])

from numpy import linalg as LA
print(np.dot(a, a))
LA.matrix_power(a, 2) #a를 두번 행렬곱으로 위 함수 적용과 같다.

[[ 2  3]
[ 6 11]]
array([[ 2,  3],
      [ 6, 11]])

대각행렬 생성

np.fill_diagonal(객체, 대각원소): 객체의 대각원소가 지정된 값으로 전환된다. np.diag_indices(n): 대각원소 수(n)를 인수로 전달하여 그 위치를 지정한다.

x=np.eye(3)
np.fill_diagonal(x, [1,2,3])
x
array([[ 1.,  0.,  0.],
      [ 0.,  2.,  0.],
      [ 0.,  0.,  3.]])

x=np.eye(3)
di=np.diag_indices(3)
x[di]=[4,2,8]
x
array([[ 4.,  0.,  0.],
      [ 0.,  2.,  0.],
      [ 0.,  0.,  8.]])

분해

1. 촐레스키분해

정방행렬이면서 대칭행렬(원행렬과 전치행렬이 동일한 행렬)인 행렬의 하삼각행렬(L)과

상삼각행렬(U)의 곱은 원행렬과 같다는 성질을 이용하여 행렬을 분해한다. np.linalg.cholesky(arr1): 배열을 촐레스키분해를 실시하여 하삼각행렬을 반환

행렬 A=LU 이므로

                                    

이 된다.

a=np.array([[2,1],[1,2]])
print(a)

[[2 1]
[1 2]]
b=np.linalg.cholesky(a)
b

array([[ 1.414,  0.   ],
      [ 0.707,  1.225]])

np.dot(np.linalg.inv(b), a) #np.linalg.inv(x) 역행렬을 계산한다.
#결과중 인덱스[1,0]의 값은 0에 대략화 된다. 그러므로 결과는 상삼각행렬이 된다.:

array([[  1.414e+00,   7.071e-01],
      [  1.110e-16,   1.225e+00]])

2) qr분해

내적이 0인 두 벡터로 구성된 행렬을 직교행렬이라 한다. 즉, 두 벡터의 내적은

 이고 위의 각은 두 벡터의 사이각이므로 이 값이 0이 된다는 것은 두 벡터의     

  사이각이 0 즉 직교한다는 것이다. 이 직교행렬의 전치행렬은 원행렬의       역행렬과 같다.

          즉, 임의의 직교행렬 Q에 대해 다음이 성립한다.

이러한 직교성은 선형적으로 독립임을 나타낸다.

즉, 직교인 두 벡터는 서로 배수의 관계가 아님을 나타낸다.

또한 이 직교성을 갖는 두 벡터의 거리가 1인 집합을 정규직교(orthonomal)이라 한다.

어떤 행렬의 정규직교 행렬을 찾는 것은 많은 선형에 관계된 모형, 대수문제에 중요한 부분이다.

행렬 A가 대칭행렬이면 다음 식을 성립하는 정규직교 행렬 Q찾는 것은 가능하다. 정규직교행렬은 그람-슈미트 방법에 의해 계산된다.

위 식을 적용하여 원행렬을 정규직교행렬과 상삼각행렬로 분해할 수 있다. np.linalg.qr(arr):배열의 qr 분해를 계산하여 정규직교행렬(q)와 상삼각행렬(r)을 반환하는 함수

a = np.random.randn(9, 6)
a
array([[ 0.211,  0.332,  0.159, -0.297, -0.449,  0.833],
      [-0.221,  0.954, -0.017,  1.534, -0.364, -1.546],
      [ 0.582,  0.819, -1.738,  0.044,  0.243,  0.197],
      [ 0.603,  0.077,  1.111,  1.009, -0.446, -0.645],
      [ 0.161,  1.347, -1.118,  0.025, -0.228, -0.715],
      [-0.422,  0.86 ,  0.614,  0.473, -0.193, -1.321],
      [-0.881,  0.059, -0.986,  0.187, -0.39 , -0.482],
      [ 1.03 ,  0.197, -0.361,  0.396,  0.738, -0.074],
      [-2.032,  0.241, -0.405,  0.397,  0.268,  0.927]])

q, r=np.linalg.qr(a)
print(q) #직교정규행렬

[[-0.08  -0.161  0.14   0.37  -0.305  0.653]
[ 0.084 -0.456  0.229 -0.573 -0.14  -0.047]
[-0.22  -0.397 -0.55  -0.098  0.03   0.223]
[-0.228 -0.041  0.465 -0.413 -0.26   0.3  ]
[-0.061 -0.647 -0.15   0.323 -0.097 -0.104]
[ 0.16  -0.409  0.477  0.244  0.214 -0.391]
[ 0.334 -0.022 -0.372 -0.27  -0.522 -0.271]
[-0.39  -0.102 -0.136 -0.321  0.585  0.045]
[ 0.77  -0.101 -0.048 -0.12   0.386  0.439]]

print(r) #상삼각행렬

[[-2.639  0.038 -0.218  0.199 -0.176  0.322]
[ 0.    -2.086  1.197 -1.005  0.293  1.447]
[ 0.     0.     2.386  0.835 -0.514 -1.026]
[ 0.     0.     0.    -1.513 -0.081  0.931]
[ 0.     0.     0.     0.     1.031  0.49 ]
[ 0.     0.     0.     0.     0.     1.593]]

np.inner(a[:,0],a[:,1])
-0.10146553997356389

3) 고유값 분해(eigenvalualbe decomposition)

정방행렬 A에서 0이 아닌 벡터 x와의 다음 관계가 성립된다면

x를 고유벡터(eigenvector), λλ를 고유값(eigenvalue)라고 한다. 고유값 분해는 정방행렬(행과 열의 수가 같은 행렬)에서만 적용가능

어떤 정방행렬에 대해 고유값(eigenvalue)와 고유벡터(eigen vector)를 계산할 수 있다.

이 고유값과 고유벡터에 의해 다음과 같이 원행렬을 나타낼 수 있다.

A: 정방행렬로서 행렬 A의 고유값(λi)들을 대각원소로 하는 대각행렬

P: 고유값에대응하는고유벡터

np.linalg.eig(arr): 정방행렬의 고유값과 고유벡터를 계산. 고유값, 고유벡터를 결과로 반환. np.linalg.eigh(arr): 에르미트행렬(자신과 켤레 전치가 같은 복속 정방행렬) 즉,

를 만족하는 복소행렬의 고유값과 고유벡터를 반환

a=np.array([[3,0],[8,-1]])
a
array([[ 3,  0],
      [ 8, -1]])

ev1, ev2=np.linalg.eig(a)
print(ev1)
print(ev2)
np.linalg.eigh(a)

[-1.  3.]
[[ 0.     0.447]
[ 1.     0.894]]
(array([-7.246,  9.246]), array([[ 0.615, -0.788],
       [-0.788, -0.615]]))

#위의 결과에서 P=ev2이므로 inv(P)AP=D 가 계산된다.
#D는 ev1을 대각원소로 하는 대각행렬이다.
np.dot(np.dot(np.linalg.inv(ev2), a), ev2)
array([[-1.,  0.],
      [ 0.,  3.]])

4) 특이값 분해(singular value decomposition)

행렬 A가 대칭해렬이 아닌 고유벡터 행렬을 가진다면 A의 고유값 분해는 계산할 수 없다. 그러나 행렬 A가 m×nm×n (m>n)의 차원의 행렬이라면

A에 특이값 분해 방법을 적용할 수 있다. 즉,

위 식에서 U,V는 직교 열들을 가진다. 즉, 다음이 성립된다.

U는 A의 왼쪽고유벡터로 이루어진 직교행렬 즉, AATAAT를 고유값분해로 직교대각화하여

얻은 직교행렬로 좌특이벡터(left singular vectors)라고 한다.

D는 대각성분이 의 고유값의 제곱근 즉,로 이루어진 대각행렬이다.

V는 오른쪽 고유백터로 이루어진 직교행렬 즉, 를 고유값 분해로 직교대각화하여

얻은 직교행렬로 우특이벡터(right singular vectors)라고 한다.

np.linalg.svd(객체): 특이값 분해를 실행하여 U, D, V를 반환

a=np.array([3,6,2,3,0,0, 0, 0]).reshape((4,2))
a
array([[3, 6],
      [2, 3],
      [0, 0],
      [0, 0]])

np.linalg.svd(a)
(array([[-0.882, -0.472,  0.   ,  0.   ],
       [-0.472,  0.882,  0.   ,  0.   ],
       [ 0.   ,  0.   ,  1.   ,  0.   ],
       [ 0.   ,  0.   ,  0.   ,  1.   ]]),
array([ 7.606,  0.394]),
array([[-0.472, -0.882],
       [ 0.882, -0.472]]))

np.linalg.norm(객체, axis=None): 행렬 또는 벡터의 놈(norm)을 반환

norm은 유클리스 거리를 나타내는 것으로 그 벡터의 길이를 나타낸 것이다.

두 벡터의 norm은 그 벡터들 사이의 거리를 나타낸다.

a=np.arange(8).reshape(4,2)
a
array([[0, 1],
      [2, 3],
      [4, 5],
      [6, 7]])

np.linalg.norm(a, axis=1)
array([ 1.   ,  3.606,  6.403,  9.22 ])

np.linalg.det(배열) : 배열의 행렬식을 계산

np.linalg.det([[2,5],[7,10]])
-15.0

a=np.arange(16).reshape(4,4)
np.linalg.det(a)
-2.9582283945787796e-30

np.linalg.matrix_rank(객체, tol=None): 객체의 rank를 반환

즉, 벡터 X, Y, Z들 사이의 다음과 같은 선형결합이 이루어진다.

cX+cY+cZ=0 이러한 선형결합을 유도할 수 있는 스칼라 c가 존재한다.

이 경우는 위의 세 벡터는 서로 의존관계에 있다.

그러나 c가 0을 제외하고 존재하지 않으면 이 세벡터는 서로 독립이 된다. 이러한 각 벡터의 독립성을 나타내는 지표가 rank이다. 다음 항등행렬은 모든 행 또는 열이 독립이다.

즉, 모든 계수가 0이외에 다른 수로 각 벡터가 만나는 지점을 발견할 수 없다.

그러므로 rank=3이 된다.

a=np.eye(3)
np.linalg.matrix_rank(a)
3

np.linalg.solve(객체1, 객체2): 객체1*계수=객체2인 방정식의 해를 계산한다.

다음 식들을 고려해보자.

3x+y=9

x+2y=8

a=np.array([[1,1],[1,2]])
b=np.array([9,8])
np.linalg.solve(a, b)
array([ 10.,  -1.])

np.linalg.lstsq(객체1(x), 객체2(y)): 최소자승법에 의한 선형 행렬방정식의 해를 계산

즉, 행렬 x,y 사이의 y=bx+e의 모형을 구축할 경우

실제 y와 추정 y 사이의 차이를 최소로 하는 조건을 만족시키는 과정을 최소자승법이라한다.

이러한 절차에 의해 선형모형의 계수행렬(b)를 계산한다.

이 함수는 계수와 잔차제곱합, rank, 특이값(singular value)을 반환한다.

x = np.array([0, 1, 2, 3])
x=np.c_[x, np.ones(len(x))]
y = np.array([-1, 0.2, 0.9, 2.1])
b, ssr, rank, sv=np.linalg.lstsq(x, y)
print(b)

[ 1.   -0.95]
print(ssr)

[ 0.05]
print(rank)

2
print(sv)
[ 4.1    1.091]

np.linalg.inv(객체): 역행렬을 계산

#기약 행 사다리꼴은 sympy 모듈을 이용하여 계산한다.
from sympy import *
a=Matrix(np.arange(3*4).reshape(3, 4)) #sympy Matrix 구조로 전환
a.rref()

(Matrix([
[1.0,   0, -1.0, -2.0],
[  0, 1.0,  2.0,  3.0],
[  0,   0,    0,    0]]), [0, 1])

행렬계산은 numpy모듈과 함께 sympy 모듈을 함계 사용할 수 있다.

from sympy import *

a=Matrix([[1,-1],[3,4], [0,2]])
a
Matrix([
[1, -1],
[3,  4],
[0,  2]])

#행렬의 차원을 나타냄
print(a.shape)

(3, 2)

#row(), col() 함수를 이용하여 지정된 열, 행으로 접근할 수 있다.
print(a.row(1))

Matrix([[3, 4]])
print(a.col(1))
Matrix([[-1], [4], [2]])

#지정된 행과 열을 제거하기 위해 col_del(), row_del()을 사용한다.
a=Matrix(np.arange(3*4).reshape(3, 4))
print(a)

Matrix([[0.0, 1.00000000000000, 2.00000000000000, 3.00000000000000], [4.00000000000000, 5.00000000000000, 6.00000000000000, 7.00000000000000], [8.00000000000000, 9.00000000000000, 10.0000000000000, 11.0000000000000]])
print(a.col_del(0))

None
print(a.row_del(0))
None

a

Matrix([
[5.0,  6.0,  7.0],
[9.0, 10.0, 11.0]])

#기약 행 사다리꼴은 sympy 모듈을 이용하여 계산한다.
from sympy import *
a=Matrix(np.arange(3*4).reshape(3, 4)) #sympy Matrix 구조로 전환
a.rref()
(Matrix([
[1.0,   0, -1.0, -2.0],
[  0, 1.0,  2.0,  3.0],
[  0,   0,    0,    0]]), [0, 1])