机器学习（4）：python基础及fft、svd、股票k线图、分形等实践

本节我们主要简单介绍机器学习常用的语言–python。楼主本身是写java的，在这之前对python并不了解，接触之后发现python比java简直要好用几千倍。这里主要通过常用的统计量、fft、股票k线图及分形等样例，介绍python的使用及各种包的加载。

1、常用的统计量

常用的统计量实践中有很多，比如均值、方差等，这里主要介绍偏度、峰度及其代码实现。

偏度：是衡量随机变量概率分布的不对称性，是相对于均值不对称程度的度量。偏度为负，则表示概率密度函数在均值左侧的尾部比在右侧的长，也即长尾在左侧；偏度为正则刚好相反。偏度为零，表示数值相对均匀的分布在均值两侧，但并不一定是对称的。

峰度：是概率密度在均值处峰值高低的特征。一般定义正太分布的峰度为0，定义峰度为四阶中心矩阵除以方差的平方减3（减3是为了让正态分布的峰度为0）。

代码样例：我们随机生成1000个数，计算其偏度、峰度，并画出对应的分布图。

import numpy as np
from scipy import stats
import math
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

def calc_statistics(x):
n = x.shape[0]  # 样本个数

# 手动计算
m = 0
m2 = 0
m3 = 0
m4 = 0
for t in x:
m += t
m2 += t*t
m3 += t**3
m4 += t**4
m /= n
m2 /= n
m3 /= n
m4 /= n

mu = m
sigma = np.sqrt(m2 - mu*mu)
skew = (m3 - 3*mu*m2 + 2*mu**3) / sigma**3
kurtosis = (m4 - 4*mu*m3 + 6*mu*mu*m2 - 4*mu**3*mu + mu**4) / sigma**4 - 3
print '手动计算均值、标准差、偏度、峰度：', mu, sigma, skew, kurtosis

# 使用系统函数验证
mu = np.mean(x, axis=0)
sigma = np.std(x, axis=0)
skew = stats.skew(x)
kurtosis = stats.kurtosis(x)
return mu, sigma, skew, kurtosis

if __name__ == '__main__':
d = np.random.randn(1000)
print d
mu, sigma, skew, kurtosis = calc_statistics(d)
print '函数库计算均值、标准差、偏度、峰度：', mu, sigma, skew, kurtosis
# 一维直方图
mpl.rcParams[u'font.sans-serif'] = 'SimHei'
mpl.rcParams[u'axes.unicode_minus'] = False
y1, x1, dummy = plt.hist(d, bins=50, normed=True, color='g', alpha=0.75)
t = np.arange(x1.min(), x1.max(), 0.05)
y = np.exp(-t**2 / 2) / math.sqrt(2*math.pi)
plt.plot(t, y, 'r-', lw=2)
plt.title(u'高斯分布，样本个数：%d' % d.shape[0])
plt.grid(True)
plt.show()

d = np.random.randn(1000, 2)
mu, sigma, skew, kurtosis = calc_statistics(d)
print '函数库计算均值、标准差、偏度、峰度：', mu, sigma, skew, kurtosis

# 二维图像
N = 30
density, edges = np.histogramdd(d, bins=[N, N])
print '样本总数：', np.sum(density)
density /= density.max()
x = y = np.arange(N)
t = np.meshgrid(x, y)
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, projection='3d')
ax.scatter(t[0], t[1], density, c='r', s=15*density, marker='o', depthshade=True)
ax.plot_surface(t[0], t[1], density, cmap=cm.Accent, rstride=2, cstride=2, alpha=0.9, lw=0.75)
ax.set_xlabel(u'X')
ax.set_ylabel(u'Y')
ax.set_zlabel(u'Z')
plt.title(u'二元高斯分布，样本个数：%d' % d.shape[0], fontsize=20)
plt.tight_layout(0.1)
plt.show()

结果输出如下：

[ 1.02372188e+00 -7.94235491e-01 -9.98698986e-01 …, -9.98193952e-01 5.28687536e-01 -1.89125334e+00]

手动计算均值、标准差、偏度、峰度： 0.0284912703006 0.99569007385 -0.0905769761327 0.175907983403

函数库计算均值、标准差、偏度、峰度： 0.0284912703006 0.99569007385 -0.0905769761327 0.175907983403

python提供了直接计算的封装包，和我们手动计算的结果一致，如果有需求直接调用封装包即可。

画出的分布图和二元高斯分布图如下：





2、fft

fft的主要功能是进行时域频域信号转换，比如三角波信号从频域恢复到时域，如下图示：



