蒙特卡洛算法

//根据上面求得的点来产生
改写后的程序（C++版）
#include <string> #include <iostream> using namespace std; int P[2] = {2, 3}; int K[2] = {63, 40}; long index; double x[2]; double q[2][100]; int d[2][100];
void HaltonSequence() { index = 1; for(int i = 0; i < 2; i++) { long k = index; x[i] = 0; for(int j = 0; j < K[i]; j++) { q[i][j] = (j == 0? 1.0: q[i][j-1])/P[i]; d[i][j] = (int)(k % P[i]); k = (k - d[i][j])/P[i]; x[i] += d[i][j] * q[i][j]; } } }
随机点的产生 伪随机算法
随机数生成算法是一类重要的算法，广泛应用于仿真技术 等场合。 伪随机数是由确定的算法生成的，其分布函数与相关性均 能通过统计测试。与真实随机数的差别在于，它们是由算法 产生的，而不是一个真实的随机过程。一般地，伪随机数的 生成方法很多，有线性同余法，直接法，逆转法等
随机点的产生
随机点的产生 准随机算法
伪随机算法都存在差异性，不均匀性。因此，不要求新 的发生器模拟真实的均匀分布，而力求任意大小的样本（尤 其是小样本）都能满足低差异性。换言之，以牺牲随机性为 代价，换来均匀性的提高，称其为准随机模拟器。 目前有3种准随机序列可用来辅助生成均匀分布随机数， 分别是Halton序列、Sobol序列、Latin超立方体序列。
随机点的产生
Halton序列
以质数为基底产生的序列 假设是以质数2为基底的Halton序列，范围在0～1之间。则从产生的列 为 1/2 1/4 3/4 1/8 5/8 3/8 7/8 1/16 9/16…….. 怎么产生的呢？ 3=1+2 5=1+4 7=3+4 9=1+8 ……….. 假设以质数3为基底，Halton序列如下 1/3 2/3 1/9 4/9 7/9 2/9 5/9 8/9 1/27 10/27 19/27……..
void main() { int i; HaltonSequence(); cout<<"("<<x[0]<<","<<x[1]<<")"<<endl; for(i=0;i<500;i++) { nextPoint(); cout<<"("<<x[0]<<","<<x[1]<<")"<<endl; } }
//给q，d赋值 //q[0]=63,q[1]=40 //d[0]=63,d[1]=40
for(int i = 0; i < K.length; i++) { long k = index; x[i] = 0; for(int j = 0; j < K[i]; j++) { 来演示 q[i][j] = (j == 0? 1.0: q[i][j-1])/P[i]; d[i][j] = (int)(k % P[i]); k = (k - d[i][j])/P[i]; x[i] += d[i][j] * q[i][j]; } } }
500个点的产生
蒙特卡洛算法的介绍 算法描述
以概率和统计理论方法为基础的一种计算方法。将所求解的问题 同一定的概率模型相联系，用计算机实现统计模拟或抽样，以获得问 题的近似解。例如pi的计算：在一个面积为一的正方形里面画一个圆， 半径0.5，与正方形四边相切，在正方形中生成一系列随机点，统计 单位圆内的点数与总点数，（圆面积和正方形面积之比为pi/4=m/n）， 当随机点取得越多（但即使取10的9次方个随机点时，其结果也仅在 前4位与圆周率吻合）时，其结果越准确。
蒙特卡洛算法
一.蒙特卡洛算法的介绍 二.随机点的产生 1.伪随机数的产生 2.准随机数的产生 2. 3.随机点产生程序分析
蒙特卡洛算法的介绍 算法简介
蒙特卡洛算法，也称统计模拟方法，是二十世纪四十年代中期 由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统 计理论为指导的一类非常重要 数值计算方法，蒙特·卡罗方法在金融 工程学，宏观经济学，计算物理学(如粒子输运计算、量子热力学计 算、空气动力学计算)等领域应用广泛。
随机点的产生 在二维中，0～1之间产生的点的序列就是（1/2， 1/3）（1/4，2/3）（3/4，1/9）（1/8，4/9） （5/8，7/9）（3/8，2/9）（7/8，5/9） （1/16，8/9）（9/16，1/27）….
核心代码分析
private static class HaltonSequence { static final int[] P = {2, 3}; static final int[] K = {63, 40}; private long index; private double[] x; private double[][] q; private int[][] d; HaltonSequence(long startindex) { index = startindex; x = new double[K.length]; q = new double[K.length][]; d = new int[K.length][]; for(int i = 0; i < K.length; i++) { q[i] = new double[K[i]]; d[i] = new int[K[i]]; } //Halton序列的产生 //以2 、3为基底产生序列 //使2和3产生的序列中元素的个数相同 //索引项 //x数组定义 //二维数组q定义 //二维数组d定义 //输入参数值，得到序列中的元素 // K.length=2
C/C++语言中伪随机数生成算法实际上是采用了“线性同余法“。 具体的计算如下： Xi = (Xi-1 * A + C ) mod M 其中A,C,M都是常数（一般会取质数）。当C=0时，叫做乘同余法。引出 一个概念叫seed，它会被作为X0被代入上式中，然后每次调用rand()函 seed X0 rand() 数都会用上一次产生的随机值来生成新的随机值。可以看出实际上用 rand()函数生成的是一个递推的序列，一切值都来源于最初的 seed。所 以当初始的seed取一样的时候，得到的序列都相同。C语言里面有 RAND_MAX这样一个宏，定义了rand()所能得到的随机值的范围。在C 里可以看到RAND_MAX被定义成0x7fff，也就是32767。rand()函数里递 推式中M的值就是32767。线性同余法生成的是伪随机数，粗略符合均 匀分布。
//把参数传给k
//这部分通过运算
wenku.baidu.com
double[] nextPoint() { 第二个点 index++; for(int i = 0; i < K.length; i++) { for(int j = 0; j < K[i]; j++) { d[i][j]++; x[i] += q[i][j]; if (d[i][j] < P[i]) { break; } d[i][j] = 0; x[i] -= (j == 0? 1.0: q[i][j-1]); } } return x; } }
void nextPoint() { index++; for(int i = 0; i < 2; i++) { for(int j = 0; j < K[i]; j++) { d[i][j]++; x[i] += q[i][j]; if (d[i][j] < P[i]) { break; } d[i][j] = 0; x[i] -= (j == 0? 1.0: q[i][j-1]); } } }