计算物理专题----蒙特卡洛积分实战

2023-09-21 01:48:31

Part one 蒙特卡洛积分计算案例

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.stats import norm, kstest

np.random.seed(0)
def integrate(a,b,n=100):
    x = np.random.uniform(a,b,n)
    total = sum(np.exp(x))
    return (b - a) * total / n



Num = 10000
F = np.zeros(Num)
for i in range(Num):
    result = integrate(0,1) 
    F[i] = result    
df = pd.DataFrame({'score':F})
bins = np.arange(1.61,1.83,0.001)  # 箱的边界值
df['score_bin'] = pd.cut(df['score'], bins)
count = df['score_bin'].value_counts()
x = [(i.left + i.right) / 2 for i in count.index]
y = count.values/sum(count.values)
# 绘制散点图
plt.scatter(x,y,s=1,c="green",label="10000 I(100) calcs")

F = np.zeros(Num)
for i in range(Num):
    result = integrate(0,1,n=500) 
    F[i] = result    
df = pd.DataFrame({'score':F})
bins = np.arange(1.61,1.83,0.001)  # 箱的边界值
df['score_bin'] = pd.cut(df['score'], bins)
count = df['score_bin'].value_counts()
x = [(i.left + i.right) / 2 for i in count.index]
y = count.values/sum(count.values)
# 绘制散点图
plt.scatter(x,y,s=1,c="red",label="10000 I(500) calcs")

F = np.zeros(Num)
for i in range(Num):
    result = integrate(0,1,n=1000) 
    F[i] = result    
df = pd.DataFrame({'score':F})
bins = np.arange(1.61,1.83,0.001)  # 箱的边界值
df['score_bin'] = pd.cut(df['score'], bins)
count = df['score_bin'].value_counts()
x = [(i.left + i.right) / 2 for i in count.index]
y = count.values/sum(count.values)
# 绘制散点图
plt.scatter(x,y,s=1,c="blue",label="10000 I(1000) calcs")

########

##Num = 100000
##F = np.zeros(Num)
##for i in range(Num):
##    result = integrate(0,1) 
##    F[i] = result    
##df = pd.DataFrame({'score':F})
##bins = np.arange(1.61,1.83,0.001)  # 箱的边界值
##df['score_bin'] = pd.cut(df['score'], bins)
##count = df['score_bin'].value_counts()
##x = [(i.left + i.right) / 2 for i in count.index]
##y = count.values/sum(count.values)
##plt.scatter(x,y,s=1,c="red",label="10000 I calcs")
####


plt.legend()
plt.savefig("hist-plot-4.jpg")
plt.show()
##print(count)
##import csv
##header = ["fit"]
##rows = [[i] for i in F]
##
##with open('1.csv','w',newline="") as file:
##    writer = csv.writer(file)
##    writer.writerow(header)
##    writer.writerows(rows)

Final Result：1.718

Real Result： 1.71828

ERROR EASTIMATE：10^-4

  • 蒙特卡洛方法的误差通常可以通过增加模拟次数来减小。具体来说，当我们重复运行模拟并计算出相应的输出后，我们可以计算这些输出的均值和方差，并根据这些数据来估计误差值。如果我们希望减小误差，可以增加模拟次数，这样均值和方差都会更加准确，从而使误差更小。

Part 2 重要抽样 分层抽样

        重要抽样积分和分层抽样积分是数值积分中常用的两种方法，以获得更精确的结果。

        重要性抽样积分涉及从不同于原始分布的概率分布中抽样，但对最终结果贡献较大的区域赋予更大的权重。换句话说，该算法选择更有可能对积分结果做出贡献的样本，并赋予它们比其他样本更高的权重。当原始分布难以采样或对最终结果贡献最大的区域很小时，这尤其有用。

        例如，考虑f(x)从a到b的积分。我们可以从与f(x)|成比例的分布中抽样，而不是从区间[a,b]中均匀抽样。这意味着具有较高的|f(x)|值的区域将被更频繁地采样，并在最终结果中被赋予更高的权重。

        分层采样整合则是将整合区域划分为若干子区域，并从每个子区域分别采集样本。这确保了所有区域的采样都是平等的，而不管该区域的被积数的值是多少。

        例如，考虑f(x,y)在区域[a,b] x [c,d]上的积分。我们可以把这个区域分成4个更小的区域，每个区域分别采样。这确保了我们从所有区域得到相等的样本，当被积量在积分区域上表现出高方差时特别有用。

        重要抽样积分和分层抽样积分都是数值积分的重要技术，它们的适用性取决于需要解决的问题。在被积量高度可变或难以采样的情况下，重要性抽样往往更有效，而分层抽样在被积量在整个积分区域表现出高度可变性的情况下是有用的。

        我们将a,b两问的结果用下表展示处理

