博文

[转贴]关于随机数生成(2005-12-19 11:59:00)

摘要:转自:http://www.gro.cn/dll.asp?url=3895/3895954.htm
很多学计算机的学生大学四年学下来连随机数是怎么产生的都弄不明白,更有甚者连随机函数怎么调用都搞不清楚,这不能不说是一种悲哀。嗯,现在让我们一起使这种悲哀成为历史吧^_^ 。首先我们要搞清楚一个概念,生成单个的随机数是没有意义的,我们所说的随机数生成,其实都是生成一个随机序列。鉴于人们通常对具体的做法比对分析要感兴趣,我会先给出一系列随机数生成算法,后面我会对这些算法和其它命题进行稍稍深入的讨论。
要找到真正的随机数来源很困难,象离子辐射事件的脉冲检测器,气体放电管和带泄露的电容,我们不可能给每台需要产生随机数的电脑配这么一套装置,况且这些东东产生的数值的随机性和精确性都有问题。所以我们只能考虑通过某种算法来产生随机数。算法都是确定的,因此我们无法产生真正统计随机的数值序列,但是,如果算法很好,所得的序列就可以通过许多随机性测试,这些数就是所谓的伪随机数了。
现在用得最广泛的伪随机数产生算法就是所谓的线性同余算法。其随机数序列{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&nb......

阅读全文(19007) | 评论:0

规格化的随机变量(2005-12-18 15:00:00)

摘要:规格化就是把随机变量的取值限制在[0,1)之间。 尝试过从双精的标准出发得到随机数,但那样很难保证其等概率性,于是还是从rand()的标准出发吧。 ANSI C给出了标准的随机变量计算方法和函数。其算法可以查相关资料,现只给出函数原型: int rand(void); 返回值为int型,在WIN32里是4个字节,但实际上有用部分只有前15bit,生成的随机数是一个0~0x7FFF(十进制32767)的随机正整数。这里假设由标准算法生成的随机数的这15个bit是等概率得到0或1的,那么就用它‘拼’成更高一些位数的随机数,于是给出下面的算法: double drand(void)
{
   return (double)(
            rand()<<15 | rand()
           ) /
          (double)0x40000000;
} rand()<<15 | rand() 将拼成一个30bit的随机整数,取值0~3FFFFFFF,于是相除之后将得到大于等于0小于1的小数。
使用
得到任意范围[a,b)内的随机浮点数:
double i=a+drand()*(b-a);
float  j=a+(float)(drand()*(b-a)); 得到任意范围[a,b]内的随机整数:
int i=a+(int)(drand()*(1+b-a));
测试 接合前面的算法,给出如下的测试用例:
double frand(double (*p)(double),double h)
{
  double x=drand(),y=drand()*h;
  if(y<p(x))return x;
  return frand(p,h);......

阅读全文(7138) | 评论:1

满足任意概率密度函数分布的随机变量生成算法(2005-12-18 15:03:00)

摘要:概率密度的定义
  概率简单地说就是在空间某处附近一确定的范围内发生某事件的可能性,而概率密度则是概率与空间范围之比的极限。在讨论更复杂一些的空间之前,可以从一维的线性空间出发定义它们。
  在一维空间中,用数轴表示所有的点,定义Ф(x)为:随机变量i,当i<x时的概率,因此得Ф(-∞)=0,Ф(+∞)=1(因为全概率为1),这个条件是直接由定义得出,称为归一化条件。
  由此当x取数轴上的点时,随机变量i取到[x,x+△x](△x>0)中时的概率应该为Ф(x+△x)-Ф(x),则定义i取到这个区间的平均概率密度p^=(Ф(x+△x)-Ф(x))/ △x,此时如果极限:
  lim(Ф(x+△x)-Ф(x))/ △x=p(x) 存在
  △ x->0
  则称极限值为这一点的概率密度,由数学知识这正是导数的定义,于是有p(x)=dФ/dx,于是有△Ф=∫p(x)dx,积分区域为△,虽然这是从一维情况推出,这是概率的一般形式定义。于是归一化条件就相应为∫p(v)dv=1,积分区域为全空间。
  由此概率密度函数就唯一决定了随机变量的分布,反过来我们也说随机变量服从某种概率密度函数的分布。 如何由p函数模拟随机变量的生成
  首先我们要有一个这样的随机函数,C的原型如下:
double random();/*返回一个大于等于0,小于1之间的等概率随机双精度小数*/
  -->规格化的随机变量<--
  这个算法的实现将在以后给出算法,事实上笔者最开始学习编程时就是学的这种随机函数,它是BASIC里的标准函数,这里只是为了讨论的需要,因为由它可以非常方便的得到任意范围内的随机数,包括浮点数或整数。
  然后我们再拿来这样一个任意的概率密度函数p(x),且满足归一化条件∫p(x)dx=1,积分区域为[0,1)。最后再给出这个函数在[0,1)上的上限(最大值)h。然后用下面的方法生成服从该分布的随机变量。
  首先得到一个等概率的随机的二维的点(x,y),它落在矩形区域[0,1;0,h]内,同时我们还要求这样一个条件y<p(x),如果得到的点不满足这个条件则重新取点,也就是说把那些不满足此条件的点‘去掉’。最后我说这里的x满足p(x)的分布。证明如下:
......

阅读全文(11119) | 评论:1

[共享]单双精度浮点数的IEEE标准格式(2005-12-17 14:20:00)

摘要:目前大多数高级语言(包括C)都按照IEEE-754标准来规定浮点数的存储格式,IEEE754规定,单精度浮点数用4字节存储,双精度浮点数用8字节存储,分为三个部分:符号位、阶和尾数。阶即指数,尾数即有效小数位数。单精度格式阶占8位,尾数占24位,符号位1位,双精度则为11为阶,53位尾数和1位符号位,如下图所示: 单精度浮点数存储格式 s 指数 尾数 31  30        23 22   0 双精度浮点数存储格式 s 指数 尾数 63 62       52 51   0    细心的人会发现,单双精度各部分所占字节数量比实际存储格式都了一位,的确是这样,事实是,尾数部分包括了一位隐藏位,允许只存储23位就可以表示24位尾数,默认的1位是规格化浮点数的第一位,当规格化一个浮点数时,总是调整它使其值大于等于1而小于2,亦即个位总是为1。例如1100B,对其规格化的结果为1.1乘以2的三次方,但个位1并不存储在23位尾数部分内,这个1是默认位。         阶以移码的形式存储。对于单精度浮点数,偏移量为127(7FH),而双精度的偏移量为1023(3FFH)。存储浮点数的阶码之前,偏移量要先加到阶码上。前面例子中,阶为2的三次方,在单精度浮点数中,移码后的结果为127+3即130(82H),双精度为1026(402H)。         浮点数有两个例外。数0.0存储为全零。无限大数的阶码存储为全1,尾数部分全零。符号位指示正无穷或者负无穷。 下面举几个例子: 单精度浮点数  十进制 规格化 符号 移阶码 尾数               &n......

阅读全文(15097) | 评论:0

多重积分的数值计算方法(2005-10-30 16:12:00)

摘要:只记一点点: 一般的数值积分,只能算简单的定积分,考虑到黎曼积分的多重积分应该如何算,数学上的方法是先化成上下限函数,拆成两个定积分,化成二次积分,但这不得不先进行符号运算求原函数,怎样可以避免求原函数? 其实是把数值积分的方法扩充,对其进行嵌套求解,关键的地方还是对区域D的划分。这样就可以简单的求解了,但是时间复杂度会大幅度提高,如果算一个定积分的复杂度是f(n),那算一个二重积分就是f(n)*f(m),因为区域从原来的n,变成了n*m,似乎这是不可避免的,由积分本身的复杂程度决定。......

阅读全文(8681) | 评论:1

[2005-10-22]快考试啦(2005-10-22 16:26:00)

摘要:11月5号,软考考试。。。 做最后的努力啊。 和前一次不同的是,这次信心十足,尽管难度大了很多。加油!......

阅读全文(4497) | 评论:5

最短带权路径问题的解法::Dijkstra & Floyd(2005-10-22 16:06:00)

摘要:最短带权路径问题的解法::Dijkstra & Floyd
  在一个网络中,如果两个结点之间有直接的因果关系,则这两个结点直接连通,在连接
两个结点的弧上标上它的代价或权,值得注意的是这样的代价不一定是对称的,即A到B
的代价不一定等于B到A的代价,实际问题中以行船为例,有顺水和逆水的区别。在图G
中,给出两个结点求这样一条最短的路径,使经过这条路径上的代价之和最小,这就是最短
路径问题。
  如果所有弧上的权都相等,则问题退化为求两个结点间的一条路径使经过的中间结点最
少。比如在一个城市的交通网络中,乘客关心的可能只是旅途中中转的次数,只希望转更少
次数的车。对于这样的问题,运用BFS对图进行层次遍历,可以在O(n)的时间复杂度上
完成搜索,考虑到BFS的层次遍历性质,不难理解,在此不讨论。
  考虑到带权的有向图(无向图可看成特殊的有向图),如何求最短路径呢?
1、 Dijkstra算法
  Dijkstra是最早提出这个算法的大牛:)。这是图论中非常有名的一个算法。
  图采用邻接矩阵的形式描述,M[i][j]表示结点i到结点j间的代价,如果没有直接因果
关系,则为无穷大,计算机中可以用一个很大的数据代替。后面的Floyd算法同样这样描述。
  差点忘了,dijkstra算法只能求出从结点i到其它各结点的最短路径。算法引入这样两
个集合S和T,S是那些已经确定了到I结点的最短路径的结点,在计算过程中为什么可
以确定某些结点已经到达了最短的要求了呢?后面再讲。T为全集U和S的差集,即那些还
未确定最短路径的结点。而且S的初值是{i},T的初值是U-{i}。另外再引入一个标记数组
D[N],其中在某一步D[k]表示当前从i到k的较短路径,D[k]的初值为M[i][k],整个算法
过程如下:
  1、 在T中选择一个D[k]最小的结点k,将k并入S,并从T中去掉,如果T为{}则
  转到3;
  2、 用k结点和T中其余结点进行一遍比较,如果D[i]>D[k]+M[k][i],则用D[k]+M[k][i]
  取代原来的D[i],重复1;
  3、 算法结束,此时D[k]中保存的就是从I到k结点的最......

阅读全文(21400) | 评论:5

树的计数公式推导(2005-10-20 18:12:00)

摘要:树的计数公式推导 首先提出树相似的定义,如果两棵树T1和T2相似,那么须满足下列条件之一:
(1) T1和T2同时为空树
(2) T1和T2不为空树,则T1和T2的所有子树分别依次对应相似(显然这里只讨论有
序树,树的孩子有序)
树的计数问题描述为:给定一棵树的结点总数n,由n个这样的结点能组成多少棵互不相似
的树呢?
为了解决这个问题,首先把树转换成二叉树,因为一棵树可以唯一的转换成与之对应的二叉
树,转换的方法是二叉树中的左孩子表示树中的最左孩子,二叉树中的右孩子则表示树中的
相邻右兄弟,按这样的方法可以得到唯一的一棵二叉树,且它的根结点没有右子树(因为树
的根结点没有兄弟)。
根据一一对应关系,树的计数问题可以转换成等价的二叉树的计数问题,这使讨论更加方便
容易。
设n个结点可以表示的二叉树计为bn,先从简单的开始:
如果n=0,显然b0=1
如果n=1,显然b1=1
如果n=2,b2=2
    N    N
   /      \
  N        N
如果n=3,b3=5
如果n再大一些,就很难绘出来了。
当n>1时可以得到一个递推公式,此时有一个根结点,剩下n-1个结点可以用来构成两个左
右子树,设左子树有i个结点,则右子树就有n-i-1个结点,由左右子树的不对称性(前面
提到只讨论有序树),所以有b[i]*b[n-i-1]种构成方法,当i取值不同结果也不同,所以所有
的构成方法为各种i值情况的和,所以
bn=∑b[i]*b[n-i-1],i从0取到n-1
所以整个数列可以递推的定义为:
b[0]=b[1]=1
b[n]= ∑b[i]*b[n-i-1],i从0取到n-1 为了求出bn的通项,需要引入构造函数B[n],在这之前首先提一下广义的二项式定理。 二项式是这样的形式:
(a+b)^p,a+b为一个代数和形式,p为非0的实数。
作......

阅读全文(9182) | 评论:4

带通配符*,?的模式匹配(2005-10-03 22:48:00)

摘要:  传统的模式匹配都是:串T匹配串S,如果存在S的子串S'==T。
  现在要求的是模糊匹配,串T中含有特殊符'*'和'?','*'匹配任意字符串或空串,'?'匹配任意一个字符,如'*Rick'匹配任意以'Rick'结尾的串和'Rick'串,'?ick'匹配任意以'ick'结尾且长度为4的串。
输入:
1、带特殊符的T串
2、一系列测试串S[i]
输出:
如果S[i]匹配T输出.T.,否则输出.F. 1、'?'很好处理,只要在原有的定位函数中加一点点就行:
int index(char *s,char *t,int pos)
{
    int i=pos,j=0,lens=strlen(s),lent=strlen(t);
    while(i......

阅读全文(11533) | 评论:5

最小二乘法线性回归(2005-09-23 21:28:00)

摘要:在现实中经常遇到这样的问题,一个函数并不是以某个数学表达式的形式给出,而是以一些自变量与因变量的对应表给出,老师讲课的时候举的个例子是犯罪人的身高和留下的脚印长,可以测出一些人的数据然后得到一张表,它反应的是一个函数,回归的意思就是将它还原成数学表达式,这个式子也称为经验表达式,之所以叫经验就是说它不完全是实际中的那样准确,是有一定偏差的,只是偏差很小罢了。 最小二乘法
设经验方程是y=F(x),方程中含有一些待定系数an,给出真实值{(xi,yi)|i=1,2,...n},将这些x,y值代入方程然后作差,可以描述误差:yi-F(xi),为了考虑整体的误差,可以取平方和,之所以要平方是考虑到误差可正可负直接相加可以相互抵消,所以记误差为: e=∑(yi-F(xi))^2 它是一个多元函数,有an共n个未知量,现在要求的是最小值。所以必然满足对各变量的偏导等于0,于是得到n个方程: de/da1=0
de/da2=0
...
de/dan=0 n个方程确定n个未知量为常量是理论上可以解出来的。用这种误差分析的方法进行回归方程的方法就是最小二乘法。 线性回归
如果经验方程是线性的,形如y=ax+b,就是线性回归。按上面的分析,误差函数为: e=∑(yi-axi-b)^2 各偏导为: de/da=2∑(yi-axi-b)xi=0
de/db=-2∑(yi-axi-b)=0 于是得到关于a,b的线性方程组: (∑xi^2)a+(∑xi)b=∑yixi
(∑xi)a+nb=∑yi 设A=∑xi^2,B=∑xi,C=∑yixi,D=∑yi,则方程化为: Aa+Bb=C
Ba+nb=D 解出a,b得: a=(Cn-BD)/(An-BB)
b=(AD-CB)/(An-BB) C++程序: #include
#include
void main()
{
   double x,y,A=0.0,B=0.0,C=0.0,D=0.0,delta;
   int n;
   cin>>n;
   for(int i=0;i<n;++i)
......

阅读全文(36285) | 评论:14