R和Python利用蒙特卡洛模拟求π——向量化编程，以及apply()运行时间函数探讨

MC的最根本思想就是一句话：通过生成随机数来模拟复杂的积分运算。蒙特卡洛模拟的意义在于简化高维积分运算，降低涉及高维积分的算法的时间复杂性。

1.编写代码
本次例子我们利用蒙特卡洛模拟来计算π，正方形内切的圆与该正方形比值为π/4

R代码如下

#圆心坐标
center <- c(0,0)
#该函数是用来计算，各点到圆心之间的距离
dist_from_center <- function(a){
  sqrt(sum((a-center)^2))
}
#随机种子
set.seed(2)
n <- 10000
#开始时间
start_time = proc.time()
#[0,1)均匀分布生成n行，2列的矩阵，按行填充矩阵
mx <- matrix(data = runif(2*n,min = 0,max = 1),
             nrow = n,ncol = 2,byrow = TRUE)
#利用apply函数，1表示对行进行计算，返回的distance是向量
distance <- apply(mx,1,dist_from_center)
#distance<1的向量只有TURE和FALSE两个值，
#TRUE值就是1，FALSE值就是0.对该向量求均值就是相当于
#求距离小于1的点所占的比例
number <- mean(distance<1)
#求出π
pi <- 4*number
#运行时间
t <- proc.time() - start_time
print(paste0('π is ',pi))
print(paste0('执行时间：',t[3],'秒'))

R运行结果

[1] "π is 3.1476"
[1] "执行时间：0.0299999999999727秒"

Python代码

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from time import perf_counter


np.random.seed(2)
n = 1000000
#圆心坐标
circle_center = (0,0)
start = perf_counter()
#以[0,1)均匀分布生成n行，2列的数据框
df = pd.DataFrame(np.random.uniform(low=0.0,high=1.0,size=(n,2)))
#该函数是用来计算，各点到圆心之间的距离
distance_from_center = lambda x:np.sqrt(((x)**2).sum())
#利用apply函数，1表示对行进行计算，返回的distance是Series
distance = df.apply(distance_from_center,axis=1)
#distance<1的向量只有TURE和FALSE两个值，
#TRUE值就是1，FALSE值就是0.对该向量求均值就是相当于
#求距离小于1的点所占的比例
number = np.mean(distance<1)
#求出π
pi = number*4

runningtime = perf_counter()-start
print(f'π is {pi}')
print(f'执行时间：{runningtime}秒')

Python运行结果

π is 3.1456
执行时间：1.885769400000072秒

结果分析
R语言就是向量化编程的语言，基本的运算单位就是向量，直接对向量进行算术运算，对数据分析特别方便。
Python编程的基本单位是标量，但是Python也可以向量化编程，需要借助第三方库，numpy，pandas，数据展示需要用到matplotlib.pyplot，语法的书写和R语言特别相似
结果中可以看到，R语言和Python中基本相同的代码来模拟蒙特卡洛来算π，R语言只用了大约0.03秒，而Python竟然用了1.88秒，足足相差了62倍！
我将模拟的点数调整为1000,000个，再进行模拟，时间就更加直观了
1000,000个点的R运行结果

[1] "π is 3.140352"
[1] "执行时间：2.72秒"

1000,000个点的Python运行结果

π is 3.14156
执行时间：170.67063510000025秒

当数量到达1000,000个点时，Python竟然运行了170秒，这也太慢了吧，难道面对数量非常大的数据框，Python这么不堪一击？
进一步研究发现Python和R都是在apply()那一行代码花费了绝大部分时间，因为apply函数将其他函数应用到矩阵的每一行，或数据框的每一行，Python第三方库pandas的apply函数返回的是serie或者dataframe，当应用于非常大的数据框，有很大的IO开销。

我尝试对二者都不使用apply函数，python直接改变sum()的参数，对行求和

np.random.seed(2)
n = 1000000
circle_center = (0,0)
start = perf_counter()

df = pd.DataFrame(np.random.uniform(low=0.0,high=1.0,size=(n,2)))
#此处不使用apply函数，改用sum(axis=1)
distance = np.sqrt(((df-circle_center)**2).sum(axis=1))
number = np.mean(distance<1)
pi = number*4

runningtime = perf_counter()-start
print(f'π is {pi}')
print(f'执行时间：{runningtime}秒')

运行结果

π is 3.14156
执行时间：0.09306700000024648秒

抛弃apply函数后，竟然时间从170秒直接变成0.09秒，这差别也太大了吧

R语言不使用apply函数，而使用内置的rowSums()函数就可以对矩阵的行进行求和

center <- c(0,0)
set.seed(2)
n <- 1000000
start_time = proc.time()
mx <- matrix(data = runif(2*n,min = 0,max = 1),
             nrow = n,ncol = 2,byrow = TRUE)
#此处不使用apply函数,改用rowSums()             
distance <- sqrt(rowSums((mx-center)^2))
number <- mean(distance<1)
pi <- 4*number
t <- proc.time() - start_time
print(paste0('π is ',pi))
print(paste0('执行时间：',t[3],'秒'))

运行结果

[1] "π is 3.140352"
[1] "执行时间：0.139999999999986秒"

R的运行时间竟然也比用apply函数的运行时间缩短了，并且抛弃apply()函数时，Python运行的反而更快了一些

2.个人结论心得

R和Python用apply函数作用于大量数据的矩阵或者数据框时，apply函数运行时间变慢，尤其在Python中用pandas库的apply()函数时间慢到崩溃，所以最好还是不要使用apply()函数了，可能因为apply函数将sum函数应用每一行数据，每一行数据都要调用一次sum函数，调用100万次sum函数会耗费大量时间，所以导致apply函数运行很慢，而直接用sum函数，只会作用于数据数据框一次，所以运行很快，如果我的理解有错误，希望有人可以指出，感激不尽！

3.数据展示

在展示数据方面，Python和R也极为相似，对于简单的图形展示，R语言使用内置函数就能成图，Python需要借助于matplotlib.pyplot

R数据展示代码

par(bg="beige")
#aspect设为1才能使x轴与y轴比例相同，这样画圆才能正确显示，否则变成椭圆
plot(mx,col="azure3",asp = 1)
abline(h=0,col="red",lty="dotdash",lwd=2)
abline(h=1,col="red",lty="dotdash",lwd=2)
abline(v=0,col="red",lty="dotdash",lwd=2)
abline(v=1,col="red",lty="dotdash",lwd=2)


points(mx[distance<1,],col="green")

library(plotrix)
draw.circle(0,0,1,border="coral2",lty="dashed",lwd=2)

图形

Python数据展示代码

fig, ax = plt.subplots(figsize=(15,10),dpi=80,facecolor='beige')

ax.axhline(y=0,ls='--',c='red',lw=2)
ax.axhline(y=1,ls='--',c='red',lw=2)
ax.axvline(x=0,ls='--',c='red',lw=2)
ax.axvline(x=1,ls='--',c='red',lw=2)
circle = plt.Circle(xy=(0,0),radius=1,fill= False,edgecolor='red',ls='--',lw=2)

color = pd.Series(data=range(n))
color[distance<1] = 'green'
color[distance>=1] = 'blue'
ax.scatter(x=df[0], y=df[1],c=color)
#aspect设为1才能使x轴与y轴比例相同，这样画圆才能正确显示，否则变成椭圆
ax.set_aspect(1)
ax.add_artist(circle)
plt.show()

图形