Num = 5e4,k=5e3,n=10

Simple MC

Stratified sampling

Round 1

1.3741162994281273

1.3741162994281273

Stats(round for 200)

figure

mean

1.3873987598712592

1.387267505078644

std

0.022234190779851584

0.02282746853161396

find the optimal n

Optimal

RESULT

ERROR

N=18

1.3868879367871603

1e-4

import numpy as np
np.random.seed(0)

#exact value 1.38629
def MC(a=0,b=1,n=50000):
    f = lambda x:1/(x+x**0.5)
    x = np.random.uniform(a,b,n)
    return sum(f(x))*(b-a)/n

#stratified sampling
def SS(a=0,b=1,n=10,k=5000):
    x1 = np.linspace(a,b,k+1)[:k]
    x2 = np.linspace(a,b,k+1)[1:]
    SUM = 0
    for i in range(k):
        SUM += MC(a=x1[i],b=x2[i],n=n)
    return SUM

import matplotlib.pyplot as plt
##S = np.array([SS() for i in range(200)])
##M = np.array([MC() for i in range(200)])
##
##print(np.mean(S))
##print(np.mean(M))
##
##print(np.std(S))
##print(np.std(M))
##
##plt.hist(S,label="S")
##
##plt.legend()
##plt.savefig("P3-1.jpg")
##plt.pause(0.01)
##plt.cla()
##plt.hist(M,label="M")
##plt.legend()
##plt.savefig("P3-2.jpg")
##plt.pause(0.01)
S = []
for n in range(5,20):
    S.append(SS(n=n,k=int(50000/n)))
S = np.array(S)

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)


def f(x):
    return 1 / (np.sqrt(x) + x)

# 积分区域划分
num_strata = 5000
strata_size = 1 / num_strata

# 每个子区域采样点数
num_samples_per_stratum = 10

# 计算每个子区域的采样点
samples = []
for i in range(num_strata):
    stratum_start = i * strata_size
    stratum_end = (i + 1) * strata_size
    stratum_samples = np.random.uniform(stratum_start, stratum_end, num_samples_per_stratum)
    samples.extend(stratum_samples)

# 计算每个采样点的函数值
values = f(samples)

# 计算积分近似值
integral_approx = np.mean(values)

print("Approximate integral value:", integral_approx)

Part 3 积分实战

CHOOSE

ERROR

N^(-1/2)

Best when alpha = 0.5

0.00167

0.3536

Num of Nodes is 8

Best when alpha = 0.6

0.04752

0.3015

Num of Nodes is 8

Best when alpha = 0.7

0.00167

0.3536

Num of Nodes is 8

Best when alpha = 0.8

0.00167

0.3536

Num of Nodes is 8

Part 6 高维积分

        Now,we will implement a MC method to find the volume of the N-dimensional sphere for the arbitrary dimension N and the sphere radius R.

        让我们估算一下一个25D的球体需要多少次实验才能得到精确的结果

1/(((1/25)**0.5)**25)=2.980232238769527e+17

        所以我们不可能从0-1开始取值，这样的运算量是我们不能接受的。

        利用C语言，我们创建了多个线程尽可能得将数值计算的极限推进到了13阶左右（计算用时在一分钟左右）

        同时我们希望使用Python语言进行编程，方便我们进行并行运算：并行运算的程序我们已在附件中给出了，遗憾的是，该方法仅能支持13阶左右的运行，远远达不到25阶的需求。（计算用时在一分钟左右）

        我们需要解决维度灾难问题。应该要从数学角度解决他。

                首先我们改用L1距离运算。稍微加快计算的速度

                其次，我们使用一个低差异序列（Sobol）序列生成更加均匀的随机数

运算结果

Dim

17

18

19

20

21

22

Theory

0.1409

0.0821

0.0466

0.0258

0.0139

0.0073

Calculate

0.1409

0.0821

0.0466

0.0258

0.0139

0.0073

23

24

25

26

27

0.0038

0.0019

0.0009

0.000466

0.000223

0.0038

