蒙特卡洛方法

'''
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]))