科學計算與數據可視化1

科學計算與數據可視化1

@(程序設計)

numpy模塊

• Numpy最重要的一個特色就是其N維數組對象（即ndarray）該對象是一個快速而靈活的大數據集容器。
• 使用Numpy，開發人員能夠執行如下操做：

一、數組的算數和邏輯運算。

二、傅立葉變換和用於圖形操做的例程。

三、與線性代數有關的操做。

numpy.ndarray()

• ndarray是Numpy的數組類，其中的全部元素必須是相同的數據類型。ndarray類的重要對象屬性有：

• 利用array函數，能夠將序列類型的對象（元組、列表和其餘數組）轉換成數組類型ndarray。

  • ndarray.ndim：數組維度。兩行及以上爲二維
  • ndarray.shape：表示數組各維度大小的元組。
  • ndarray.size：數組元素的總個數，等於shape屬性中元組元素的乘積。
  • ndarray.dtype：數組中元素的數據類型
  • numpy.ndarray()就是numpy的構造函數，咱們可使用這個函數建立一個ndarray對象。構造函數有以下幾個可選參數：

np.savetxt()

• 將array保存到txt文件，並保持原格式

np.loadtxt()

• loadtxt(fname, dtype=<class 'float'>, comments='#', delimiter=None, converters=None, skiprows=0, usecols=None, unpack=False, ndmin=0)
• np.loadtxt()用於從文本加載數據。文本文件中的每一行必須含有相同的數據。

arange()

• numpy.arange([start, ]stop, [step, ]dtype=None)

arange函數是numpy內置的相似range的函數，其返回的是數組對象，而不是列表。

• range()與arrange()的區別

一、range()和np.arange()的返回類型不一樣，range()返回的是range；
二、object，而np.arange()返回的是ndarray類型；
三、range()不支持步長爲小數，而np.arange()支持步長(step)爲小數；
四、range()和np.arange()均可用於迭代；
五、range()和np.arange()都有三個參數，以第一個參數爲起點，第三個參數爲步長，截止到第二個參數以前的不包括第二個參數的數據序列。
五、range()可用於迭代，而np.arange做用遠不止於此，它是一個序列，可被當作向量使用。

array()

• array(...)
  array(object, dtype=None, copy=True, order='K', subok=False, ndmin=0)
• 數組的特色：

一、數組是相同數據類型的元素的集合。
二、數組中的各元素的存儲是有前後順序的，它們在內存中按照這個前後順序連續存放在一塊兒。
三、數組元素用整個數組的名字和它本身在數組中的順序位置來表示。例如，a[0]表示名字爲a的數組中的第一個元素，a[1]表明數組a的第二個元素，以此類推。

Array[0:]　　——>切片從前面序號「0」開始到結尾，包括「0」位
　　[2, 3, 9, 1, 4, 7, 6, 8]
Array[:-1]　　——>切片從後面序號「-1」到最前，不包括「-1」位
　　[2, 3, 9, 1, 4, 7, 6]
Array[3:-2]　　——>切從前面序號「3」開始(包括)到從後面序號「-2」結束(不包括)
　　[1, 4, 7]
Array[3::2]　　——>從前面序號「3」(包括)到最後，其中分隔爲「2」
　　[1, 7, 8]
Array[::2]　　——>從整列表中切出，分隔爲「2」
　　[2, 9, 4, 6]
Array[3::]　　——>從前面序號「3」開始到最後，沒有分隔
　　[1, 4, 7, 6, 8]
Array[3::-2]　　——>從前面序號「3」開始，往回數第二個，由於分隔爲「-2」
　　[1, 3]　　
Array[-1]　　　　——>此爲切出最後一個
　　8
Array[::-1]　　　　——>此爲倒序
　　[8, 6, 7, 4, 1, 9, 3, 2]

ndarray()

建立二維數組、訪問數組對象屬性

import numpy as np
list1=[5,6.5,9,2,3,7.8,5.6,4.9]
arr1=np.array(list1)
print(arr1)
print(arr1.dtype)   #數組中元素的數據類型
print(arr1.ndim)    #數組維度。兩行及以上爲二維
print(arr1.shape)   #表示數組各維度大小的元組。行*列
print(arr1.size)    #數組元素的總個數，等於shape屬性中元組元素的乘積。
list2=[[1,2,3,4,5],[6,7,8,9,10]]
arr2=np.array(list2)
print(arr2)
print(arr2.dtype)   #數組中元素的數據類型
print(arr2.ndim)    #數組維度。兩行及以上爲二維
print(arr2.shape)   #表示數組各維度大小的元組。行*列
print(arr2.size)    #數組元素的總個數，等於shape屬性中元組元素的乘積。

