r语言的程序用python实现_R语言和Python跑模拟，哪个快？-CSDN博客

今天在复习线性回归模型的时候，做到了一道题问回归平方和MSR的期望值的题，一时做不出来，遂去R语言写了一个程序，循环一万次，随机生成X与Y，求MSR。然后求一万个MSR的样本均值作为MSR的总体均值的估计。原题如下：

In a small-scale regression study, five observations on

were obtained corresponding to

and

. Assume that

. What are the expected values of MSR and MSE here?

MSE的均值是

，然而MSR的均值是多少呢？我的程序伪代码大致如下：

循环{

生成5个N(0, 0.36)正态随机数作为误差项error

X = [1, 4, 10, 11, 14]

Y=5+3X+error

Y对X回归

获取ANOVA表

从ANOVA表中提取MSR，并且保存在列表中

}

R语言代码如下：

msrs

for(i in 1:10000){

error

reg

ssr

msrs[i]

}

msr_means

for(i in 1:10000){

meani

msr_means[i]

}

plot(msr_means, type='l')

print(t2-t1)

}

mc()

这段代码不仅把模拟经历的时间打印了出来，还画了收敛图。输出的模拟时间为：

Time difference of 16.21569 secs

样本均值向着总体均值收敛的收敛图如下：

横坐标：第i个MSR样本，纵坐标：第1到i个样本的均值

好家伙，我跑一万次就浪费我16秒，那如果我跑一亿次呢？我算了一下，要花44小时:-)……我想了想，Python跑一次同样的代码要多久呢？话不多说，用Python写了一段一模一样的代码。需要注意的是，Python中生成随机数、处理数据、进行线性回归、获取ANOVA表这些，依赖于一些第三方库numpy、pandas、statsmodels等。不像R，直接一行代码完成一项工作。Python代码如下：

import time

import numpy as np

import pandas as pd

import statsmodels.api as sm

from statsmodels.formula.api import ols

import matplotlib.pyplot as plt

def mc():

t1 = time.time()

msrs = []

for i in range(10000):

error = np.random.normal(scale=0.6, size=5)

X = np.array([1, 4, 10, 11, 14])

Y = 5 + 3 * X + error

df = pd.DataFrame({'X': X, 'Y': Y})

lm = ols('Y~X', data=df).fit()

table = sm.stats.anova_lm(lm)

msr = table.loc['X', 'sum_sq']

msrs.append(msr)

msr_means = []

for i in range(10000):

meani = np.mean(msrs[:i + 1])

msr_means.append(meani)

plt.plot(msr_means)

plt.show()

t2 = time.time()

print(t2 - t1)

mc()

收敛图如下：

1000次模拟-Python

天呐，Python循环了一万次居然花了我122秒！这比R慢多了。如果我想要跑一亿次，那么跑这个程序就要花整整14天(╯‵□′)╯︵┻━┻。

太不像样了，我思考了一下，很有可能是在计算OLS估计和获取ANOVA表的过程中花费了很多时间，statsmodels库速度可能太慢了。我自编了一个函数实现了线性回归OLS估计以及求MSR，其中OLS估计利用矩阵运算：

注意本例子中只有一个预测变量，所以MSR和SSR相等。该函数如下：

def simple_ols(X, y):

n = len(y)

y = np.matrix(y.reshape(n, 1))

design = np.matrix(np.insert(X.reshape(n, 1), 0, values=1, axis=1))

b = np.linalg.inv(design.T * design) * design.T * y

y_pred = design * b

y_mean = np.mean(y)

ssr = (y_pred - y_mean).T * (y_pred - y_mean)

return float(ssr)

写好了，来跑一下看看要多久吧！

def mc():

t1 = time.time()

msrs = []

for i in range(10000):

error = np.random.normal(scale=0.6, size=5)

X = np.array([1, 4, 10, 11, 14])

Y = 5 + 3 * X + error

msr = simple_ols(X, Y)

msrs.append(msr)

msr_means = []

for i in range(10000):

meani = np.mean(msrs[:i + 1])

msr_means.append(meani)

plt.plot(msr_means)

plt.show()

t2 = time.time()

print(t2 - t1)

mc()

收敛图我就不放了，最后结果是3.69秒。不错不错，这速度进步比原来调包的Python程序和R程序都要快多了。我算了算，如果要搞一亿次，那么就要10小时。效率仍然很低，再想想看看哪里可以优化的。

我的目的在于对MSR的均值进行估计，因此画收敛图不是必要选项。那么计算那一万次样本均值也不必要了。把这段删掉：

def mc():

t1 = time.time()

msrs = []

for i in range(10000):

error = np.random.normal(scale=0.6, size=5)

X = np.array([1, 4, 10, 11, 14])

Y = 5 + 3 * X + error

msr = simple_ols(X, Y)

msrs.append(msr)

t2 = time.time()

print(t2 - t1)

print(np.mean(msrs))

mc()

这段程序跑了1.4501秒。这就意味着，跑一亿次4个小时。进步显著。再想想，看看还能不能忘哪个方向进步~这时我想到了矩阵求逆可能是一个耗时比较久的运算，话不多说，又自编一个函数进行二阶矩阵求逆：

import time

import numpy as np

def inv_2(X):

# 求矩阵X的逆

a, b, c, d = X[0, 0], X[0, 1], X[1, 0], X[1, 1]

det = a * d - b * c

return np.matrix([[d, -b], [-c, a]]) / det

def simple_ols(X, y):

# y对X回归，输出MSR

n = len(y)

y = np.matrix(y.reshape(n, 1))

design = np.matrix(np.insert(X.reshape(n, 1), 0, values=1, axis=1))

b = inv_2(design.T * design) * design.T * y

y_pred = design * b

y_mean = np.mean(y)

ssr = (y_pred - y_mean).T * (y_pred - y_mean)

return float(ssr)

def mc():

# MC模拟一万次

t1 = time.time()

msrs = []

for i in range(10000):

error = np.random.normal(scale=0.6, size=5)

X = np.array([1, 4, 10, 11, 14])

Y = 5 + 3 * X + error

msr = simple_ols(X, Y)

msrs.append(msr)

t2 = time.time()

print(t2 - t1)

print(np.mean(msrs))

mc()

最后结果是1.5628秒，好吧这还变慢了一丢丢。我没法了。不会DSA，不会C++的编程渣渣已经黔驴技穷了。我的电脑有4个CPU，其实我还可以写个四核并行运算，再把运算时间降为原来的1/4，可是我要复习了，不搞了。