抛硬币的问题

kuanguang

一个均匀硬币不断抛啊抛，一直抛到出现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 />

我觉得不对，请各位看看。

easttiger

设总共等待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 />

所以HTT早一些

dclong

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

easttiger

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

kuanguang

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

这句话，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)))))

可以产生运行时间，但是没有结果输出了，请问，如何做到，既有结果输出，又有运行时间？

easttiger

you are not honest in this criticism...

yanlinlin82

回复第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 本身。

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

nan.xiao

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

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

yanlinlin82

回复第8楼的肖楠：果然是高手呀[s:17]

cloud_wei

回复第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 />

代码中 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 />

Microhimalaya

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