!/usr/bin/python
-- coding:utf-8 --
import numpy as np
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import time
from scipy.optimize import leastsq
from scipy import stats
import scipy.optimize as opt
import scipy
import matplotlib.pyplot as plt
from scipy.stats import norm, poisson
from scipy.interpolate import BarycentricInterpolator
from scipy.interpolate import CubicSpline
import math
import seaborn
def residual(t, x, y):
return y - (t[0] * x ** 2 + t[1] * x + t[2])
def residual2(t, x, y):
print(t[0], t[1])
return y - t[0]np.sin(t[1]x)
x ** x x > 0
(-x) ** (-x) x < 0
def f(x):
y = np.ones_like(x)
i = x > 0
y[i] = np.power(x[i], x[i])
i = x < 0
y[i] = np.power(-x[i], -x[i])
return y
if name == "main":
# # 開場白：
# numpy是非常好用的數據包，如：可以這樣得到這個二維數組
# [[ 0 1 2 3 4 5]
# [10 11 12 13 14 15]
# [20 21 22 23 24 25]
# [30 31 32 33 34 35]
# [40 41 42 43 44 45]
# [50 51 52 53 54 55]]
# a = np.arange(0, 60, 10).reshape((-1, 1)) + np.arange(6)
# print a
# 正式開始 -:)
# 標准Python的列表(list)中，元素本質是對象。
# 如：L = [1, 2, 3]，需要3個指針和三個整數對象，對於數值運算比較浪費內存和CPU。
# 因此，Numpy提供了ndarray(N-dimensional array object)對象：存儲單一數據類型的多維數組。
# # 1.使用array創建
# 通過array函數傳遞list對象
# L = [1, 2, 3, 4, 5, 6]
# print("L = ", L)
# a = np.array(L)
# print("a = ", a)
# print(type(a))
# # # 若傳遞的是多層嵌套的list，將創建多維數組
# b = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])
# print(b)
# # #
# # # # 數組大小可以通過其shape屬性獲得
# print(a.shape)
# print(b.shape)
# # #
# # # 也可以強制修改shape
# b.shape = 4, 3
# print(b)
# # # # 注：從(3,4)改為(4,3)並不是對數組進行轉置，而只是改變每個軸的大小，數組元素在內存中的位置並沒有改變
# # #
# # # 當某個軸為-1時，將根據數組元素的個數自動計算此軸的長度
# b.shape = 2, -1
# print(b)
# print(b.shape)
# # #
# b.shape = 3, 4
# # # 使用reshape方法，可以創建改變了尺寸的新數組，原數組的shape保持不變
# c = b.reshape((4, -1))
# print("b = \n", b
# print('c = \n', c)
# #
# # # 數組b和c共享內存，修改任意一個將影響另外一個
# b[0][1] = 20
# print("b = \n", b)
# print("c = \n", c)
# # #
# # # 數組的元素類型可以通過dtype屬性獲得
# print(a.dtype)
# print(b.dtype)
# # #
# # # # 可以通過dtype參數在創建時指定元素類型
# d = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]], dtype=np.float)
# f = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]], dtype=np.complex)
# print(d)
# print(f)
# # #
# # # 如果更改元素類型，可以使用astype安全的轉換
# f = d.astype(np.int)
# print(f)
# #
# # # 但不要強制僅修改元素類型，如下面這句，將會以int來解釋單精度float類型
# d.dtype = np.int
# print(d)
#
# 2.使用函數創建
# 如果生成一定規則的數據，可以使用NumPy提供的專門函數
# arange函數類似於python的range函數：指定起始值、終止值和步長來創建數組
# 和Python的range類似，arange同樣不包括終值；但arange可以生成浮點類型，而range只能是整數類型
# a = np.arange(1, 10, 0.5)
# print(a)
# # #
# # # linspace函數通過指定起始值、終止值和元素個數來創建數組，缺省包括終止值
# b = np.linspace(1, 10, 10)
# print('b = ', b)
# # #
# # 可以通過endpoint關鍵字指定是否包括終值
# c = np.linspace(1, 10, 10, endpoint=False)
# print('c = ', c)
# # #
# # 和linspace類似，logspace可以創建等比數列
# # 下面函數創建起始值為10^1，終止值為10^2，有20個數的等比數列
# d = np.logspace(1, 2, 10, endpoint=True)
# print(d)
# # #
# # # 下面創建起始值為2^0，終止值為2^10(包括)，有10個數的等比數列
# f = np.logspace(0, 10, 11, endpoint=True, base=2)
# print(f)
# # #
# # # 使用 frombuffer, fromstring, fromfile等函數可以從字節序列創建數組
# s = 'abcd'
# g = np.fromstring(s, dtype=np.int8)
# print(g)
# #
# 3.存取
# 3.1常規辦法：數組元素的存取方法和Python的標准方法相同
# a = np.arange(10)
# print(a)
# # # 獲取某個元素
# print(a[3])
# # # # 切片[3,6)，左閉右開
# print(a[3:6])
# # # # 省略開始下標，表示從0開始
# print(a[:5])
# # # # 下標為負表示從后向前數
# print(a[3:])
# # # # 步長為2
# print(a[1:9:2])
# # # # 步長為-1，即翻轉
# print(a[::-1])
# # # # 切片數據是原數組的一個視圖，與原數組共享內容空間，可以直接修改元素值
# a[1:4] = 10, 20, 30
# print(a)
# # # # 因此，在實踐中，切實注意原始數據是否被破壞，如：
# b = a[2:5]
# b[0] = 200
# print(a)
#
# # 3.2 整數/布爾數組存取
# # 3.2.1
# 根據整數數組存取：當使用整數序列對數組元素進行存取時，
# 將使用整數序列中的每個元素作為下標，整數序列可以是列表(list)或者數組(ndarray)。
# 使用整數序列作為下標獲得的數組不和原始數組共享數據空間。
# a = np.logspace(0, 9, 10, base=2)
# print(a)
# i = np.arange(0, 10, 2)
# print i
# # # # 利用i取a中的元素
# b = a[i]
# print b
# # # b的元素更改，a中元素不受影響
# b[2] = 1.6
# print b
# print a
# # 3.2.2
# 使用布爾數組i作為下標存取數組a中的元素：返回數組a中所有在數組b中對應下標為True的元素
# # 生成10個滿足[0,1)中均勻分布的隨機數
# a = np.random.rand(10)
# print a
# # # 大於0.5的元素索引
# print a > 0.5
# # # # 大於0.5的元素
# b = a[a > 0.5]
# print b
# # # # 將原數組中大於0.5的元素截取成0.5
# a[a > 0.5] = 0.5
# print a
# # # # b不受影響
# print b
# 3.3 二維數組的切片
# [[ 0 1 2 3 4 5]
# [10 11 12 13 14 15]
# [20 21 22 23 24 25]
# [30 31 32 33 34 35]
# [40 41 42 43 44 45]
# [50 51 52 53 54 55]]
# a = np.arange(0, 60, 10) # 行向量
# print 'a = ', a
# b = a.reshape((-1, 1)) # 轉換成列向量
# print b
# c = np.arange(6)
# print c
# f = b + c # 行 + 列
# print f
# # 合並上述代碼：
# a = np.arange(0, 60, 10).reshape((-1, 1)) + np.arange(6)
# print a
# # 二維數組的切片
# print a[[0, 1, 2], [2 ,3, 4]]
# print a[4, [2, 3, 4]]
# print a[4:, [2, 3, 4]]
# i = np.array([True, False, True, False, False, True])
# print a[i]
# print a[i, 3]
# # 4.1 numpy與Python數學庫的時間比較
# for j in np.logspace(0, 7, 10):
# j = int(j)
# x = np.linspace(0, 10, j)
# start = time.clock()
# y = np.sin(x)
# t1 = time.clock() - start
#
# x = x.tolist()
# start = time.clock()
# for i, t in enumerate(x):
# x[i] = math.sin(t)
# t2 = time.clock() - start
# print j, ": ", t1, t2, t2/t1
# 4.2 元素去重
# 4.2.1直接使用庫函數
# a = np.array((1, 2, 3, 4, 5, 5, 7, 3, 2, 2, 8, 8))
# print '原始數組：', a
# # 使用庫函數unique
# b = np.unique(a)
# print '去重后：', b
# # 4.2.2 二維數組的去重，結果會是預期的么？
# c = np.array(((1, 2), (3, 4), (5, 6), (1, 3), (3, 4), (7, 6)))
# print '二維數組', c
# print '去重后：', np.unique(c)
# # 4.2.3 方案1：轉換為虛數
# # r, i = np.split(c, (1, ), axis=1)
# # x = r + i * 1j
# x = c[:, 0] + c[:, 1] * 1j
# print '轉換成虛數：', x
# print '虛數去重后：', np.unique(x)
# print np.unique(x, return_index=True) # 思考return_index的意義
# idx = np.unique(x, return_index=True)[1]
# print '二維數組去重：\n', c[idx]
# # 4.2.3 方案2：利用set
# print '去重方案2：\n', np.array(list(set([tuple(t) for t in c])))
# # 4.3 stack and axis
# a = np.arange(1, 10).reshape((3, 3))
# b = np.arange(11, 20).reshape((3, 3))
# c = np.arange(101, 110).reshape((3, 3))
# print 'a = \n', a
# print 'b = \n', b
# print 'c = \n', c
# print 'axis = 0 \n', np.stack((a, b, c), axis=0)
# print 'axis = 1 \n', np.stack((a, b, c), axis=1)
# print 'axis = 2 \n', np.stack((a, b, c), axis=2)
# 5.繪圖
# 5.1 繪制正態分布概率密度函數
# mu = 0
# sigma = 1
# x = np.linspace(mu - 3 * sigma, mu + 3 * sigma, 50)
# y = np.exp(-(x - mu) ** 2 / (2 * sigma ** 2)) / (math.sqrt(2 * math.pi) * sigma)
# print x.shape
# print 'x = \n', x
# print y.shape
# print 'y = \n', y
# # plt.plot(x, y, 'ro-', linewidth=2)
# plt.plot(x, y, 'r-', x, y, 'go', linewidth=2, markersize=8)
# plt.grid(True)
# plt.show()
mpl.rcParams['font.sans-serif'] = [u'Bold'] #FangSong/黑體 FangSong/KaiTi
mpl.rcParams['axes.unicode_minus'] = False
# # 5.2 損失函數：Logistic損失(-1,1)/SVM Hinge損失/ 0/1損失
=============================================================================
x = np.array(np.linspace(start=-2, stop=3, num=1001, dtype=np.float))
y_logit = np.log(1 + np.exp(-x)) / math.log(2)
y_boost = np.exp(-x)
y_01 = x < 0
y_hinge = 1.0 - x
y_hinge[y_hinge < 0] = 0
plt.plot(x, y_logit, 'r-', label='Logistic Loss', linewidth=2)
plt.plot(x, y_01, 'g-', label='0/1 Loss', linewidth=2)
plt.plot(x, y_hinge, 'b-', label='Hinge Loss', linewidth=2)
plt.plot(x, y_boost, 'm--', label='Adaboost Loss', linewidth=2)
plt.grid()
plt.legend(loc='upper right')
plt.title(u"SVM loss 函數圖像")
plt.savefig('1.png')
plt.show()
=============================================================================
# # 5.3 x^x
=============================================================================
x = np.linspace(-1.3, 1.3, 101)
y = f(x)
plt.plot(x, y, 'g-', label='x^x', linewidth=2)
plt.grid()
plt.legend(loc='upper left')
plt.show()
=============================================================================
# # 5.4 胸型線
=============================================================================
x = np.arange(1, 0, -0.001)
y = (-3 * x * np.log(x) + np.exp(-(40 * (x - 1 / np.e)) ** 4) / 25) / 2
plt.figure(figsize=(5,7))
plt.plot(y, x, 'r-', linewidth=2)
plt.grid(True)
plt.show()
=============================================================================
# 5.5 心形線
=============================================================================
t = np.linspace(0, 7, 100)
x = 16 * np.sin(t) ** 3
y = 13 * np.cos(t) - 5 * np.cos(2t) - 2 * np.cos(3t) - np.cos(4*t)
plt.plot(x, y, 'r-', linewidth=2)
plt.grid(True)
plt.show()
=============================================================================
# # 5.6 漸開線
=============================================================================
t = np.linspace(0, 50, num=1000)
x = t*np.sin(t) + np.cos(t)
y = np.sin(t) - t*np.cos(t)
plt.plot(x, y, 'y+', linewidth=2)
plt.grid()
plt.show()
=============================================================================
# # Bar
=============================================================================
mpl.rcParams['font.sans-serif'] = [u'SimHei'] #黑體 FangSong/KaiTi
mpl.rcParams['axes.unicode_minus'] = False
x = np.arange(0, 10, 0.05)
y = np.sin(x)
plt.bar(x, y, width=0.04, linewidth=0.2)
plt.plot(x, y, 'r--', linewidth=2)
plt.title(u'Sin曲線')
plt.xticks(rotation=-60)
plt.xlabel('X')
plt.ylabel('Y')
plt.grid()
plt.show()
=============================================================================
# # 6. 概率分布
# # 6.1 均勻分布
=============================================================================
x = 100*np.random.rand(10000)
t = np.arange(len(x))
plt.hist(x, 30, color='m', alpha=0.5, label=u'均勻分布')
plt.plot(t, x, 'r-', label=u'均勻分布')
plt.legend(loc='upper left')
plt.grid()
plt.show()
=============================================================================
# # 6.2 驗證中心極限定理
=============================================================================
t = 1000
a = np.zeros(10000)
for i in range(t):
a += np.random.uniform(-5, 5, 10000)
a /= t
plt.hist(a, bins=30, color='g', alpha=0.5, normed=True, label=u'均勻分布疊加')
plt.legend(loc='upper left')
plt.grid()
plt.show()
=============================================================================
# 6.21 其他分布的中心極限定理
=============================================================================
lamda = 10
p = stats.poisson(lamda)
y = p.rvs(size=1000)
mx = 30
r = (0, mx)
bins = r[1] - r[0]
plt.figure(figsize=(10, 8), facecolor='w')
plt.subplot(121)
plt.hist(y, bins=bins, range=r, color='g', alpha=0.8, normed=True)
t = np.arange(0, mx+1)
plt.plot(t, p.pmf(t), 'ro-', lw=2)
plt.grid(True)
=============================================================================
#
=============================================================================
lamda=10
N = 1000
M = 10000
plt.subplot(111)
a = np.zeros(M, dtype=np.float)
p = stats.poisson(lamda)
for i in np.arange(N):
y = p.rvs(size=M)
a += y
a /= N
plt.hist(a, bins=20, color='g', alpha=0.8, normed=True)
plt.grid(b=True)
plt.show()
=============================================================================
# # 6.3 Poisson分布
=============================================================================
x = np.random.poisson(lam=5, size=10000)
print(x)
pillar = 15
a = plt.hist(x, bins=pillar, normed=True, range=[0, pillar], color='g', alpha=0.5)
plt.grid()
plt.show()
print(a)
print(a[0].sum())
=============================================================================
# # 6.4 直方圖的使用
=============================================================================
mu = 2
sigma = 3
data = mu + sigma * np.random.randn(1000)
h = plt.hist(data, 30, normed=1, color='#a0a0ff')
x = h[1]
y = norm.pdf(x, loc=mu, scale=sigma)
plt.plot(x, y, 'r--', x, y, 'ro', linewidth=2, markersize=4)
plt.grid()
plt.show()
=============================================================================
# # 6.5 插值
=============================================================================
x = np.random.poisson(lam=5, size=10000)
print(x)
pillar = 15
a = plt.hist(x, bins=pillar, normed=True, range=[0, pillar], color='g', alpha=0.5)
rv = poisson(5)
x1 = a[1]
y1 = rv.pmf(x1)
itp = BarycentricInterpolator(x1, y1) # 重心插值
x2 = np.linspace(x1.min(), x1.max(), 50)
y2 = itp(x2)
cs = scipy.interpolate.CubicSpline(x1, y1) # 三次樣條插值
plt.plot(x2, cs(x2), 'm--', linewidth=5, label='CubicSpine') # 三次樣條插值
plt.plot(x2, y2, 'g-', linewidth=3, label='BarycentricInterpolator') # 重心插值
plt.plot(x1, y1, 'r-', linewidth=1, label='Actural Value') # 原始值
plt.legend(loc='upper right')
plt.grid()
plt.show()
=============================================================================
# 7. 繪制三維圖像
=============================================================================
x, y = np.ogrid[-3:3:100j, -3:3:100j]
u = np.linspace(-3, 3, 101)
x, y = np.meshgrid(u, u)
z = xynp.exp(-(x2 + y2)/2) / math.sqrt(2*math.pi)
z = xynp.exp(-(x2 + y2)/2) / math.sqrt(2*math.pi)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z, rstride=5, cstride=5, cmap=cm.coolwarm, linewidth=0.1)
ax.plot_surface(x, y, z, rstride=5, cstride=5, cmap=cm.Accent, linewidth=0.5)
plt.show()
cmaps = [('Perceptually Uniform Sequential',
['viridis', 'inferno', 'plasma', 'magma']),
('Sequential', ['Blues', 'BuGn', 'BuPu',
'GnBu', 'Greens', 'Greys', 'Oranges', 'OrRd',
'PuBu', 'PuBuGn', 'PuRd', 'Purples', 'RdPu',
'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd']),
('Sequential (2)', ['afmhot', 'autumn', 'bone', 'cool',
'copper', 'gist_heat', 'gray', 'hot',
'pink', 'spring', 'summer', 'winter']),
('Diverging', ['BrBG', 'bwr', 'coolwarm', 'PiYG', 'PRGn', 'PuOr',
'RdBu', 'RdGy', 'RdYlBu', 'RdYlGn', 'Spectral',
'seismic']),
('Qualitative', ['Accent', 'Dark2', 'Paired', 'Pastel1',
'Pastel2', 'Set1', 'Set2', 'Set3']),
('Miscellaneous', ['gist_earth', 'terrain', 'ocean', 'gist_stern',
'brg', 'CMRmap', 'cubehelix',
'gnuplot', 'gnuplot2', 'gist_ncar',
'nipy_spectral', 'jet', 'rainbow',
'gist_rainbow', 'hsv', 'flag', 'prism'])]
=============================================================================
# 8.1 scipy
# 線性回歸例1
=============================================================================
x = np.linspace(-2, 2, 50)
A, B, C = 2, 3, -1
y = (A * x ** 2 + B * x + C) + np.random.rand(len(x))*0.75
t = leastsq(residual, [0, 0, 0], args=(x, y))
theta = t[0]
print('真實值：', A, B, C)
print('預測值：', theta)
y_hat = theta[0] * x ** 2 + theta[1] * x + theta[2]
plt.plot(x, y, 'r-', linewidth=2, label=u'Actual')
plt.plot(x, y_hat, 'g-', linewidth=2, label=u'Predict')
plt.legend(loc='upper left')
plt.grid()
plt.show()
=============================================================================
# # 線性回歸例2
x = np.linspace(0, 5, 100)
A = 5
w = 1.5
y = A * np.sin(w*x) + np.random.rand(len(x)) - 0.5
t = leastsq(residual2, [3, 1], args=(x, y))
theta = t[0]
print('真實值：', A, w)
print('預測值：', theta)
y_hat = theta[0] * np.sin(theta[1] * x)
plt.plot(x, y, 'r-', linewidth=2, label='Actual')
plt.plot(x, y_hat, 'g-', linewidth=2, label='Predict')
plt.legend(loc='lower left')
plt.grid()
plt.show()
# # 8.2 使用scipy計算函數極值
# a = opt.fmin(f, 1)
# b = opt.fmin_cg(f, 1)
# c = opt.fmin_bfgs(f, 1)
# print a, 1/a, math.e
# print b
# print c
# marker description
# ”.” point
# ”,” pixel
# “o” circle
# “v” triangle_down
# “^” triangle_up
# “<” triangle_left
# “>” triangle_right
# “1” tri_down
# “2” tri_up
# “3” tri_left
# “4” tri_right
# “8” octagon
# “s” square
# “p” pentagon
# “*” star
# “h” hexagon1
# “H” hexagon2
# “+” plus
# “x” x
# “D” diamond
# “d” thin_diamond
# “|” vline
# “_” hline
# TICKLEFT tickleft
# TICKRIGHT tickright
# TICKUP tickup
# TICKDOWN tickdown
# CARETLEFT caretleft
# CARETRIGHT caretright
# CARETUP caretup
# CARETDOWN caretdown