[5. 6.5 9. 2. 3. 7.8 5.6 4.9]
float64
1
(8,)
8
[[ 1 2 3 4 5]
[ 6 7 8 9 10]]
int32
2
(2, 5)
10

建立指定數據類型的數組對象

import numpy as np
list2=[[1,2,3,4,5],[6,7,8,9,10]]
arr2=np.array(list2)
arr3=np.array([10,20,30,40],dtype=np.float64)
print(arr3)
arr4=arr2.astype(np.float64) #轉換數據類型 float->int
print(arr4)
print(arr4.dtype)

[10. 20. 30. 40.]
[[ 1. 2. 3. 4. 5.]
[ 6. 7. 8. 9. 10.]]
float64

練習eg13_array1.py()

產生3個數組l一、l二、l3，分別爲[0 1 2 3]、[0 2 4 6]、[0 3 6 9]，輸出l一、l二、l3以及(l1,l2,l3)，並將(l1,l2,l3)的內容經過np.savetxt方法存入aa.txt文件中，而後用np.loadtxt讀出數據並顯示。程序保存爲eg13_array1.py

import numpy as np
l1=np.arange(4)  #易錯
l2,l3=l1*2,l1*3
print('l1：',l1)
print('l2：',l1)
print('l3：',l1)
print((l1,l2,l3))
np.savetxt('aa.txt',(l1,l2,l3)) #易錯
aa=np.loadtxt('aa.txt')
print(aa)

l1 [0 1 2 3]
l2 [0 2 4 6]
l3 [0 3 6 9]
(array([0, 1, 2, 3]), array([0, 2, 4, 6]), array([0, 3, 6, 9]))
aa: [[0. 1. 2. 3.]
[0. 2. 4. 6.]
[0. 3. 6. 9.]]

思考：繼續eg13_array1.py

import numpy as np
l1=np.arange(4)
l2,l3=l1*2,l1*3
print(l1)
print(l2)
print(l3)
aa=np.loadtxt('aa.txt')
print('aa的類型：',aa.dtype)
print('aa的維度：',aa.ndim)
print('aa的總個數：',aa.size)
print('aa的形狀：',aa.shape) #reshape()數組對象中的方法，用於改變數組的形狀
print('bb:',aa.reshape(4,3))  #reshape()改成一個四維數組

[0 1 2 3]
[0 2 4 6]
[0 3 6 9]
aa的類型： float64
aa的維度： 2
aa的總個數： 12
aa的形狀： (3, 4)
bb: [[0. 1. 2.]
[3. 0. 2.]
[4. 6. 0.]
[3. 6. 9.]]

定義數組的類型

import numpy as np
d1=np.loadtxt('bb.txt')
print('d1:',d1,'sum=',sum(d1))
d2=np.loadtxt('bb.txt',dtype=int)
print('d2:',d2,'sum=',sum(d2))
d3=np.loadtxt('bb.txt',dtype=str)
print('d3:',d3,'sum=',sum([eval(i) for i in d3]))  #這邊不是很明白

d1: [1. 2. 3. 4. 5.] sum= 15.0
d2: [1 2 3 4 5] sum= 15
d3: ['1' '2' '3' '4' '5'] sum= 15

astype()顯示轉換類型

import numpy as np
aa=np.array(['1.25','-9.6','42'],dtype=np.str)
bb=aa.astype(float)
print(bb)

[ 1.25 -9.6 42. ]

數組和標量之間的運算

import numpy as np
arr=np.array([[1,2,3],[4,5,6]])
print(arr)
print(arr*2)
print(1/arr)
print(arr-arr*2)

[[1 2 3]
[4 5 6]]
[[ 2 4 6]
[ 8 10 12]]
[[1. 0.5 0.33333333]
[0.25 0.2 0.16666667]]
[[-1 -2 -3]
[-4 -5 -6]]

數組的元素級運算與函數

• 大小相等的數組之間的任何算數運算都會應用到元素級

import numpy as np
arr1=np.arange(1,16).reshape(3,5)
print(arr1)
arr2=np.arange(1,30,2).reshape(3,5)
print(arr2)
print(arr1*arr2)
print(arr1/arr2)
print(arr1+arr2)
print(arr1-arr2)

[[ 1 2 3 4 5]
[ 6 7 8 9 10]
[11 12 13 14 15]]
[[ 1 3 5 7 9]
[11 13 15 17 19]
[21 23 25 27 29]]
[[ 1 6 15 28 45]
[ 66 91 120 153 190]
[231 276 325 378 435]]
[[1. 0.66666667 0.6 0.57142857 0.55555556]
[0.54545455 0.53846154 0.53333333 0.52941176 0.52631579]
[0.52380952 0.52173913 0.52 0.51851852 0.51724138]]
[[ 2 5 8 11 14]
[17 20 23 26 29]
[32 35 38 41 44]]
[[ 0 -1 -2 -3 -4]
[ -5 -6 -7 -8 -9]
[-10 -11 -12 -13 -14]]