对应代码如下（同时提供锯齿波等各种转换）：

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

def triangle_wave(size, T):
t = np.linspace(-1, 1, size, endpoint=False)
# where
# y = np.where(t < 0, -t, 0)
# y = np.where(t >= 0, t, y)
y = np.abs(t)
y = np.tile(y, T) - 0.5
x = np.linspace(0, 2*np.pi*T, size*T, endpoint=False)
return x, y

def sawtooth_wave(size, T):
t = np.linspace(-1, 1, size)
y = np.tile(t, T)
x = np.linspace(0, 2*np.pi*T, size*T, endpoint=False)
return x, y

def triangle_wave2(size, T):
x, y = sawtooth_wave(size, T)
return x, np.abs(y)

def non_zero(f):
f1 = np.real(f)
f2 = np.imag(f)
eps = 1e-4
return f1[(f1 > eps) | (f1 < -eps)], f2[(f2 > eps) | (f2 < -eps)]

if __name__ == "__main__":
mpl.rcParams['font.sans-serif'] = [u'simHei']
mpl.rcParams['axes.unicode_minus'] = False
np.set_printoptions(suppress=True)

x = np.linspace(0, 2*np.pi, 16, endpoint=False)
print '时域采样值：', x
y = np.sin(2*x) + np.sin(3*x + np.pi/4) + np.sin(5*x)
# y = np.sin(x)

N = len(x)
print '采样点个数：', N
print '\n原始信号：', y
f = np.fft.fft(y)
print '\n频域信号：', f/N
a = np.abs(f/N)
print '\n频率强度：', a

iy = np.fft.ifft(f)
print '\n逆傅里叶变换恢复信号：', iy
print '\n虚部：', np.imag(iy)
print '\n实部：', np.real(iy)
print '\n恢复信号与原始信号是否相同：', np.allclose(np.real(iy), y)

plt.subplot(211)
plt.plot(x, y, 'go-', lw=2)
plt.title(u'时域信号', fontsize=15)
plt.grid(True)
plt.subplot(212)
w = np.arange(N) * 2*np.pi / N
print u'频率采样值：', w
plt.stem(w, a, linefmt='r-', markerfmt='ro')
plt.title(u'频域信号', fontsize=15)
plt.grid(True)
plt.show()

# 三角/锯齿波
x, y = triangle_wave(20, 5)
# x, y = sawtooth_wave(20, 5)
N = len(y)
f = np.fft.fft(y)
# print '原始频域信号：', np.real(f), np.imag(f)
print '原始频域信号：', non_zero(f)
a = np.abs(f / N)

# np.real_if_close
f_real = np.real(f)
eps = 0.3 * f_real.max()
print eps
f_real[(f_real < eps) & (f_real > -eps)] = 0
f_imag = np.imag(f)
eps = 0.3 * f_imag.max()
print eps
f_imag[(f_imag < eps) & (f_imag > -eps)] = 0
f1 = f_real + f_imag * 1j
y1 = np.fft.ifft(f1)
y1 = np.real(y1)
# print '恢复频域信号：', np.real(f1), np.imag(f1)
print '恢复频域信号：', non_zero(f1)

plt.figure(figsize=(8, 8), facecolor='w')
plt.subplot(311)
plt.plot(x, y, 'g-', lw=2)
plt.title(u'三角波', fontsize=15)
plt.grid(True)
plt.subplot(312)
w = np.arange(N) * 2*np.pi / N
plt.stem(w, a, linefmt='r-', markerfmt='ro')
plt.title(u'频域信号', fontsize=15)
plt.grid(True)
plt.subplot(313)
plt.plot(x, y1, 'b-', lw=2, markersize=4)
plt.title(u'三角波恢复信号', fontsize=15)
plt.grid(True)
plt.tight_layout(1.5, rect=[0, 0.04, 1, 0.96])
plt.suptitle(u'快速傅里叶变换FFT与频域滤波', fontsize=17)
plt.show()

3、SVD

在上一节已经简单介绍过SVD，可以说svd是处理图像重要特征的一个重要手段。给出一个图片，通过svd找出其最主要的特征，并基于特征对图像进行还原。这里主要给出svd的代码实现。代码如下：

import numpy as np
import os
from PIL import Image
import matplotlib.pyplot as plt
import matplotlib as mpl
from pprint import pprint

def restore1(sigma, u, v, K):  # 奇异值、左特征向量、右特征向量
m = len(u)
n = len(v[0])
a = np.zeros((m, n))
for k in range(K):
uk = u[:, k].reshape(m, 1)
vk = v[k].reshape(1, n)
a += sigma[k] * np.dot(uk, vk)
a[a < 0] = 0
a[a > 255] = 255
# a = a.clip(0, 255)
return np.rint(a).astype('uint8')

