标签: montecarlo

任何人都可以推荐一本关于c ++中蒙特卡罗算法的好书吗？优选地,应用于物理学,甚至更优选地,物理学的类型是量子力学.

谢谢!

我有一个方法,它使用随机样本来近似计算.这种方法被称为数百万次,因此选择随机数的过程非常有效.

我不确定javas Random().nextInt真的有多快,但我的程序似乎并没有像我想的那样受益.

选择随机数时,我会执行以下操作(半伪代码):

// Repeat this 300000 times
Set set = new Set();
while(set.length != 5)
    set.add(randomNumber(MIN,MAX));

现在,这显然有一个糟糕的最坏情况运行时间,因为理论上的随机函数可以为永恒添加重复数字,从而永远保持在while循环中.但是,数字是从{0..45}中选择的,因此重复的值大部分都不太可能.

当我使用上面的方法时,它只比我的其他方法快40%,这不是近似的,但会产生正确的结果.这大约跑了100万次,所以我期待这种新方法至少快50%.

您对更快的方法有什么建议吗？或许你知道一种更有效的方法来生成一组随机数.

澄清一下,这是两种方法:

// Run through all combinations (1 million). This takes 5 seconds
 for(int c1 = 0; c1 < deck.length; c1++){
    for(int c2 = c1+1; c2 < deck.length; c2++){
     for(int c3 = c2+1; c3 < deck.length; c3++){
        for(int c4 = c3+1; c4 < deck.length; c4++){
         for(int c5 = c4+1; c5 < deck.length; c5++){
             enumeration(hands, cards, deck, c1, c2, c3, …

我正在研究蛋白质结构的蒙特卡罗模拟脚本.我从来没有在蒙特卡洛脚本之前做过.我将大规模扩展这个计划.根据蛋白质xyz坐标我必须定义盒子大小.此框将被划分为大小为0.5 A的网格.根据距离和角度标准,我必须根据Boltzmann概率分布指定点.

我的程序应该通过取0.5 A的网格在每个方向上移动并生成随机点并检查距离和角度的条件.如果满足条件,则根据玻尔兹曼概率分布丢弃该点.

这是我生成随机点的代码

from __future__ import division    
import math as mean    
from numpy import *   
import numpy as np   
from string import *    
from random import *    

def euDist(cd1, cd2):# calculate distance
    d2 = ((cd1[0]-cd2[0])**2 + (cd1[1]-cd2[1])**2 + (cd1[2]-cd2[2])**2)
    d1 = d2 ** 0.5
    return round(d1, 2)

def euvector(c2,c1):# generate vector
    x_vec = (c2[0] - c1[0])
    y_vec = (c2[1] - c1[1])
    z_vec = (c2[2] - c1[2])
    return (x_vec, y_vec, z_vec)


 for arang in range(1000):  # generate random point …

我有一组数字~100,我希望在这个集合上执行MC模拟,基本思路是我完全随机化集合,对前20个值做一些比较/检查,存储结果并重复.

现在实际的比较/检查算法非常快,它实际上在大约50个CPU周期内完成.考虑到这一点,为了优化这些模拟,我需要尽可能快地生成随机集.

目前我正在使用George Marsaglia的Multiply With Carry算法,该算法为我提供了17个CPU周期内的随机整数,非常快.但是,使用Fisher-Yates混洗算法,我必须生成100个随机整数,~1700个CPU周期.这远远超过了我的比较时间.

所以我的问题是有没有其他众所周知/强大的技术来进行这种类型的MC模拟,我可以避免长的随机集生成时间？

我想只是从集合中随机选择20个值,但我必须进行碰撞检查以确保选择了20个唯一条目.

更新:

谢谢你的回复.关于我在帖子之后提出的方法,我还有另外一个问题.问题是,这是否会提供真实的(假设RNG是好的)随机输出.基本上我的方法是设置一个与输入数组长度相同的整数值数组,将每个值设置为零.现在我开始从输入集中随机选择20个值,如下所示:

int pcfast[100];
memset(pcfast,0,sizeof(int)*100);
int nchosen = 0;
while (nchosen<20)
{
    int k = rand(100); //[0,100]
    if ( pcfast[k] == 0 )
    {
        pcfast[k] = 1;
        r[nchosen++] = s[k]; // r is the length 20 output, s the input set.
    }
}

基本上我上面提到的,随机选择20个值,除了它似乎是一种确保没有碰撞的优化方式.这会提供良好的随机输出吗？它很快.

我有一个关于使用Python动态构建直方图的概念性问题.我想弄清楚是否有一个好的算法或可能是现有的包.

我编写了一个运行蒙特卡罗模拟的函数,调用1,000,000,000次,并在每次运行结束时返回64位浮点数.以下是上述功能:

def MonteCarlo(df,head,span):
    # Pick initial truck
    rnd_truck = np.random.randint(0,len(df))
    full_length = df['length'][rnd_truck]
    full_weight = df['gvw'][rnd_truck]

    # Loop using other random trucks until the bridge is full
    while True:
        rnd_truck = np.random.randint(0,len(df))
        full_length += head + df['length'][rnd_truck]
        if full_length > span:
            break
        else:
            full_weight += df['gvw'][rnd_truck]

    # Return average weight per feet on the bridge
    return(full_weight/span)

df是一个Pandas数据帧对象,其列标记为'length'和'gvw',分别是卡车长度和重量.head是两个连续卡车之间的距离,span是桥长.只要卡车列车的总长度小于桥梁长度,该功能就会随意将卡车放在桥上.最后,计算每英尺桥上存在的卡车的平均重量(桥上存在的总重量除以桥长度).

因此,我想构建一个表格直方图,显示返回值的分布,可以在以后绘制.我有一些想法:

1. 继续在numpy向量中收集返回的值,然后在MonteCarlo分析完成后使用现有的直方图函数.这是不可行的,因为如果我的计算是正确的,我只需要7.5 GB的内存(1,000,000,000 64位浮点数~7.5 GB)

2. 初始化具有给定范围和数量的bin的numpy数组.每次运行结束时,将匹配区域中的项目数增加1.问题是,我不知道我会得到的价值范围.设置具有范围和适当的箱尺寸的直方图是未知的.我还必须弄清楚如何将值分配给正确的箱子,但我认为它是可行的.

3. 以某种方式做它在飞行中.每次函数返回一个数字时,修改范围和bin大小.我认为这从头开始写起来太棘手了.

好吧,我打赌可能有更好的方法来处理这个问题.任何想法都会受到欢迎!

在第二个注释中,我测试运行上述函数1,000,000,000次只是为了获得计算的最大值(下面的代码片段).这需要大约一个小时的时间span = 200.如果我运行更长的跨度,计算时间会增加(while循环运行时间越长,用卡车填充桥).你认为有没有办法优化这个？ …

假设我想用蒙特卡罗模拟计算Pi作为练习.

我正在编写一个函数,它(0, 1), (1, 0)随机选取一个正方形中的一个点并测试该点是否在圆内.

import scala.math._
import scala.util.Random

def circleTest() = {
  val (x, y) = (Random.nextDouble, Random.nextDouble)
  sqrt(x*x + y*y) <= 1
}

然后我正在编写一个函数,它将测试函数和试验次数作为参数,并返回测试结果为真的试验部分.

def monteCarlo(trials: Int, test: () => Boolean) =
  (1 to trials).map(_ => if (test()) 1 else 0).sum * 1.0 / trials

......我可以计算Pi

monteCarlo(100000, circleTest) * 4

现在我想知道monteCarlo功能是否可以改进.你怎么写monteCarlo高效和可读？

例如,由于试验的数量很大,是否值得使用view或iterator代替Range(1, trials)而reduce代替map和sum？

Ocaml进程只能使用一个核心,为了使用多个核心,我必须运行多个进程.

是否有任何Ocaml框架可用于并行化蒙特卡罗模拟？

我需要为我的宠物射线追踪器生成随机值的算法进行澄清.
我从一个点发射光线.我有这些光线分布的问题:我需要分布均匀,但它不是......

我现在面临的问题是,在我对结果空间的扭曲之后,最初的均匀分布是不均匀的.

因此,例如,如果极坐标系我生成r和t角.分布不均匀且不均匀:靠近每个极点的空间比靠近赤道的结果密度大得多.原因很清楚:我将均匀分布的点从圆柱空间转换为球形.而且我扭曲了结果.同样的问题是如果我规范化立方体中随机生成的点.

我现在的想法是这样的:我想创建一个四面体,对其顶点进行标准化,将每个面(三角形)与中间的点分开,对其进行标准化并递归重复,直到我有足够的点.然后我稍微"扭曲"这些点.然后我再次将它们标准化.而已.

据我所知,这种方法本身并不是纯粹的数学蒙特卡罗方法,因为除了最后一步之外,我不会在任何步骤中使用随机分布.我不喜欢这种复杂性的解决方案.

任何人都可以建议更简单但仍然

• 随机
• 制服
• 快速
• 简单

谢谢!

编辑:
我需要一个快速的方法,而不仅仅是正确的方法.这就是我问蒙特卡罗的原因.提供的答案是正确的,但不是很快.四面体的方法很快,但不是很"随机"=>不正确.
我真的需要一些更合适的东西.

我一直在尝试使用Python创建一个脚本,让我生成大量的点,用于蒙特卡罗方法,以计算对Pi的估计.我到目前为止的脚本是这样的:

import math
import random
random.seed()

n = 10000

for i in range(n):
    x = random.random()
    y = random.random()
    z = (x,y)

    if x**2+y**2 <= 1:
        print z
    else:
        del z

到目前为止,我能够生成我需要的所有点,但我想得到的是运行脚本以便在以后的计算中使用时产生的点数.我不是在寻找令人难以置信的精确结果,只是一个足够好的估计.任何建议将不胜感激.

我试图将一些python代码翻译成C++.代码的作用是运行蒙特卡罗模拟.我认为Python和C++的结果可能非常接近,但似乎发生了一些有趣的事情.

这是我在Python中所做的:

self.__length = 100
self.__monte_carlo_array=np.random.uniform(0.0, 1.0, self.__length)

这是我在C++中所做的:

int length = 100;
std::random_device rd;
std::mt19937_64 mt(rd());
std::uniform_real_distribution<double> distribution(0, 1);

for(int i = 0; i < length; i++)
{
    double d = distribution(mt);
    monte_carlo_array[i] = d;
}

我在Python和C++中以100x5的随机数运行,然后使用这些随机数进行蒙特卡罗模拟.

在蒙特卡罗模拟中,我将阈值设置为0.5,因此我可以轻松验证结果是否均匀分布.

这是monte carlo模拟的概念草案:

for(i = 0; i < length; i++)
{
    if(monte_carlo_array[i] > threshold)    // threshold = 0.5
        monte_carlo_output[i] = 1;
    else
        monte_carlo_output[i] = 0;
}

由于monte carlo数组的长度为120,我希望1在Python和C++中看到60 秒.我计算了1s 的平均数,并发现,尽管C++和Python中的平均数约为60,但趋势是高度相关的.而且,Python中的平均数总是高于 C++.

我是否知道这是因为我做错了什么,还是仅仅因为C++和Python中随机生成机制的区别？

[编辑]请注意,Python中的RNG也是Mersenne …