import numpy as np
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 matplotlib.pyplot as plt
from scipy.stats import norm, poisson
import math
import scipy
from scipy.interpolate import BarycentricInterpolator
from scipy.interpolate import CubicSpline

a = np.arange(0, 60, 10).reshape((-1, 1)) + np.arange(6)
print (a)

1.使用array創(chuàng)建

標(biāo)準(zhǔn)Python的列表(list)中，元素本質(zhì)是對象洋侨。
如：L = [1, 2,
3]舍扰，需要3個(gè)指針和三個(gè)整數(shù)對象，對于數(shù)值運(yùn)算比較浪費(fèi)內(nèi)存和CPU希坚。
因此边苹，Numpy提供了ndarray(N-dimensional array
object)對象：存儲(chǔ)單一數(shù)據(jù)類型的多維數(shù)組。

# 通過array函數(shù)傳遞list對象
L = [1, 2, 3, 4, 5, 6]
print ("L = ", L)
a = np.array(L)  # 數(shù)組沒有逗號
print ("a = ", a)
print (type(a), type(L))

# 若傳遞的是多層嵌套的list裁僧，將創(chuàng)建多維數(shù)組
b = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])
print (b)

# # # # # 數(shù)組大小可以通過其shape屬性獲得
print (a.shape)
print (b.shape)

# # 也可以強(qiáng)制修改shape
b.shape = 4, 3
print (b)
# # # 注：從(3,4)改為(4,3)并不是對數(shù)組進(jìn)行轉(zhuǎn)置个束，而只是改變每個(gè)軸的大小慕购，數(shù)組元素在內(nèi)存中的位置并沒有改變

# # 當(dāng)某個(gè)軸為-1時(shí)，將根據(jù)數(shù)組元素的個(gè)數(shù)自動(dòng)計(jì)算此軸的長度
b.shape = 2, -1
print (b)
print (b.shape)
# # #
b.shape = 3, 4
print (b)

# # # 使用reshape方法茬底，可以創(chuàng)建改變了尺寸的新數(shù)組沪悲，原數(shù)組的shape保持不變
c = b.reshape((4, -1))
print ("b = \n", b)
print ('c = \n', c)

# # # 數(shù)組b和c共享內(nèi)存，修改任意一個(gè)將影響另外一個(gè)
b[0][1] = 20
print ("b = \n", b)
print ("c = \n", c)

# # # 數(shù)組的元素類型可以通過dtype屬性獲得
print (a.dtype)
print (b.dtype)

# # # # 可以通過dtype參數(shù)在創(chuàng)建時(shí)指定元素類型
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安全的轉(zhuǎn)換
f = d.astype(np.int)
print (f)

# # 但不要強(qiáng)制僅修改元素類型殿如，如下面這句，將會(huì)以int來解釋單精度float類型
d.dtype = np.int
print (d)

2.使用函數(shù)創(chuàng)建

如果生成一定規(guī)則的數(shù)據(jù)捶枢，可以使用NumPy提供的專門函數(shù)
arange函數(shù)類似于python的range函數(shù)：指定起始值握截、終止值和步長來創(chuàng)建數(shù)組
和Python的range類似，arange同樣不包括終值烂叔；但arange可以生成浮點(diǎn)類型谨胞，而range只能是整數(shù)類型

a = np.arange(1, 10, 0.5)
print (a)

# # # linspace函數(shù)通過指定起始值、終止值和元素個(gè)數(shù)來創(chuàng)建數(shù)組蒜鸡，缺省包括終止值
b = np.linspace(1, 10, 10)
print ('b = ', b)

# 可以通過endpoint關(guān)鍵字指定是否包括終值
c = np.linspace(1, 10, 10, endpoint=False)
print ('c = ', c)

# # 和linspace類似胯努，logspace可以創(chuàng)建等比數(shù)列
# 下面函數(shù)創(chuàng)建起始值為10^1，終止值為10^2逢防，有10個(gè)數(shù)的等比數(shù)列
d = np.logspace(1, 2, 9, endpoint=True)
print (d)

# # # 下面創(chuàng)建起始值為2^0叶沛，終止值為2^10(包括)，有10個(gè)數(shù)的等比數(shù)列
f = np.logspace(0, 10, 10, endpoint=True, base=2)
print (f)

# # # 使用 frombuffer, fromstring, fromfile等函數(shù)可以從字節(jié)序列創(chuàng)建數(shù)組
s = 'abcdz'
g = np.fromstring(s, dtype=np.int8) # 復(fù)制對應(yīng)的ASII碼
print (g)

3.存取

# 3.1常規(guī)辦法：數(shù)組元素的存取方法和Python的標(biāo)準(zhǔn)方法相同
a = np.arange(10)
print (a)

# # 獲取某個(gè)元素
print (a[3])
# # # # 切片[3,6)忘朝，左閉右開
print (a[3:6])
# 省略開始下標(biāo)灰署，表示從0開始
print (a[:5])
# 下標(biāo)為負(fù)表示從后向前數(shù)
print (a[3:])
# 步長為2
print (a[1:9:2])

# # # # 步長為-1，即翻轉(zhuǎn)
print (a[::-1])

# # # # 切片數(shù)據(jù)是原數(shù)組的一個(gè)視圖局嘁，與原數(shù)組共享內(nèi)容空間溉箕，可以直接修改元素值
a[1:4] = 10, 20, 30
print (a)

# 因此，在實(shí)踐中悦昵，切實(shí)注意原始數(shù)據(jù)是否被破壞肴茄，如：
b = a[2:5]
b[0] = 200
print (a)

