使用Python和GNU Octave繪制數(shù)據(jù)
了解如何使用 Python 和 GNU Octave 完成一項(xiàng)常見(jiàn)的數(shù)據(jù)科學(xué)任務(wù)。
數(shù)據(jù)科學(xué)是跨越編程語(yǔ)言的知識(shí)領(lǐng)域。有些語(yǔ)言以解決這一領(lǐng)域的問(wèn)題而聞名，而另一些則鮮為人知。這篇文章將幫助你熟悉用一些流行的語(yǔ)言完成數(shù)據(jù)科學(xué)的工作。
選擇 Python 和 GNU Octave 做數(shù)據(jù)科學(xué)工作
我經(jīng)常嘗試學(xué)習(xí)一種新的編程語(yǔ)言。為什么？這既有對(duì)舊方式的厭倦，也有對(duì)新方式的好奇。當(dāng)我開(kāi)始學(xué)習(xí)編程時(shí)，我唯一知道的語(yǔ)言是 C 語(yǔ)言。那些年的編程生涯既艱難又危險(xiǎn)，因?yàn)槲冶仨毷謩?dòng)分配內(nèi)存、管理指針、并記得釋放內(nèi)存。
后來(lái)一個(gè)朋友建議我試試 Python，現(xiàn)在我的編程生活變得輕松多了。雖然程序運(yùn)行變得慢多了，但我不必通過(guò)編寫(xiě)分析軟件來(lái)受苦了。然而，我很快就意識(shí)到每種語(yǔ)言都有比其它語(yǔ)言更適合自己的應(yīng)用場(chǎng)景。后來(lái)我學(xué)習(xí)了一些其它語(yǔ)言，每種語(yǔ)言都給我?guī)?lái)了一些新的啟發(fā)。發(fā)現(xiàn)新的編程風(fēng)格讓我可以將一些解決方案移植到其他語(yǔ)言中，這樣一切都變得有趣多了。
為了對(duì)一種新的編程語(yǔ)言（及其文檔）有所了解，我總是從編寫(xiě)一些執(zhí)行我熟悉的任務(wù)的示例程序開(kāi)始。為此，我將解釋如何用 Python 和 GNU Octave 編寫(xiě)一個(gè)程序來(lái)完成一個(gè)你可以歸類為數(shù)據(jù)科學(xué)的特殊任務(wù)。如果你已經(jīng)熟悉其中一種語(yǔ)言，從它開(kāi)始，然后通過(guò)其他語(yǔ)言尋找相似之處和不同之處。這篇文章并不是對(duì)編程語(yǔ)言的詳盡比較，只是一個(gè)小小的展示。
所有的程序都應(yīng)該在命令行上運(yùn)行，而不是用圖形用戶界面（GUI）。完整的例子可以在 polyglot_fit 存儲(chǔ)庫(kù)中找到。
編程任務(wù)
你將在本系列中編寫(xiě)的程序:
- 從 CSV 文件中讀取數(shù)據(jù)
- 用直線插入數(shù)據(jù)（例如
f(x)=m ⋅ x + q） - 將結(jié)果生成圖像文件
這是許多數(shù)據(jù)科學(xué)家遇到的常見(jiàn)情況。示例數(shù)據(jù)是 Anscombe 的四重奏的第一組，如下表所示。這是一組人工構(gòu)建的數(shù)據(jù)，當(dāng)用直線擬合時(shí)會(huì)給出相同的結(jié)果，但是它們的曲線非常不同。數(shù)據(jù)文件是一個(gè)文本文件，以制表符作為列分隔符，開(kāi)頭幾行作為標(biāo)題。此任務(wù)將僅使用第一組（即前兩列）。
Python 方式
Python 是一種通用編程語(yǔ)言，是當(dāng)今最流行的語(yǔ)言之一（依據(jù) TIOBE 指數(shù)、RedMonk 編程語(yǔ)言排名、編程語(yǔ)言流行指數(shù)、GitHub Octoverse 狀態(tài)和其他來(lái)源的調(diào)查結(jié)果）。它是一種解釋型語(yǔ)言；因此，源代碼由執(zhí)行該指令的程序讀取和評(píng)估。它有一個(gè)全面的標(biāo)準(zhǔn)庫(kù)并且總體上非常好用（我對(duì)這最后一句話沒(méi)有證據(jù)；這只是我的拙見(jiàn)）。
安裝
要使用 Python 開(kāi)發(fā)，你需要解釋器和一些庫(kù)。最低要求是：
- NumPy 用于簡(jiǎn)化數(shù)組和矩陣的操作
- SciPy 用于數(shù)據(jù)科學(xué)
- Matplotlib 用于繪圖
在 Fedora 安裝它們是很容易的：
sudo dnf install python3 python3-numpy python3-scipy python3-matplotlib
代碼注釋
在 Python中，注釋是通過(guò)在行首添加一個(gè) # 來(lái)實(shí)現(xiàn)的，該行的其余部分將被解釋器丟棄：
# 這是被解釋器忽略的注釋。
fitting_python.py 示例使用注釋在源代碼中插入許可證信息，第一行是特殊注釋，它允許該腳本在命令行上執(zhí)行:
#!/usr/bin/env python3
這一行通知命令行解釋器，該腳本需要由程序 python3 執(zhí)行。
需要的庫(kù)
在 Python 中，庫(kù)和模塊可以作為一個(gè)對(duì)象導(dǎo)入（如示例中的第一行），其中包含庫(kù)的所有函數(shù)和成員。可以通過(guò)使用 as 方式用自定義標(biāo)簽重命名它們：
import numpy as npfrom scipy import statsimport matplotlib.pyplot as plt
你也可以決定只導(dǎo)入一個(gè)子模塊（如第二行和第三行）。語(yǔ)法有兩個(gè)（基本上）等效的方式：import module.submodule 和 from module import submodule。
定義變量
Python 的變量是在第一次賦值時(shí)被聲明的：
input_file_name = "anscombe.csv"delimiter = "\t"skip_header = 3column_x = 0column_y = 1
變量類型由分配給變量的值推斷。沒(méi)有具有常量值的變量，除非它們?cè)谀K中聲明并且只能被讀取。習(xí)慣上，不應(yīng)被修改的變量應(yīng)該用大寫(xiě)字母命名。
打印輸出
通過(guò)命令行運(yùn)行程序意味著輸出只能打印在終端上。Python 有 print() 函數(shù)，默認(rèn)情況下，該函數(shù)打印其參數(shù)，并在輸出的末尾添加一個(gè)換行符：
print("#### Anscombe's first set with Python ####")
在 Python 中，可以將 print() 函數(shù)與字符串類的格式化能力相結(jié)合。字符串具有format 方法，可用于向字符串本身添加一些格式化文本。例如，可以添加格式化的浮點(diǎn)數(shù)，例如:
print("Slope: {:f}".format(slope))
讀取數(shù)據(jù)
使用 NumPy 和函數(shù) genfromtxt() 讀取 CSV 文件非常容易，該函數(shù)生成 NumPy 數(shù)組：
data = np.genfromtxt(input_file_name, delimiter = delimiter, skip_header = skip_header)
在 Python 中，一個(gè)函數(shù)可以有數(shù)量可變的參數(shù)，你可以通過(guò)指定所需的參數(shù)來(lái)傳遞一個(gè)參數(shù)的子集。數(shù)組是非常強(qiáng)大的矩陣狀對(duì)象，可以很容易地分割成更小的數(shù)組：
x = data[:, column_x]y = data[:, column_y]
冒號(hào)選擇整個(gè)范圍，也可以用來(lái)選擇子范圍。例如，要選擇數(shù)組的前兩行，可以使用：
first_two_rows = data[0:1, :]
擬合數(shù)據(jù)
SciPy 提供了方便的數(shù)據(jù)擬合功能，例如 linregress() 功能。該函數(shù)提供了一些與擬合相關(guān)的重要值，如斜率、截距和兩個(gè)數(shù)據(jù)集的相關(guān)系數(shù):
slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)print("Slope: {:f}".format(slope))print("Intercept: {:f}".format(intercept))print("Correlation coefficient: {:f}".format(r_value))
因?yàn)?linregress() 提供了幾條信息，所以結(jié)果可以同時(shí)保存到幾個(gè)變量中。
繪圖
Matplotlib 庫(kù)僅僅繪制數(shù)據(jù)點(diǎn)，因此，你應(yīng)該定義要繪制的點(diǎn)的坐標(biāo)。已經(jīng)定義了 x 和 y 數(shù)組，所以你可以直接繪制它們，但是你還需要代表直線的數(shù)據(jù)點(diǎn)。
fit_x = np.linspace(x.min() - 1, x.max() + 1, 100)
linspace() 函數(shù)可以方便地在兩個(gè)值之間生成一組等距值。利用強(qiáng)大的 NumPy 數(shù)組可以輕松計(jì)算縱坐標(biāo)，該數(shù)組可以像普通數(shù)值變量一樣在公式中使用：
fit_y = slope * fit_x + intercept
該公式在數(shù)組中逐元素應(yīng)用；因此，結(jié)果在初始數(shù)組中具有相同數(shù)量的條目。
要繪圖，首先，定義一個(gè)包含所有圖形的圖形對(duì)象：
fig_width = 7 #inchfig_height = fig_width / 16 * 9 #inchfig_dpi = 100fig = plt.figure(figsize = (fig_width, fig_height), dpi = fig_dpi)
一個(gè)圖形可以畫(huà)幾個(gè)圖；在 Matplotlib 中，這些圖被稱為軸。本示例定義一個(gè)單軸對(duì)象來(lái)繪制數(shù)據(jù)點(diǎn)：
ax = fig.add_subplot(111)ax.plot(fit_x, fit_y, label = "Fit", linestyle = '-')ax.plot(x, y, label = "Data", marker = '.', linestyle = '')ax.legend()ax.set_xlim(min(x) - 1, max(x) + 1)ax.set_ylim(min(y) - 1, max(y) + 1)ax.set_xlabel('x')ax.set_ylabel('y')
將該圖保存到 PNG 圖形文件中，有:
fig.savefig('fit_python.png')
如果要顯示（而不是保存）該繪圖，請(qǐng)調(diào)用：
plt.show()
此示例引用了繪圖部分中使用的所有對(duì)象：它定義了對(duì)象 fig 和對(duì)象 ax。這在技術(shù)上是不必要的，因?yàn)?plt 對(duì)象可以直接用于繪制數(shù)據(jù)集。《Matplotlib 教程》展示了這樣一個(gè)接口：
plt.plot(fit_x, fit_y)
坦率地說(shuō)，我不喜歡這種方法，因?yàn)樗[藏了各種對(duì)象之間發(fā)生的重要交互。不幸的是，有時(shí)官方的例子有點(diǎn)令人困惑，因?yàn)樗麄儍A向于使用不同的方法。在這個(gè)簡(jiǎn)單的例子中，引用圖形對(duì)象是不必要的，但是在更復(fù)雜的例子中（例如在圖形用戶界面中嵌入圖形時(shí)），引用圖形對(duì)象就變得很重要了。
結(jié)果
命令行輸入：
#### Anscombe's first set with Python ####Slope: 0.500091Intercept: 3.000091Correlation coefficient: 0.816421
這是 Matplotlib 產(chǎn)生的圖像：
Plot and fit of the dataset obtained with Python
GNU Octave 方式
GNU Octave 語(yǔ)言主要用于數(shù)值計(jì)算。它提供了一個(gè)簡(jiǎn)單的操作向量和矩陣的語(yǔ)法，并且有一些強(qiáng)大的繪圖工具。這是一種像 Python 一樣的解釋語(yǔ)言。由于 Octave 的語(yǔ)法幾乎兼容 MATLAB，它經(jīng)常被描述為一個(gè)替代 MATLAB 的免費(fèi)方案。Octave 沒(méi)有被列為最流行的編程語(yǔ)言，而 MATLAB 則是，所以 Octave 在某種意義上是相當(dāng)流行的。MATLAB 早于 NumPy，我覺(jué)得它是受到了前者的啟發(fā)。當(dāng)你看這個(gè)例子時(shí)，你會(huì)看到相似之處。
安裝
fitting_octave.m 的例子只需要基本的 Octave 包，在 Fedora 中安裝相當(dāng)簡(jiǎn)單：
sudo dnf install octave
代碼注釋
在 Octave 中，你可以用百分比符號(hào)（%）為代碼添加注釋，如果不需要與 MATLAB 兼容，你也可以使用 #。使用 # 的選項(xiàng)允許你編寫(xiě)像 Python 示例一樣的特殊注釋行，以便直接在命令行上執(zhí)行腳本。
必要的庫(kù)
本例中使用的所有內(nèi)容都包含在基本包中，因此你不需要加載任何新的庫(kù)。如果你需要一個(gè)庫(kù)，語(yǔ)法是 pkg load module。該命令將模塊的功能添加到可用功能列表中。在這方面，Python 具有更大的靈活性。
定義變量
變量的定義與 Python 的語(yǔ)法基本相同：
input_file_name = "anscombe.csv";delimiter = "\t";skip_header = 3;column_x = 1;column_y = 2;
請(qǐng)注意，行尾有一個(gè)分號(hào)；這不是必需的，但是它會(huì)抑制該行結(jié)果的輸出。如果沒(méi)有分號(hào)，解釋器將打印表達(dá)式的結(jié)果：
octave:1> input_file_name = "anscombe.csv"input_file_name = anscombe.csvoctave:2> sqrt(2)ans = 1.4142
打印輸出結(jié)果
強(qiáng)大的函數(shù) printf() 是用來(lái)在終端上打印的。與 Python 不同，printf() 函數(shù)不會(huì)自動(dòng)在打印字符串的末尾添加換行，因此你必須添加它。第一個(gè)參數(shù)是一個(gè)字符串，可以包含要傳遞給函數(shù)的其他參數(shù)的格式信息，例如：
printf("Slope: %f\n", slope);
在 Python 中，格式是內(nèi)置在字符串本身中的，但是在 Octave 中，它是特定于 printf() 函數(shù)。
讀取數(shù)據(jù)
dlmread() 函數(shù)可以讀取類似 CSV 文件的文本內(nèi)容：
data = dlmread(input_file_name, delimiter, skip_header, 0);
結(jié)果是一個(gè)矩陣對(duì)象，這是 Octave 中的基本數(shù)據(jù)類型之一。矩陣可以用類似于 Python 的語(yǔ)法進(jìn)行切片：
x = data(:, column_x);y = data(:, column_y);
根本的區(qū)別是索引從 1 開(kāi)始，而不是從 0 開(kāi)始。因此，在該示例中，x 列是第一列。
擬合數(shù)據(jù)
要用直線擬合數(shù)據(jù)，可以使用 polyfit() 函數(shù)。它用一個(gè)多項(xiàng)式擬合輸入數(shù)據(jù)，所以你只需要使用一階多項(xiàng)式：
p = polyfit(x, y, 1);slope = p(1);intercept = p(2);
結(jié)果是具有多項(xiàng)式系數(shù)的矩陣；因此，它選擇前兩個(gè)索引。要確定相關(guān)系數(shù)，請(qǐng)使用 corr() 函數(shù)：
r_value = corr(x, y);
最后，使用 printf() 函數(shù)打印結(jié)果：
printf("Slope: %f\n", slope);printf("Intercept: %f\n", intercept);printf("Correlation coefficient: %f\n", r_value);
繪圖
與 Matplotlib 示例一樣，首先需要?jiǎng)?chuàng)建一個(gè)表示擬合直線的數(shù)據(jù)集:
fit_x = linspace(min(x) - 1, max(x) + 1, 100);fit_y = slope * fit_x + intercept;
與 NumPy 的相似性也很明顯，因?yàn)樗褂昧?linspace() 函數(shù)，其行為就像 Python 的等效版本一樣。
同樣，與 Matplotlib 一樣，首先創(chuàng)建一個(gè)圖對(duì)象，然后創(chuàng)建一個(gè)軸對(duì)象來(lái)保存這些圖：
fig_width = 7; %inchfig_height = fig_width / 16 * 9; %inchfig_dpi = 100;fig = figure("units", "inches","position", [1, 1, fig_width, fig_height]);ax = axes("parent", fig);set(ax, "fontsize", 14);set(ax, "linewidth", 2);
要設(shè)置軸對(duì)象的屬性，請(qǐng)使用 set() 函數(shù)。然而，該接口相當(dāng)混亂，因?yàn)樵摵瘮?shù)需要一個(gè)逗號(hào)分隔的屬性和值對(duì)列表。這些對(duì)只是代表屬性名的一個(gè)字符串和代表該屬性值的第二個(gè)對(duì)象的連續(xù)。還有其他設(shè)置各種屬性的函數(shù)：
xlim(ax, [min(x) - 1, max(x) + 1]);ylim(ax, [min(y) - 1, max(y) + 1]);xlabel(ax, 'x');ylabel(ax, 'y');
繪圖是用 plot() 功能實(shí)現(xiàn)的。默認(rèn)行為是每次調(diào)用都會(huì)重置坐標(biāo)軸，因此需要使用函數(shù) hold()。
hold(ax, "on");plot(ax, fit_x, fit_y,"marker", "none","linestyle", "-","linewidth", 2);plot(ax, x, y,"marker", ".","markersize", 20,"linestyle", "none");hold(ax, "off");
此外，還可以在 plot() 函數(shù)中添加屬性和值對(duì)。legend 必須單獨(dú)創(chuàng)建，標(biāo)簽應(yīng)手動(dòng)聲明：
lg = legend(ax, "Fit", "Data");set(lg, "location", "northwest");
最后，將輸出保存到 PNG 圖像：
image_size = sprintf("-S%f,%f", fig_width * fig_dpi, fig_height * fig_dpi);image_resolution = sprintf("-r%f,%f", fig_dpi);print(fig, 'fit_octave.png','-dpng',image_size,image_resolution);
令人困惑的是，在這種情況下，選項(xiàng)被作為一個(gè)字符串傳遞，帶有屬性名和值。因?yàn)樵?Octave 字符串中沒(méi)有 Python 的格式化工具，所以必須使用 sprintf() 函數(shù)。它的行為就像 printf() 函數(shù)，但是它的結(jié)果不是打印出來(lái)的，而是作為字符串返回的。
在這個(gè)例子中，就像在 Python 中一樣，圖形對(duì)象很明顯被引用以保持它們之間的交互。如果說(shuō) Python 在這方面的文檔有點(diǎn)混亂，那么 Octave 的文檔就更糟糕了。我發(fā)現(xiàn)的大多數(shù)例子都不關(guān)心引用對(duì)象；相反，它們依賴于繪圖命令作用于當(dāng)前活動(dòng)圖形。全局根圖形對(duì)象跟蹤現(xiàn)有的圖形和軸。
結(jié)果
命令行上的結(jié)果輸出是：
#### Anscombe's first set with Octave ####Slope: 0.500091Intercept: 3.000091Correlation coefficient: 0.816421
它顯示了用 Octave 生成的結(jié)果圖像。
Plot and fit of the dataset obtained with Octave
接下來(lái)
Python 和 GNU Octave 都可以繪制出相同的信息，盡管它們的實(shí)現(xiàn)方式不同。如果你想探索其他語(yǔ)言來(lái)完成類似的任務(wù)，我強(qiáng)烈建議你看看 Rosetta Code。這是一個(gè)了不起的資源，可以看到如何用多種語(yǔ)言解決同樣的問(wèn)題。
