蒙特卡洛方法

'''
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]))
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值