关于数值概率算法及随机数

一.关于随机数你了解多少？

           要找到真正的随机数来源很困难，象离子辐射事件的脉冲检测器，气体放电管和带泄露的电容，我们不可能给每台需要产生随机数的电脑配这么一套装置，况且这些东东产生的数值的随机性和精确性都有问题。所以我们只能考虑通过某种算法来产生随机数。算法都是确定的，因此我们无法产生真正统计随机的数值序列，但是，如果算法很好，所得的序列就可以通过许多随机性测试，这些数就是所谓的伪随机数了。          

线性同余法是产生伪随机数的最常用的方法。由线性同余法产生的随机序列a0,a1,…,an满足

其随机数序列{Xn}由方程：Xn+1 = ( aXn + c ) mod m

得到，其中m>0称为模数，0≤ a <m称为乘数，0≤c <m称为增量，0≤X0<m称为初始值或种子，当m、a、c、X0都是整数时，通过这个方程就能产生一系列[0，m)范围内的整数了。

很明显，对一个线性同余随机产生算法来说，最重要的是m、a、c的选择。我们希望产生的序列足够长，能包含[0，m)内所有的数，并且产生的数是随机的，最好能用32bit算术高效实现。于是乎为了有足够的空间产生足够长的队列，我们会使m足够大；为了使序列足够完整，足够象随机序列，我们会发现当m是素数，c为0时就可以使序列的周期（即序列不重复的长度）对于一定的a达到m-1；要能用32bit算术高效实现，我们的m当然不能超过32bit，幸好刚好有一个满足上述条件的很方便的素数2^31-1给我们用，于是产生函数就变成了Xn+1 = ( aXn) mod ( 2^31 – 1 )。在超过20亿个a的可选取值中，只有很少的几个可以满足上述所有条件，其中之一就是a = 7^5 = 16807。

unsigned myrand()    {        return (seed = (seed * 16807L) & 0x7fffffffL);    }

这种产生器得到了广泛的使用和比其它产生器都彻底的测验，其产生的随机数有在分布上有很好的特性，但是却不具备良好的不可预测性。所以我们在使用时往往会以系统时钟（模m）或系统时钟与随机数的和（模m）作为种子来重新启动序列，这也是现在绝大多数随机数产生器的做法。

另一种产生伪随机数序列的方法就是使用加密逻辑，通过加密算法的雪崩效应，输入的每一位改变就会造成输出的巨大改变，并且假如密钥被妥善保存的话从序列中前面的数推出后面的数在计算上是不可行的，于是我们可以使用一个简单的累加器来作为输入，得到的随机序列可以直接使用（当然也可以把系统时间等等因素考虑进去作为输入来产生更随机的序列，不过这通常是没有必要的），这种方法常常用来产生会话密钥或是现时。具体采用的加密逻辑和运算方式有很多，这里就不一一介绍了，大家感兴趣可以找本密码学的书来看看。
        在密码学领域里比较流行的一种产生安全的伪随机数的方法是BBS产生器，它用其研制者的名字命名（Blum Blum Shub，不是我们经常灌水的那个BBS），它的密码编码强度据说具有最强的公开证明。首先，我们找两个大素数p和q，要求它们被4除都余3，令n = p×q，再选择一个随机数s，使s与n互素，然后我们通过如下算法产生一系列比特Bi：

X0 = (s^2)mod n，for i = 1 to ∞Xi = (Xi-1^2)mod nBi = Xi mod 2

每次迭代都只取出最低位的bit。可以看出BBS产生随机数的过程非常复杂，运算量很大，嗯，这就是为安全性付出的代价。

二.数值概率算法

参考博客： http://www.cnblogs.com/chinazhangjie/archive/2010/11/11/1874924.html

1、用随机投点法计算pi值

　　设有一半径为r的圆及其外切四边形。向该正方形随机地投掷n个点。设落入圆内的点数为k。由于所投入的点在正方形上均匀分布，因而所投入的点落入圆内的概率为(PI * pow(r,2)) / (4 * pow(r,2)) = PI / 4 。所以当n足够大时，k与n之比就逼近这一概率。从而，PI 约等于 (4*k)/n.如下图：

实现：

#include <iostream>#include <cstdlib>#include <limits>using namespace std;// 获得0-1之间的随机数double get_random_num (){    return (double)rand () / RAND_MAX ;}// 用随机投点法计算 PIdouble darts (int n){    int k = 0 ;    for (int i = 0; i < n; ++ i) {        double x = get_random_num() ;        double y = get_random_num() ;        if ((x * x + y * y) <= 1.0) {            ++ k ;        }    }    return (4 * k) / (double)n ;}int main(){    cout << darts (200000000) << endl ;}