二元函數及說明

一維數組的索引和切片（同列表）

• 數組切片是原始數組的視圖，數據並不會被複制，即視圖上的任何修改都會直接反映到源數組上。

讀csv文件

• c,v=np.loadtxt('data.csv', delimiter=',', usecols=(6,7), unpack=True,skiprows=1)

  • fname：讀取文件的文件名。例如C:/Dataset/iris.txt。
  • delimiter：數據之間的分隔符。如使用逗號","。
  • dtype：數據類型。如float，str等。
  • usecols：選取數據的列。
  • skiprows : list-like or integer, default None須要忽略的行數（從文件開始處算起），或須要跳過的行號列表（從0開始）。

平均值和加權平均值

加權平均值的概念

• np.mean()
• np.average()
• 對應的權值列表：weights=[]
• 將各數值乘以相應的權數，而後加總求和獲得整體值，再除以總的單位數。

import numpy as np
a=(70,80,60)
print(np.mean(a)) #平均值
print(np.average(a,weights=[3,3,4]))  #加權平均值  average()

70.0
69.0

VWAP(Volume-Weighted Average Price)：成交量加權平均價格

• 表明金融資產的「平均」價格，某個價格的成交量越高，該價格所佔的權重就越大
• VWAP就是以成交量爲權重計算出來的加權平均值

TWAP(Time-Weighted Average Price)：時間加權平均價格

• 只是一個變種，基本的思想就是最近的價格重要性大一些,僅僅爲了說明問題，並不必定很是正確
• 選用arange函數建立一個從0開始依次增加的天然數序列，天然數的個數即爲收盤價的個數

例題：eg3_vwap_twap.py

利用data.csv文件，讀出收盤價、成交量，計算成交量加權平均價格VWAP、算術平均值和時間加權平均價格TWAP，程序保存爲eg13_vwap_twap3.py

import numpy as np
c,v=np.loadtxt('data.csv', delimiter=',', \
               usecols=(6,7), unpack=True,skiprows=1)   #usecols列數是從0開始的
print('收盤價：',c)
print('成交量：',v)
vwap=np.average(c,weights=v)
print('成交量加權平均價格VWAP=%f'%vwap)  #這邊看不太懂？成交量加權平均價格
print('算術平均值mean1=',np.mean(c)) #算術平均值
print('算術平均值mean1=',c.mean())   #算術平均值
t=np.arange(len(c))
print('時間加權平均價格=',np.average(c,weights=t))  #時間加權平均價格

收盤價： [336.1 339.32 345.03 344.32 343.44 346.5 351.88 355.2 358.16 354.54
356.85 359.18 359.9 363.13 358.3 350.56 338.61 342.62 342.88 348.16
353.21 349.31 352.12 359.56 360. 355.36 355.76 352.47 346.67 351.99]
成交量： [21144800. 13473000. 15236800. 9242600. 14064100. 11494200. 17322100.

1. 1. 1. 1. 1. 1. 17184100.
2. 1. 1. 1. 1. 1. 14395400.
3. 1. 1. 1. 1. 1. 16192700.
4. 16824200.]

成交量加權平均價格VWAP=350.589549
算術平均值mean1= 351.0376666666667
算術平均值mean1= 351.0376666666667
時間加權平均價格= 352.4283218390804

ptp函數

• ptp函數能夠計算數組的取值範圍，返回數組元素的最大值和最小值之間的差值，即返回值=max(array)-min(array)

示例4：eg13_max_min4.py

利用data.csv，計算最高價中的最大值和最小值以及最高價和最低價中最大值和最小值之間的差值，程序保存爲eg13_max_min4.py

import numpy as np
hst,lst=np.loadtxt('data.csv',delimiter=',',usecols=(4,5),\
                   unpack=True,skiprows=1)
hh=np.max(hst)
ll=np.min(hst)
print('最高價中最大值=',hh)
print('最高價中最小值=',ll)  #np.ptp()
print('最高價中最大值和最小值之間的差值=','%.2f'%np.ptp(hst))    
print('最低價中最大值和最小值之間的差值=','%.2f'%np.ptp(lst))

最高價中最大值= 364.9
最高價中最小值= 340.04
最高價中最大值和最小值之間的差值= 24.86
最低價中最大值和最小值之間的差值= 26.97

數組的基本統計分析函數及說明

argmax、argmin、max、min

import numpy as np
b=np.array([1,5,8,9,334])
print(np.argmax(b))  #np.argmax();np.argmin
print(np.argmin(b))    #reshape()的乘積恰好是arange()的列*行
print(np.arange(7,17).reshape(2,5))