# 3.2 整數(shù)/布爾數(shù)組存取
# 3.2.1
# 根據(jù)整數(shù)數(shù)組存取：當(dāng)使用整數(shù)序列對數(shù)組元素進(jìn)行存取時(shí)但指，
# 將使用整數(shù)序列中的每個(gè)元素作為下標(biāo)寡痰，整數(shù)序列可以是列表(list)或者數(shù)組(ndarray)。
# 使用整數(shù)序列作為下標(biāo)獲得的數(shù)組不和原始數(shù)組共享數(shù)據(jù)空間棋凳。
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
# 使用布爾數(shù)組i作為下標(biāo)存取數(shù)組a中的元素：返回?cái)?shù)組a中所有在數(shù)組b中對應(yīng)下標(biāo)為True的元素
# 生成10個(gè)滿足[0,1)中均勻分布的隨機(jī)數(shù)
a = np.random.rand(10)
print (a)
# 大于0.5的元素索引
print (a > 0.5)
# 大于0.5的元素
b = a[a > 0.5]
print (b)
# 將原數(shù)組中大于0.5的元素截取成0.5
a[a > 0.5] = 0.5
print (a)
# # # b不受影響
print (b)

# 3.3 二維數(shù)組的切片
a = np.arange(0, 60, 10)    # 行向量
print ('a = ', a)
b = a.reshape((-1, 1))      # 轉(zhuǎn)換成列向量
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)

# 二維數(shù)組的切片
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數(shù)學(xué)庫的時(shí)間比較
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直接使用庫函數(shù)
a = np.array((1, 2, 3, 4, 5, 5, 7, 3, 2, 2, 8, 8))
print ('原始數(shù)組：', a)
# # 使用庫函數(shù)unique
b = np.unique(a)
print ('去重后：', b)

# # 4.2.2 二維數(shù)組的去重，結(jié)果會(huì)是預(yù)期的么剩岳？
c = np.array(((1, 2), (3, 4), (5, 6), (1, 3), (3, 4), (7, 6)))
print (u'二維數(shù)組：\n', c)
print ('去重后：', np.unique(c))

# # 4.2.3 方案1：轉(zhuǎn)換為虛數(shù)
# r, i = np.split(c, (1, ), axis=1)
# x = r + i * 1j
x = c[:, 0] + c[:, 1] * 1j
print ('轉(zhuǎn)換成虛數(shù)：', x)
print ('虛數(shù)去重后：', np.unique(x))
print (np.unique(x, return_index=True))   # 思考return_index的意義
idx = np.unique(x, return_index=True)[1]
print ('二維數(shù)組去重：\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))

a = np.arange(1, 10).reshape(3,3)
print (a)
b = a + 10
print (b)
print (np.dot(a, b)) # dot 正常的矩陣乘法
print (a * b) # * 對應(yīng)元素的相乘

a = np.arange(1, 10)
print (a)
b = np.arange(20,25)
print (b)
print (np.concatenate((a, b)))

5.繪圖

# 5.1 繪制正態(tài)分布概率密度函數(shù)
# 自定義字體使中文正常顯示
mpl.rcParams['font.sans-serif'] = [u'SimHei']  #FangSong/黑體 FangSong/KaiTi
mpl.rcParams['axes.unicode_minus'] = False

mu = 0
sigma = 1
x = np.linspace(mu - 3 * sigma, mu + 3 * sigma, 51)
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.figure(facecolor='w') # 背景色設(shè)置為白色
# plt.plot(x, y, 'ro-', linewidth=2)
plt.plot(x, y, 'r-', x, y, 'go', linewidth=2, markersize=8)
plt.xlabel('X', fontsize=15)
plt.ylabel('Y', fontsize=15)
plt.title(u'高斯分布函數(shù)', fontsize=18)
plt.grid(True)
plt.show()

# 5.2 損失函數(shù)：Logistic損失(-1,1)/SVM Hinge損失/ 0/1損失
plt.figure(figsize=(10,8),dpi=100) # 指定圖像尺寸和dpi
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.savefig('1.png')
plt.show()

# 5.3 x^x

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

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), facecolor='w')
plt.plot(y, x, 'r-', linewidth=2)
plt.grid(True)
plt.title(u'胸型線', fontsize=20)
# plt.savefig('breast.png')
plt.show()

# 5.5 心形線
t = np.linspace(0, 2*np.pi, 100)
x = 16 * np.sin(t) ** 3
y = 13 * np.cos(t) - 5 * np.cos(2*t) - 2 * np.cos(3*t) - 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, 'r-', linewidth=2)
plt.grid()
plt.show()

# Bar
x = np.arange(0, 10, 0.1)
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 = np.random.rand(10000)
t = np.arange(len(x))
#plt.hist(x, 30, color='#000000', alpha=0.5, label=u'均勻分布') # alpha透明度
plt.plot(t, x, 'g.', label=u'均勻分布') # .小點(diǎn), o 圈, - 線
plt.legend(loc='upper left')
plt.grid()
plt.show()

# 6.2 驗(yàn)證中心極限定理
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)

N = 1000
M = 10000
plt.subplot(122)
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 插值
rv = poisson(5)
x1 = a[1]
y1 = rv.pmf(x1)
itp = BarycentricInterpolator(x1, y1)  # 重心插值
x2 = np.linspace(x.min(), x.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]
# print(x,y)
u = np.linspace(-3, 3, 101)
x, y = np.meshgrid(u, u)
z = x*y*np.exp(-(x**2 + y**2)/2) / math.sqrt(2*math.pi)
# z = x*y*np.exp(-(x**2 + y**2)/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'])]