2.计算定积分

设f(x)是[0，1]上的连续函数，且0 <= f(x) <= 1。

需要计算的积分为，积分I等于图中的面积G。

在图所示单位正方形内均匀地作投点试验，则随机点落在曲线下面的概率为

假设向单位正方形内随机地投入n个点(xi,yi)。如果有m个点落入

G内，则随机点落入G内的概率

假定 f(x) = x * x (0 <= x <= 1)计算定积分

实现：

代码#include <iostream>#include <cstdlib>using namespace std;// 获得0-1之间的随机数double get_random_num (){    return (double)rand () / RAND_MAX ;}// 用随机投点法计算 y = pow(x,2)定积分double darts (int n){    int k = 0 ;    for (int i = 0; i < n; ++ i) {        double x = get_random_num() ;        double y = get_random_num() ;        if (y <= x * x) {            ++ k ;        }    }    return k / (double)n ;}int main(){    cout << darts (10000000) << endl ;    return 0;}

3.解非线性方程组

求解下面的非线性方程组

其中，x1,x2,…,xn是实变量，fi是未知量x1,x2,…,xn的非线性实函数。要求确定上述方程组在指定求根范围内的一组解

在指定求根区域D内，选定一个随机点x₀作为随机搜索的出发点。在算法的搜索过程中，假设第j步随机搜索得到的随机搜索点为x_j。在第j+1步，计算出下一步的随机搜索增量Dx_j。从当前点x_j依Dx_j得到第j+1步的随机搜索点。当x<e时，取为所求非线性方程组的近似解。否则进行下一步新的随机搜索过程。

三、舍伍德(Sherwood)算法

设A是一个确定性算法，当它的输入实例为x时所需的计算时间记为tA(x)。设Xn是算法A的输入规模为n的实例的全体，则当问题的输入规模为n时，算法A所需的平均时间为

 

这就是舍伍德算法设计的基本思想。当s(n)与tA(n)相比可忽略时，舍伍德算法可获得很好的平均性能。

例如：线性时间选择

// 在数组a从l到r之间的子序列上，寻找第k小元素

select(Type a[], int l, int r, int k){

    // 从l到r任意选取一个元素，放在l上，用 pivot 表示

    // 从l + 1到r之间，找到元素j，使得j左边的元素都小于pivot;j右边的元素都大于pivot

    // 那么，或者pivot是要找的元素(k=j-l+1)

    // 或者要找的元素在l,j之间(k<j-l+1)

    // 或者要找的元素在j+1,r之间(k>j-l+1)

}

有时也会遇到这样的情况，即所给的确定性算法无法直接改造成舍伍德型算法。此时可借助于随机预处理技术，不改变原有的确定性算法，仅对其输入进行随机洗牌，同样可收到舍伍德算法的效果。例如，对于确定性选择算法，可以用下面的洗牌算法shuffle将数组a中元素随机排列，然后用确定性选择算法求解。这样做所收到的效果与舍伍德型算法的效果是一样的。

洗牌算法实现：

代码/* 主题：随机洗牌算法* 作者：chinazhangjie* 邮箱：[e-mail]* 开发语言：C++* 开发环境：Code::Blocks 10.05* 时间: 2010.11.11*/#include <iostream>#include <cstdlib>#include <iterator>#include <algorithm>using namespace std;template <class T>void my_shuffle (T* arr, int len){    for (int i = 0; i < len; ++ i) {        int r = rand() % len ;        swap (arr[i], arr[r]) ;    }}int main(){    const int len = 10 ;    int arr[len] ;    for (int i = 0; i < len; ++ i) {        arr[i] = i ;    }    my_shuffle (arr, len) ;    cout << "my_shuffle: ";    copy (arr, arr + len, ostream_iterator<int> (cout, " ")) ;    cout << endl ;    for (int i = 0; i < len; ++ i) {        arr[i] = i ;    }    random_shuffle (arr, arr + len) ;    cout << "stl->random_shuffle: ";    copy (arr, arr + len, ostream_iterator<int> (cout, " ")) ;    cout << endl ;    return 0;}

四、拉斯维加斯(Las Vegas)算法

