本節我們主要簡單介紹機器學習常用的語言–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有了一定的瞭解，各種封裝包特別好用有木有，下一節將正式進入機器學習。