14连号问题引起了很多人的关注，该问题的精确概率很难算，不过用蒙特卡罗模拟到比较简单，美中不足的是精度不好提高。下面的代码仅仅计算存在14连号的概率，不管有几串14连号。代码本身很容易懂，大家看看：
<br />
cont <- function(n = 50000){   <br />
    x <- 0<br />
    for(i in 1:n){<br />
        nums <- sort(sample(1138,514))<br />
        nums1 <- nums[1:501]<br />
        nums2 <- nums[14:514]<br />
        flag <- nums2-nums1<br />
       if(any(flag==13)) x <- x + 1 #这里不区分14连号的串数<br />
        <br />
        nums3 <- nums[1:500]<br />
        nums4 <- nums[15:514]<br />
        flag <- nums4-nums3<br />
        if(any(flag==14)) x <- x - 1 #去掉连号超过15的<br />
    }<br />
x/n<br />
} <br />




我得到的答案是0.0045左右，谁的计算机速度快一点，不妨把模拟次数加大点，如果是服务器，那让它多跑一会看看。

如果限定只有一串14连号或者其他数目的连号，那么代码简化为：
<br />
con2 <- function(n = 50000, con = 14){   <br />
    x <- 0<br />
    for(i in 1:n){<br />
        nums <- sort(sample(1138,514))<br />
        nums1 <- nums[1:(514-con + 1)]<br />
        nums2 <- nums[con:514]<br />
        flag <- nums2-nums1<br />
        if(sum(flag==(con-1))==1) x <- x + 1 ##注意这一句<br />
    }<br />
x/n<br />
} <br />




得到的答案稍小一点。



如果有什么问题，请指出。
其实14连号有很多要明确的地方 比如  这14个连号的  抽出来的顺序是不是连着的
有意思，这个连号判断的语句写得很巧妙。



我把末尾一句改成：

cat(sprintf("%.0f / %.0f = %.8f", x, n, x / n), "\n")

好把结果打印出来。



跑了三个不同数量级的模拟，得到如下结果：
<br />
> system.time(con2())<br />
245 / 50000 = 0.00490000 <br />
   user  system elapsed <br />
   9.69    0.00    9.78 <br />
> system.time(con2(500000))<br />
2334 / 500000 = 0.00466800 <br />
   user  system elapsed <br />
  96.11    0.02   97.19 <br />
> system.time(con2(5000000))<br />
23098 / 5000000 = 0.00461960 <br />
   user  system elapsed <br />
 958.69    0.44  972.29 <br />




结果都不到 0.005，而新闻上的专家算出的 1.4% 到底是怎么算出来的呢？



另外，这种算法没考虑抽出顺序，但就算考虑抽出顺序，那结果也只会比这个更小。
<br />
cont <- function(n = 50000){   <br />
    x <- 0<br />
    for(i in 1:n){<br />
        nums <- sort(sample(1138,514))<br />
        nums1 <- nums[1:501]<br />
        nums2 <- nums[14:514]<br />
        flag <- nums2-nums1<br />
       if(any(flag == 13)) x <- x + 1 <br />
       }<br />
       x/n<br />
       }<br />


应该不受抽出顺序的影响，因为结果公布的时候是按照号的大小顺序。如果把14及以上的都考虑进来，模拟了几次，结果在0.0076~0.0080之间，离0.014还有一段距离，呵呵。专家果然是专家啊。
昨晚人品爆发时上了cos，之后就一直forbidden了，今天也这样。我在博客中也谈了谈：

http://taiyun.cos.name/2009/08/14-consecutive-numbers-probability/



若要算14++连号的，代码差不多，算下来是0.0082左右，这个很多人都得出来了。代码：
<br />
con <- function(n = 50000, con = 14){   <br />
    x <- 0<br />
    for(i in 1:n){<br />
        nums <- sort(sample(1138,514))<br />
        nums1 <- nums[1:(514-con + 1)]<br />
        nums2 <- nums[con:514]<br />
        flag <- nums2-nums1<br />
        if(any(flag==(con-1))) x <- x + 1<br />
    }<br />
x/n<br />
} <br />
con()             ## 有14++连号的概率<br />
con()-con(con=15) ## 只有14连号的概率<br />




当然P(14++)-P(15++)就是只有14连号的概率了。



专家的0.0149的确粗糙了点，肯定是错误的，呵呵。



至于抽出顺序，我觉得应该公布的是最终结果吧，如果14连号是依次连续抽出的，那概率可真的小的可怜。并且我认为这不



应该是一个问题，因为一般都是公布结果而已，没有说是第一个第二个之类的。
我用捆绑法给出一种思路，大家看行不行？

首先1138个数里14个连号的共有1125个，而这14个号在514个号里占了14各位，就有501个位置任其选，剩下的500个号来自剩下的1124个数字。按这种思路，即先从1125个数任选一个放在514个位置中的任意一个位置，然后再从1124个数中任选500个数全排剩下的500个位置，则运用排列组合计算式应该是：

这里没有编辑器，又复制不过来，不好意思！有了这个思想写起来很简单，各位自己写吧！
补充说明一下，分子应该是从1125中任选一个，乘以从1124中任选500个，再乘以501的阶乘（即捆绑后的全排列）。

分母应该是从1138个中任选514进行全排列。
a=1125;b=501;

> c=b/a

> for(i in 1:13){

+ c=c*(b+i)/(a+i)}

> d=1125*c

> d

[1] 0.01499478

>

这个应该是正确答案！
首先要感谢 daxia1015 对这种粗略计算方法进行了解释，让我们了解到新闻中的专家可能采用的计算方法及其存在的错误。



这个方法的问题在于：在选取14个连数位置后，再选择剩下的500个数时，没有考虑后来选择的数字，与前面的14个数再连起来的情况。于是存在很多情况被重复计数（比如1:15 连号的数据，在1:14连号时会被选取，在2:15连号时也会被选取），所以结果比正确结果大。



另外，对于这个错误的算法，虽然这个程序设计的数相对较大，但在 R 中也可以不必用循环做乘一次除一次的方式的计算。而直接可以：

> 1125 * prod(501:514) / prod(1125:1138)

[1] 0.01499478

看一下分子：

> 1125 * prod(501:514)

[1] 8.45406e+40

双精度的数表示10的40次方的数还是没有问题的。
本人又模拟了一下，代码如下：

m=100000

c<-rep(0,m)

for(j in 1:m)

{

x<-1:1138

y<-sample(x,514)

z<-sort(y)



for(i in 1:501){ U<-z[(i+1):(i+13)] ;

                 V<-z[i:(i+12)]+1 ;

                 temp<-(U==V)

                 if(sum(temp)==13) {c[j]=1;break}

               }

}

sum(c)/m

这个计算要花点时间，运行结果是 0.00793，好像模拟次数越大，结果越小！（m=10000时为0.0088），大家可以试一下。