void obstinate(Object x, Object y){    // 反复调用拉斯维加斯算法LV(x,y)，直到找到问题的一个解y    bool success= false;    while (!success) success = lv(x,y);}

设p(x)是对输入x调用拉斯维加斯算法获得问题的一个解的概率。一个正确的拉斯维加斯算法应该对所有输入x均有 p(x) >0。设t(x)是算法obstinate找到具体实例x的一个解所需的平均时间 ,s(x)和e(x)分别是算法对于具体实例x求解成功或求解失败所需的平均时间，则有 t(x) = p(x) * s(x) + (1-p(x))(e(x) + t(x))

解此方程可得：

n后问题 

　　对于n后问题的任何一个解而言，每一个皇后在棋盘上的位置无任何规律，不具有系统性，而更象是随机放置的。由此容易想到下面的拉斯维加斯算法。

在棋盘上相继的各行中随机地放置皇后，并注意使新放置的皇后与已放置的皇后互不攻击，直至n个皇后均已相容地放置好，或已没有下一个皇后的可放置位置时为止。

　　如果将上述随机放置策略与回溯法相结合，可能会获得更好的效果。可以先在棋盘的若干行中随机地放置皇后，然后在后继行中用回溯法继续放置，直至找到一个解或宣告失败。随机放置的皇后越多，后继回溯搜索所需的时间就越少，但失败的概率也就越大。

12皇后问题，使用带回溯的拉斯维加斯算法

五、蒙特卡罗(Monte Carlo)算法

　　在实际应用中常会遇到一些问题，不论采用确定性算法或概率算法都无法保证每次都能得到正确的解答。蒙特卡罗算法则在一般情况下可以保证对问题的所有实例都以高概率给出正确解，但是通常无法判定一个具体解是否正确。

　　设p是一个实数，且1/2 <p <1。如果一个蒙特卡罗算法对于问题的任一实例得到正确解的概率不小于p，则称该蒙特卡罗算法是p正确的，且称p - 1/2是该算法的优势。

如果对于同一实例，蒙特卡罗算法不会给出2个不同的正确解答，则称该蒙特卡罗算法是一致的。

有些蒙特卡罗算法除了具有描述问题实例的输入参数外，还具有描述错误解可接受概率的参数。这类算法的计算时间复杂性通常由问题的实例规模以及错误解可接受概率的函数来描述。

对于一个一致的p正确蒙特卡罗算法，要提高获得正确解的概率，只要执行该算法若干次，并选择出现频次最高的解即可。

如果重复调用一个一致的(1/2 + e)正确的蒙特卡罗算法2m-1次，得到正确解的概率至少为1-d，其中，

对于一个解所给问题的蒙特卡罗算法MC(x)，如果存在问题实例的子集X使得：

(1)当x不属于X时，MC(x)返回的解是正确的；

(2)当x属于X时，正确解是y0，但MC(x)返回的解未必是y0。

称上述算法MC(x)是偏y0的算法。

重复调用一个一致的，p正确偏y0蒙特卡罗算法k次，可得到一个O(1-(1-p)k)正确的蒙特卡罗算法，且所得算法仍是一个一致的偏y0蒙特卡罗算法。

主元素问题 

设T[1:n]是一个含有n个元素的数组。当|{i|T[i]=x}|>n/2时，称元素x是数组T的主元素。

template<class Type>

bool Majority(Type *T, int n)

{

    // 判定主元素的蒙特卡罗算法

    int i = rnd.Random(n) + 1;

    Type x = T[i];    // 随机选择数组元素

    int k = 0;

    for (int j = 1;j <= n;j ++) {

       if (T[j] == x)

           k ++;

    }

    return (k > n/2);  // k > n / 2 时T含有主元素

}

如果含有主元素，则以>1/2的概率返回true

如果没有主元素，则一定返回false

template <typename Type>

bool Majority2(Type *T, int n)

{

    if (Majority(T,n))

       return true;

    else

       return Majority(T,n);

}

p+(1-p)p=1-(1-p)^2>=3/4

 

template<class Type>

bool MajorityMC(Type *T, int n, double e)

{

    // 重复调用算法Majority

    int k = ceil(log(1/e)/log(2));

    for (int i = 1;i <= k;i ++)

    if (Majority(T,n))

       return true;

    return false;

}

对于任何给定的e>0，算法majorityMC重复调用élog(1/e)ù 次算法majority。它是一个偏真蒙特卡罗算法，且其错误概率小于e。算法majorityMC所需的计算时间显然是O(nlog(1/ e))。