Simple Linear Regression: R vs. Python

有時候聽到一些人在講 Python 比 R 「更適合」作資料分析或科學計算,我就覺得他們在胡說八道。

買 iPhone 人通常都被別人認為比較有品味,但如果把 iPhone 拿來當鐵鎚,固然能把釘子敲進去,但終究無法證明 iPhone 是適合當鐵鎚,更無法證明 iPhone 比鐵鎚更像鐵鎚。要敲釘子的話,還是拿鐵鎚來當工具比較適當,也比較直接容易。

在這裡,我們不使用 R 跟 Python 的 Regression package/module,純粹比較兩種程式語言「內建(built-in)」的數學與向量計算功能、從理論撰寫出程式的直覺跟簡潔、與程式碼的多寡,計算簡單線性迴歸 (Simple Linear Regression)

後記:

有專精 Python 的朋友向我抗議,認為我這篇文章用的比較程式不正確的例子,他並寫了 只有 10 行的 Python 程式(外部連結) 來對照。

我的看法是,如果純粹只是要比誰程式短的話,我還可以更精簡成只要 2 行:

d = read.csv(“d:/data/insurance.csv”)
c(beta0=nrow(d)*sum(d$x*d$y) – sum(d$x)*sum(d$y))/(n*sum(d$x^2)-sum(d$x)^2), beta1=mean(d$y)-beta1*mean(d$x))

我是覺得我舉的例子已經很正確。因為底下的 R 程式是只要學過一學期中前一兩週的 R 課程學生都可以寫得出這樣的程式碼。但 Python 朋友寫的程式卻不是一般初學者可以寫得出來的。一方派出輕量級,另一方派出重量級,結果後者程式還是比較長,這樣好向並沒有比較「正確」:)

況且,如果我們再從簡單線性迴歸問題擴充成複迴歸(Multiple Regression),並把解釋變數數目改成 1000 個,轉為矩陣求解的話,則 R 更可以輕鬆地在 3 行之內不用任何外掛套件、不用 lm 函數就可以搞定,那上面的 Python 要改成幾行呢?

d = read.csv(“DataWith1001Variables.csv”)
X=cbind(1,as.matrix(subset(d,select=-y))) # 扣除 y 變數
# 或 X=cbind(1,do.call(cbind,subset(d,select=-y))) 
( betas = solve(t(X)%*%X)%*%t(X)%*%d$y )

此外,Python 朋友改寫後的 Python 程式,還是彰顯了 Python 並沒有內建 mean 函數、向量、與向量化運算的窘態。光看這一點,要宣稱 Python 比 R 更適合做資料分析是很勉強的。

如果換個說法「Python 『目前』比 R 適合做深度學習」,那我完全同意。不過,一來,「深度學習」從來就不是「資料分析」的全部。二來,Python 在深度學習過程只是擔任一個仲介者的角色,後端的運算 Libraries 根本不是用 Python 寫的。如果只是中間人角色,那 R 軟體上面遲早也都會有這類的介接 packages 出現。

 

(1) R : 7 lines

R 語言內建 CSV 檔讀寫功能、平均數與許多統計機率相關函數、與許多向量化運算函數 (Vectorized Functions)

?View Code LANGUAGE
xdata = read.csv("d:/data/insurance.csv")   # x, y
attach(xdata)
n = length(x)
beta1 = (n*sum(x*y) - sum(x)*sum(y))/(n*sum(x^2)-sum(x)^2)
beta0 = mean(y)-beta1*mean(x)
options(digits=12)
cat("beta0 = ",beta0,",beta1 = ",beta1,"\n")

Output: beta0 = 19.9944857591 ,beta1 = 3.41382356007

 

(2) Python : 36 lines

單純不加特殊外掛的 Python 並沒有內建 CSV 檔讀寫功能,沒有向量化運算函數,甚至連平均數 (mean) 都得自己寫函數去計算

(modified from Python code in  How To Implement Simple Linear Regression From Scratch With Python )

?View Code LANGUAGE
from csv import reader
from math import sqrt
def load_csv(filename):
	dataset = list()
	with open(filename, 'r') as file:
		csv_reader = reader(file)
		for row in csv_reader:
			if not row:
				continue
			dataset.append(row)
	return dataset
def str_column_to_float(dataset, column):
	for row in dataset:
		row[column] = float(row[column].strip())
def mean(values):
	return sum(values) / float(len(values))
def covariance(x, mean_x, y, mean_y):
	covar = 0.0
	for i in range(len(x)):
		covar += (x[i] - mean_x) * (y[i] - mean_y)
	return covar
def variance(values, mean):
	return sum([(x-mean)**2 for x in values])
def coefficients(dataset):
	x = [row[0] for row in dataset]
	y = [row[1] for row in dataset]
	x_mean, y_mean = mean(x), mean(y)
	b1 = covariance(x, x_mean, y, y_mean) / variance(x, x_mean)
	b0 = y_mean - b1 * x_mean
	return [b0, b1]
filename = 'insurance.csv'
dataset = load_csv(filename)
for i in range(len(dataset[0])):
	str_column_to_float(dataset, i)	
betas = coefficients(dataset)
print "beta0 = ",betas[0],",beta1 = ",betas[1]

Output: beta0 = 19.9944857591 ,beta1 = 3.41382356007