'''
auther:一亩三分地~
function:蒙特卡洛方法实例
'''
import os
import numpy as np
import matplotlib.pyplot as plt
# 蒙特卡洛方法求解定积分1-3范围内,x**3
def cac(max_num,x_low,x_high,func_method):
count = 0
# for i in range(max_num):
# x_range = np.arange(x_low,x_high,0.001) # 用别的随机方法是不是更好????????????
# y_range = np.arange(0,func_method(x_high),0.001)
# # 用随机值代表在矩形面积内投掷
# x_point = np.random.choice(x_range)
# y_point = np.random.choice(y_range)
# # 判断投掷点是否在x**3曲线下面
# if func_method(x_point) >= y_point: # 这一块用np.array对象会不会更快,尤其是含有加法的时候
# count += 1
# 将投掷点用均匀分布形成
x_point = np.random.uniform(x_low, x_high, max_num)
y_point = np.random.uniform(0, func_method(x_high), max_num)
count = sum(func_method(x_point) >= y_point)
# 用投中次数/总的投掷 作为x**3的面积,即K=count/max_num
k = count/max_num
area_total = (x_high-x_low) * (func_method(x_high)-0)
area_x3 = k * area_total
return area_x3
area = cac(10000,1,2,lambda x:x**3)
print('定积分面积为:{:}'.format(area))
def monte_carlo(total_points, x_range, power):
""" 幂函数曲线下方的点个数占总点数的比例就是面积比
Args:
total_points: 生成点总个数
x_range: 矩形x的取值范围
power: 幂次数
Returns:
定积分面积,并以list形式返回所有的x和y点
"""
# 1:定义曲线下方点个数的计数器与点容器
num_curve_points = 0
x_list = []
y_list = []
color_list = []
# 2:生成total_points个随机点
for _ in range(total_points):
# 3:在长宽均为1的矩形内生成随机点(x,y)
rand_x = np.random.uniform(x_range[0], x_range[1], size=1)
# y起始值为0,因为曲线在x轴上方,这里之前误以为y起始值是x起始值的power次方
rand_y = np.random.uniform(0, x_range[1] ** power, size=1)
x_list.append(rand_x)
y_list.append(rand_y)
# 4:判断随机点是否在曲线下方,如果在曲线下方则计数器加1
if rand_x ** power >= rand_y:
num_curve_points += 1
color_list.append("orange")
else:
color_list.append("y")
# 5:根据曲线下方点个数和总点数的比例就是面积比
rectangle_area = (x_range[1] - x_range[0]) * (x_range[1] ** power - 0)
curve_area = (float(num_curve_points) / float(total_points)) * rectangle_area
return curve_area, x_list, y_list, color_list
area_com = monte_carlo(1000,[1,2],3)
print('网上计算结果:{:}'.format(area_com[0]))
蒙特卡洛方法
最新推荐文章于 2022-01-16 19:11:17 发布