0.0019

0.0009

0.000456

0.000198

 

import numpy as np
import scipy.special as spi
import matplotlib.pyplot as plt
np.random.seed(0)
import time


def prod(n):
    p = 1
    for i in range(1,n+1):
        p *= i
    return p

def f(n):
    if n%2==0:
        k = int(n/2)
        return np.pi**(k)/prod(k)
    elif n%2==1:
        k = int((n-1)/2)
        return np.pi**k*prod(k)*2**n/prod(n)

def main(Num=1e8,dim=4):
    assert(dim>=2),"dim should greater than 2!!!"
    Sum = 0
    #start = time.time()
    for i in range(int(Num)): 
        x = np.random.random(dim)
        p = sum(x**2)
        if p<=1:
            Sum += 1
    #end = time.time()
    #print("Used:",end-start)
    #print("Volume:",Sum/Num*(2**dim))
    #print("Theory:",f(dim))
    print(Sum)
    return Sum/Num*(2**dim),f(dim)

Dim = [20,21,22,23,24,25,26]
V = [0 for i in range(len(Dim))]
F = [0 for i in range(len(Dim))]


for i in range(len(Dim)):
    s = time.time()
    V[i],F[i] = main(dim=Dim[i])
    e = time.time()
    print(e-s)

    
plt.scatter(range(len(Dim)),V,c="red")
plt.scatter(range(len(Dim)),F,c="blue")
plt.plot(range(len(Dim)),V,c="red",label="volume")
plt.plot(range(len(Dim)),F,c="blue",label="Theory")
plt.legend()
plt.pause(0.01)


def main(Num=1e8,dim=4):
    Sum = 0
    for i in range(int(Num)): 
        x = np.random.random(dim)
        p = sum(x**2)
        if p<=1:
            Sum += 1
    print(Sum)
    return Sum/Num*(2**dim),f(dim)

main.py 

import datetime
import numpy as np

import multiprocessing as mp
 
    
def final_fun(name,param):
    dim = param
    Num = int(5*1e4)
    Sum = np.random.random((dim,Num))**2
    Sum = sum(Sum[:,:])
    Sum = len(np.where(Sum<=1)[0])
    
    return {name: Sum}
 
 
if __name__ == '__main__':

    dim = 20
    
    start_time = datetime.datetime.now()
 
    num_cores = int(mp.cpu_count())
    
    pool = mp.Pool(num_cores)
    param_dict = {'task1': dim,
                  'task2': dim,
                  'task3': dim,
                  'task4': dim,
                  'task5': dim,
                  'task6': dim,
                  'task7': dim,
                  'task8': dim,
                  'task9': dim,
                  'task10': dim,
                  'task11': dim,
                  'task12': dim,
                  'task13': dim,
                  'task14': dim,
                  'task15': dim,
                  'task16': dim}
    
    results = [pool.apply_async(final_fun, args=(name, param)) for name, param in param_dict.items()]
    
    results = [p.get() for p in results]
 
    end_time = datetime.datetime.now()
    use_time = (end_time - start_time).total_seconds()
 
    print("多进程计算 共消耗: " + "{:.2f}".format(use_time) + " 秒")
    print(results)

利用拉丁超立方体采样生成在高维更加均匀的随机数

拉丁超立方体采样的步骤如下：

        1.对于每个维度，将采样区间均匀地划分为相等的子区间。例如，如果要生成100个样本点，那么将采样区间[0, 1]划分为100个子区间，每个子区间的长度为1/100。

        2.在每个维度的每个子区间内，随机选择一个点作为采样点。可以使用均匀分布的随机数生成器来生成这些随机点。

        3.重复步骤2，直到在每个维度上都选择了一个采样点。

import numpy as np
np.random.seed(0)

def latin_hypercube_sampling(dimensions,num_samples):
    samples = np.zeros((num_samples,dimensions))

    intervals = np.linspace(0,1,num_samples+1)
    
    for i in range(dimensions):
        offsets = np.array([np.random.uniform(intervals[np.random.randint(0,num_samples+1)],intervals[np.random.randint(0,num_samples+1)]) for i in range(num_samples)])
        samples[:,i] = offsets
    
    return samples

if __name__ == '__main__':
    dimensions = 25
    num_samples = 1000
    
    # Perform Latin Hypercube Sampling
    samples = latin_hypercube_sampling(dimensions, num_samples)
    
    # Calculate mean and variance for each dimension
    means = np.mean(samples, axis=0)
    variances = np.var(samples, axis=0)
    
    # Print results
    print("Means:", means)
    print("Variances:", variances)