4
0
[[ 7 8 9 10 11]
[12 13 14 15 16]]

axis

axis=0 跨行
axis=1 跨列

import numpy as np
a=np.arange(7,17).reshape(2,5)
print(a)
print(np.argmin(a,axis=0))
print(np.argmin(a,axis=1))

[[ 7 8 9 10 11]
[12 13 14 15 16]]
[0 0 0 0 0]
[0 0]

NumPy經常使用函數：計算中位數和方差

思考1：si13_var1.py

• 中位數：np.median()
• 方差：np.var()

import numpy as np
price=np.loadtxt('data.csv',delimiter=',',usecols=(6,),unpack=True,skiprows=1)
print(price)
print('中位數：',np.median(price))
print('方差：', np.var(price))

[336.1 339.32 345.03 344.32 343.44 346.5 351.88 355.2 358.16 354.54
356.85 359.18 359.9 363.13 358.3 350.56 338.61 342.62 342.88 348.16
353.21 349.31 352.12 359.56 360. 355.36 355.76 352.47 346.67 351.99]
中位數： 352.055
方差： 50.126517888888884

import numpy as np
price=np.loadtxt('data.csv',delimiter=',',usecols=(6,),unpack=True,skiprows=1)
print(price)
print('中位數',np.median(price))
#驗證剛纔求的中位數是否正確
#一、將數組進行排序
sorted = np.msort(price)
print(sorted)
#二、計算數組的元素個數
n = len(sorted)
print(n)
#print('middle', '=', sorted[n - 1]//2)  #n爲奇數
print('中位數：', (sorted[n // 2] + sorted[(n - 1) // 2]) / 2)    #n爲偶數
#方差等於各個數據與平均數之差的平方和的平均數,用來度量隨機變量和其數學指望（即均值）之間的偏離程度
#((x1 - a)^2 + (x2 - a)^2 + (x3 - a)^2 ... + (xn - a)^2) / n
print('方差：', np.var(price))

[336.1 339.32 345.03 344.32 343.44 346.5 351.88 355.2 358.16 354.54
356.85 359.18 359.9 363.13 358.3 350.56 338.61 342.62 342.88 348.16
353.21 349.31 352.12 359.56 360. 355.36 355.76 352.47 346.67 351.99]
中位數 352.055
[336.1 338.61 339.32 342.62 342.88 343.44 344.32 345.03 346.5 346.67
348.16 349.31 350.56 351.88 351.99 352.12 352.47 353.21 354.54 355.2
355.36 355.76 356.85 358.16 358.3 359.18 359.56 359.9 360. 363.13]
30
中位數： 352.055
方差： 50.126517888888884

簡單收益率和對數收益率

• 簡單收益率：相鄰兩個價格之間的變化率
• 對數收益率：全部價格取對數後兩兩之間的差值，也能夠用來衡量價格的變化率

計算曆史波動率（如年波動率和月波動率）時，須要用到對數收益率。

* 若是a的x次方等於N（a>0，且a不等於1），那麼數x叫作以a爲底N的對數（logarithm），記做x=logaN。其中，a叫作對數的底數，N叫作真數。

• 年波動率 = 對數收益率的標準差std / 其均值mean，再除以交易日倒數的平方根。一般交易日取252天。

diff函數

• 返回一個由相鄰數組元素的差值構成的數組

import numpy as np
a=np.array([1,2,5,4,3,7,8,38])
print(np.diff(a))

[ 1 3 -1 -1 4 1 30]

numpy.where（暫時不看）

where(condition, [x, y])：

-- x, y不爲空: condition爲True, 返回x; False, 返回y
-- x, y爲空： 返回condition爲True的數組下標

np.where(關係表達式)：數組中知足關係表達式的元素的下標數組

numpy.take（暫時不看）

eg13_diff5.py

np.take(數組，下標數組)：數組中由下標數組所表示的元素集合
利用data.csv，對於收盤價，計算簡單收益率及標準差、對數收益率及標準差、年波動率和月波動率

import numpy as np
c=np.loadtxt('data.csv',delimiter=',',usecols=(6,),\
             unpack=True,skiprows=1)
returns=np.diff(c)/c[:-1]
print('簡單收益率',returns)
print('簡單收益率的標準差:',np.std(returns))
logreturns=np.diff(np.log(c))
print('對數收益率：',logreturns)
print('對數收益率的標準差：',np.std(logreturns))
posretindices=np.where(returns>0)
print('全部正值元素的元素值:',np.take(returns.posretindices)) #有問題
annual_volatility=np.std(logreturns)/np.mean(logreturns)
annual_volatility=annual_volatility/np.sqrt(1/252)
print('年波動率',annual_volatility)
print('月波動率',annual_volatility*np.sqrt(1/12))