def restore2(sigma, u, v, K):  # 奇异值、左特征向量、右特征向量
m = len(u)
n = len(v[0])
a = np.zeros((m, n))
for k in range(K+1):
for i in range(m):
a[i] += sigma[k] * u[i][k] * v[k]
a[a < 0] = 0
a[a > 255] = 255
return np.rint(a).astype('uint8')

if __name__ == "__main__":
A = Image.open("me.png", 'r')
print A
output_path = r'.\Pic'
if not os.path.exists(output_path):
os.mkdir(output_path)
a = np.array(A)
print a.shape
K = 50
u_r, sigma_r, v_r = np.linalg.svd(a[:, :, 0])
u_g, sigma_g, v_g = np.linalg.svd(a[:, :, 1])
u_b, sigma_b, v_b = np.linalg.svd(a[:, :, 2])
plt.figure(figsize=(10,10), facecolor='w')
mpl.rcParams['font.sans-serif'] = [u'simHei']
mpl.rcParams['axes.unicode_minus'] = False
for k in range(1, K+1):
print k
R = restore1(sigma_r, u_r, v_r, k)
G = restore1(sigma_g, u_g, v_g, k)
B = restore1(sigma_b, u_b, v_b, k)
I = np.stack((R, G, B), axis=2)
Image.fromarray(I).save('%s\\svd_%d.png' % (output_path, k))
if k <= 12:
plt.subplot(3, 4, k)
plt.imshow(I)
plt.axis('off')
plt.title(u'奇异值个数：%d' % k)
plt.suptitle(u'SVD与图像分解', fontsize=20)
plt.tight_layout(0.3, rect=(0, 0, 1, 0.92))
# plt.subplots_adjust(top=0.9)
plt.show()

使用不同个数特征值还原后的图片如下：



4、股票k线图

玩股票的人都知道，股票k线图是描述股票涨跌的一个直观体现，如果知道某只股票一定时间内的交易数据，如何汇出对应的股票k线图呢？代码如下：

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.finance import candlestick_ohlc

if __name__ == "__main__":
mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = False

np.set_printoptions(suppress=True, linewidth=100, edgeitems=5)
data = np.loadtxt('SH600000.txt', dtype=np.float, delimiter='\t', skiprows=2, usecols=(1, 2, 3, 4))
data = data[:50]
N = len(data)

t = np.arange(1, N+1).reshape((-1, 1))
data = np.hstack((t, data))

fig, ax = plt.subplots(facecolor='w')
fig.subplots_adjust(bottom=0.2)
candlestick_ohlc(ax, data, width=0.6, colorup='r', colordown='g', alpha=0.9)
plt.xlim((0, N+1))
plt.grid(b=True)
plt.title(u'股票K线图', fontsize=18)
plt.tight_layout(2)
plt.show()

5、分形

第一次知道分形，是看电影奇异博士时发现的，瞬间觉得这个东东太神奇，忍不住想尝试下。

#!/usr/bin/python
# -*- coding:utf-8 -*-
__author__ = 'xuena'

import numpy as np
import pylab as pl
import time
from matplotlib import cm

def iter_point(c):
z = c
for i in xrange(1, 100): # 最多迭代100次
if abs(z)>2: break # 半径大于2则认为逃逸
z = z*z+c
return i # 返回迭代次数

def draw_mandelbrot(cx, cy, d):
"""
绘制点(cx, cy)附近正负d的范围的Mandelbrot
"""
x0, x1, y0, y1 = cx-d, cx+d, cy-d, cy+d
y, x = np.ogrid[y0:y1:200j, x0:x1:200j]
c = x + y*1j
start = time.clock()
mandelbrot = np.frompyfunc(iter_point,1,1)(c).astype(np.float)
print "time=",time.clock() - start
pl.imshow(mandelbrot, cmap=cm.jet, extent=[x0,x1,y0,y1])
pl.gca().set_axis_off()

x,y = 0.27322626, 0.595153338

pl.subplot(231)
draw_mandelbrot(-0.5,0,1.5)
for i in range(2,7):
pl.subplot(230+i)
draw_mandelbrot(x, y, 0.2**(i-1))
pl.subplots_adjust(0.02, 0, 0.98, 1, 0.02, 0)
pl.show()

很神奇的分形图如下：



通过实践，我们对python有了一定的了解，各种封装包特别好用有木有，下一节将正式进入机器学习。