更多推荐

尚硅谷wepack课程学习笔记

为什么需要使用打包工具？开发时使用的框架、es6语法、less等浏览器无法识别。需要经过编译成浏览器能识别的css、js才可以运行。打包工具可以帮我们编译，号可以做代码压缩、兼容处理、性能优化。常见的打包工具有什么？vite、webpack、glup、gruntwebapck最基本的使用?是一个静态资源打包工具，以一个

@Valid注解的作用及@Valid注解与@Validated的区别

1.@Valid注解导入依赖<dependency><groupId>javax.validation</groupId><artifactId>validation-api</artifactId></dependency><dependency><groupId>org.hibernate.validator</g

面试官：你是怎么理解ES6中Proxy的？使用场景?

🎬岸边的风：个人主页🔥个人专栏:《VUE》《javaScript》⛺️生活的理想，就是为了理想的生活!目录一、介绍二、用法参数handler解析Reflectget()set()deleteProperty()取消代理三、使用场景一、介绍定义：用于定义基本操作的自定义行为本质：修改的是程序默认形为，就形同于在编程语

ConfigMaps-1

文章目录主要内容一.使用YAML文件创建1.在data节点创建了一些键值：代码如下（示例）:2.解释二.使用命令行创建1.创建了一个名为person的键值：代码如下（示例）:2.解释3.创建了一个index.html文件，然后用--from-file来引用代码如下（示例）:4.解释总结主要内容使用YAML文件创建使用命

Python 缓存库

文章目录缓存库缓存库的类型Python中有用的缓存库Python中的Redis缓存库Python中的lru_cache库Python中的其他缓存库总结缓存是一种可以存储数据以供快速访问的内存类型。它是一个小而快速的内存，用于保存经常访问的数据。缓存是至关重要的，因为它可以通过减少系统访问缓慢的主存储器的次数来提高系统性

【从入门到起飞】JavaAPI—System，Runtime，Object，Objects类

🎊专栏【JavaSE】🍔喜欢的诗句：更喜岷山千里雪三军过后尽开颜。🎆音乐分享【如愿】🎄欢迎并且感谢大家指出小吉的问题🥰文章目录🍔System类⭐exit()⭐currentTimeMillis()🎄用处⭐arraycopy()🍔Runtime类⭐创建对象⭐exit()⭐availableProcesso

写一篇nginx配置指南

nginx.conf配置找到Nginx的安装目录下的nginx.conf文件，该文件负责Nginx的基础功能配置。配置文件概述Nginx的主配置文件(conf/nginx.conf)按以下结构组织：配置块功能描述全局块与Nginx运行相关的全局设置events块与网络连接有关的设置http块代理、缓存、日志、虚拟主机等

【计算机网络】网络编程接口 Socket API 解读（6）

Socket是网络协议栈暴露给编程人员的API，相比复杂的计算机网络协议，API对关键操作和配置数据进行了抽象，简化了程序编程。本文讲述的socket内容源自Linuxman。本文主要对各API进行详细介绍，从而更好的理解socket编程。recvrecv()遵循POSIX.1-20081.库标准c库，libc,-lc

Ubuntu 安装 CUDA 与 CUDNN GPU加速引擎

一、NVIDIA（英伟达）显卡驱动安装NVIDIA显卡驱动可以通过指令sudoaptpurgenvidia*删除以前安装的NVIDIA驱动版本，重新安装。1.1.关闭系统自带驱动nouveau注意！在安装NVIDIA驱动以前需要禁止系统自带显卡驱动nouveau：可以先通过指令lsmod|grepnouveau查看no

嵌入式笔试面试刷题(day11)

文章目录前言一、字节流，数据报，报文二、makefile怎么引入库和模块三、多次free一块内存空间会怎么样四、字符操作函数越界会发生什么五、QT中一个信号可以连接多个槽函数吗六、QT中一个槽函数可以对应多个信号吗总结前言本篇文章继续刷题。一、字节流，数据报，报文1.数据报（Datagram）：数据报是一种独立的、特定

Linux：基础开发工具之Makefile和缓冲区的基本概念

