Dec 292012
 
This entry is part 8 of 10 in the series Matlab基础班附加材料

在Matlab中rand(), randn()等等,所有的随机数函数生成的都是所谓的"伪随机数"(pseudo-random number),特点就是它们看上去随机,但却都是用某种特殊的算法计算出来的,前一个随机数确定后,后面那个数字就与前面那个数字有关系了(当然肯定不是简单的线性关系)。

一般应用的时候,这种伪随机数已经够用,但是如果你要求比较特殊,需要"真正"的随机数,那么就需要跳出matlab自带的随机数框架。

问题来了,怎样生成真正的随机数呢?第一,可以确定的,凡是人工生成的都不是随机的,例如计算机的伪随机数。第二,真正的随机数要到自然届中追寻。最经典的就是粒子物理中的数据,例如薛定谔的猫,还有就是常见的气象数据。

当然,个人理解,气象数据的随机性有一些特别,例如,现在科技已经可以把天气预报做得比较准确了,我们可以知道明天某个时间点上面的大致温度,但是,由此我们还是可以得到一个纯随机的数字。例如,明天某时刻温度是xx.xxxxxxx,假设xx.x都已经可以预测得差不多了,但是0.0xxxxxx这些尾巴上面的数字却无法预测,纯随机。

当然,原理与实现存在很长的距离,但是互联网已经可以给我们提供"真正随机数"的服务。一种是用量子物理的数据来生成随机数,具体的网址可能要找一些,我暂时还没有找到实用的网址;第二种是用气象数据来生成随机数,这种就是我们这篇博文的主角。

一、注意事项

在讲应用之前必须先讲讲这个服务的约束。它有流量限制,不遵守规则的人要被封IP!你不能做两件事情:

1、 不能在一直很高频率地使用这个服务。克服的办法比较简单,例如你要做1000次循环,每次循环用到一个随机数,那么你不能把生成随机数的代码放在循环里面,那样你会访问1000次服务器,每次虽然只是生成一个数据,但是还是可能被封IP。避免的方法是先访问服务器获得1000个随机数,储存起来,然后再做循环。

2、每个IP的可以生成的随机数数量有一定限制,可以在Quota页面查询还有多少额度。如果额度用满,要么等一会儿额度自动恢复,要么换一个IP。

二、Matlab的实现

有人已经写了几个函数实现。http://www.mathworks.com/matlabcentral/fileexchange/27942 不过由于使用方法比较简单,我还是自己整个简单的函数 realrand.m 来做。

function [result,quotaconsumption]=realrand(N,x)
if nargin==0
	N=1000	%默认1000个
else
	N=min(1e4,ceil(abs(N)));
end
if nargin<2
	x=2^26	%eps是1e(-16),生成的随机数精度只要 1e(-8)就行。
else
	x=min(1e9,ceil(abs(x)));
end
%quota estimate
quotaestimate=N*ceil(1+log2(x));
%First get Quota
quota1=str2num(urlread('http://www.random.org/quota/?format=plain'));
if quotaestimate>quota1
	error(sprintf('Quota not enough! You get %d left, and the recent call consumes about %d',quota1, quotaestimate))
end
result=str2num(urlread(sprintf('http://www.random.org/integers/?num=%d&min=-%d&max=%d&col=1&base=10&format=plain&rnd=new',N,x,x)));
%result=(x+result)./(2*x);
if nargout>1
	quota2=str2num(urlread('http://www.random.org/quota/?format=plain'));
	quotaconsumption=quota1-quota2;
end

return

用myvar=realrand,就可以得到随机数了。不过它是[0,1]区间内的均匀分布,

三、获得最终想要的随机数

如何由网络生成的到自己想要的分布呢?例如正态分布。可以先得到[0,1]区间上的均匀分布,然后由正态分布的CDF逆函数得到想要的分布。
myrandomnumbers=norminv(myvar,0,1)

就可以得到标准正态分布了。

四、简单检验
用JBtest简单检验一下是否符合正态分布吧。

>> [h,p,jbstats]=jbtest(myrandomnumbers,0.01,2^12)
h =
     0
p =
    0.9570
jbstats =
    0.0875

检验结果当然是比较好啦。。。

No related posts.

Series Navigation<< Matlab官方文档的错误-apiext.pdf-11-4OpenCL-Toolbox的安装(Matlab2012A+nVidia 610M) >>
Bookmark/FavoritesSina WeiboGoogle+FacebookQQTwitterYahoo BookmarksBaiduDiggEmailGoogle GmailOutlook.comEvernotePrintAIMLinkedInBlogger PostKindle ItOrkutShare

Related Posts:

  One Response to “生成“真正随机”的随机数”

  1. This code can not generate uniformly distributed numbers between 0 and 1...

 Leave a Reply

(required)

(required)


*

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>