一个均匀硬币不断抛啊抛,一直抛到出现HTT(正反反)和HTH(正反正),哪个需要抛的次数多?

感觉次数应该是一样的,从理论上来说。

有人写了代码:

<br />
coin.seq = function(v) {<br />
    x = NULL<br />
    n = 0<br />
    while (!identical(x, v)) {<br />
        x = append(x[length(x) - 1:0], rbinom(1, 1, 0.5))<br />
        n = n + 1<br />
    }<br />
    n<br />
}<br />
set.seed(1024)<br />
mean(htt <- replicate(1e+05, coin.seq(c(1, 0, 0))))<br />
# [1] 8.003<br />
mean(hth <- replicate(1e+05, coin.seq(c(1, 0, 1))))<br />
# [1] 10.006<br />


我觉得不对,请各位看看。
</p>

设总共等待x步, 且当首两掷结果是HH时总共等待y步

<br />
     HTH      HTT<br />
HH   y        y<br />
HTH  3        <3+x<br />
HTT  3+x      3<br />
T    1+x      1+x<br />
</p>

所以HTT早一些

回复 第1楼 的 kuanguang:http://dclong.github.com/en/2010/05/how-long-observe-pattern/

2楼的推理是从假设两个x相等开始推出矛盾

但如果x不一样, y也不会一样

具体计算:

HTH:

y/4 + 3/8 + (3+x)/8 + (1+x)/8 = x

又发现y-1实际上是首掷为H条件下的等待时间, 所以还有方程

(y-1)/2 + (x+1)/2 = x

解得y=x=10

HTT:

当首3掷为HTH条件下的等待时间为2+y-1 = y+1, 所以有方程

y/4 + (y+1)/8 + 3/8 + (1+x)/2 = x

首掷为H条件下的等待时间有方程

(y-1)/2 + (x+1)/2 = x

解得y=x=8

上面的代码可读性不好,感觉有错误,

这句话,x = append(x[length(x) - 1:0], rbinom(1, 1, 0.5)),不大对劲,我修改了代码,可以运行,结果也差不多。

<br />
toss.coin = function(v) {<br />
    x = NULL<br />
    n = 0<br />
    while (!identical(x, v)) {<br />
        n=n+1<br />
        if(length(x)  != 3)<br />
           x=append(x,rbinom(1, 1, 0.5))<br />
        else<br />
           x = append(x[2:3],rbinom(1, 1, 0.5))<br />
    }<br />
    return (n)<br />
}<br />
mean(htt <- replicate(1e+05, toss.coin(c(1, 0, 0))))<br />
mean(hth <- replicate(1e+05, toss.coin(c(1, 0, 1))))<br />


还有个问题,我想看运行时间,

system.time(mean(htt <- replicate(1e+05, toss.coin(c(1, 0, 0)))))

可以产生运行时间,但是没有结果输出了,请问,如何做到,既有结果输出,又有运行时间?
</p>

回复 第5楼 的 kuanguang:原来的代码是没错的,而且比你改写的要高效些。

“length(x) - 1:0”的意思是“c(length(x)-1, length(x))”,也就是取向量的最后两个元素。值得注意的是当length(x) < 2的情况:

当 length(x) == 1 时,就成为 x[c(0,1)],而下标0取不到任何元素,得到的只有第一个元素,即结果为 x 本身;

当 length(x) == 0 时,就成为 x[c(-1,0)],此时 x 本身为 NULL,不管下标如何,结果都是 NULL,仍为 x 本身。

个人觉得写这段代码的应是个高手。

回复 第7楼 的 yanlinlin82: 分明是魏大版主的大作啊 [s:11]

http://cos.name/2011/01/different-ways-to-solve-a-tossing-problem/

1 个月 后

回复 第9楼 的 yanlinlin82:是益辉写的。

我的那个pdf中有一个离散序列概率的通用程序,根据ROSS书总结写的:

<br />
coinE <- function(res.seq, p.seq){<br />
E <- 1/prod(p.seq)<br />
n <- length(res.seq)<br />
for(i in 1:(n-1)){<br />
f <- 1*!sum(abs(diff(res.seq,i)))<br />
if(f==1) E <- E + 1/prod(p.seq[1:(n-i)])<br />
}<br />
return(E)<br />
}</p>
<p>> coinE (c(1,0,0), rep(0.5,3)) ## HTT 首次发生时间期望<br />
[1] 8<br />
> coinE (c(1,0,1), rep(0.5,3)) ## HTH 首次发生时间期望<br />
[1] 10<br />
</p>

代码中 res.seq 是输入花样,需用整数表示;p.seq 和res.seq 等长,表示对应花样元素的单次概率.

用该程序计算 HTT 和 HTH 的首次发生时间期望,得到期望值分别为 8、10.

这里想说的是能得到解析解是很幸福的,模拟在数量大的时候速度太慢。

<br />
> ## 抛掷均匀骰子,得到 123456123456 首次时间期望<br />
> coinE (c(1:6,1:6), rep(1/6,12))<br />
[1] 2176828992<br />
> ## 抛掷不均匀硬币 (正反面概率分别为 0.4,0.6),得到 HTHTHTHT 首次时间期望<br />
> coinE (rep(1:2,4), rep(c(0.4,0.6),4))<br />
[1] 395.2739<br />
</p>

这个作个8个状态的马尔科夫应该可以算出来吧。