文章目录动静态库自动化构建代码缓冲区原理实现具体实现动静态库首先要知道什么是链接：C程序中，并没有定义printf的函数实现,且在预编译中包含的stdio.h中也只有该函数的声明,而没有定义函数的实现系统把这些函数实现都被做到名为libc.so.6的库文件中去了,在没有特别指定时,gcc会到系统默认的搜索路径“/usr

热文推荐
  1. 面向过程与面向对象、面向对象三大特性的介绍和示例

    面向过程：将问题分解成一个个详细的步骤，然后通过函数实现每一个步骤，并依次调用特点：1、适合解决简单的问题，不需要过多的协作和抽象2、关注问题的解决步骤而不是问题的本质3、代码复用性低，扩展性差，不易维护4、只有封装，没有继承和多态面向对象：通过分析问题，分解出一个个对象，然后通过不同对象之间的调用和相互协作来解决问题

  2. 论如何在Android中还原设计稿中的阴影

    每当设计稿上注明需要添加阴影时，Android上总是显得比较棘手，因为Android的阴影实现方式与Web和iOS有所区别。一般来说阴影通常格式是有：X:在X轴的偏移度Y:在Y轴偏移度Blur:阴影的模糊半径Color:阴影的颜色何为阴影但是在Android中却比较单一，只有一个度量单位：Elevation，作为在An

  3. ELK 企业级日志分析系统 ELFK

    一、概述ELK平台是一套完整的日志集中处理解决方案，将ElasticSearch、Logstash和Kiabana三个开源工具配合使用，完成更强大的用户对日志的查询、排序、统计需求。使用ELK日志文件系统的原因：日志主要包括系统日志、应用程序日志和安全日志。系统运维和开发人员可以通过日志了解服务器软硬件信息、检查配置过

  4. ELK集群部署

    目录ELKElasticsearch集群部署部署Elasticsearch软件安装Elasticsearch-head插件ELKLogstash部署ELKKiabana部署Filebeat+ELK部署ELKElasticsearch集群部署Node1节点（2C/4G）：node1/192.168.233.30Elast

  5. 浅谈C++|STL之list+forward_list篇

    一.list基本概念功能:将数据进行链式存储链表（list)是一种物理存储单元上非连续的存储结构，数据元素的逻辑顺序是通过链表中的指针链接实现的链表的组成:链表由—系列结点组成结点的组成:一个是存储数据元素的数据域，另一个是存储下一个结点地址的指针域STL中的链表是一个双向循环链表由于链表的存储方式并不是连续的内存空间

  6. 阿里云无影云电脑详细介绍：价格、使用和功能优势说明

    什么是阿里云无影云电脑？无影云电脑（原云桌面）是一种快速构建、高效管理桌面办公环境，无影云电脑可用于远程办公、多分支机构、安全OA、短期使用、专业制图等使用场景，阿里云百科分享无影云桌面的详细介绍、租用价格、云电脑的优势、使用场景、网络架构、无影云电脑与云服务器的区别以及关于无影云电脑的常见问题解答FAQ：目录阿里云无

  7. 华为云云耀云服务器L实例评测｜用PHP从数据库到后端到前端完整实现一个中秋节祝福语项目

    🏆作者简介，黑夜开发者，CSDN领军人物，全栈领域优质创作者✌，CSDN博客专家，阿里云社区专家博主，2023年6月CSDN上海赛道top4。🏆数年电商行业从业经验，历任核心研发工程师，项目技术负责人。🏆本文已收录于PHP专栏：PHP进阶实战教程，评测专区。🎉欢迎👍点赞✍评论⭐收藏文章目录🚀一、前言🚀二、

  8. TCP详解之重传机制

    TCP详解之重传机制TCP实现可靠传输的方式之一，是通过序列号与确认应答。在TCP中，当发送端的数据到达接收主机时，接收端主机会返回一个确认应答消息，表示已收到消息。但在错综复杂的网络，并不一定能如上图那么顺利能正常的数据传输，万一数据在传输过程中丢失了呢？所以TCP针对数据包丢失的情况，会用重传机制解决。接下来说说常

  9. 【C++】C++11（列表初始化和右值引用）

    前言：C++的发展其实并不是一蹴而就的，他经历了很多阶段，本章开始，我们将进入C++11的学习。我们大致来看一下C++的发展历程：当然在这些之中还发行了其他的版本，C++还在不断的向后发展。但是：现在公司主流使用还是C++98和C++11。相比于C++98/03，C++11则带来了数量可观的变化，其中包含了约140个新