くじ引きの順番で当たる確率は変わるのだろうか?

ポテト=ルシフェル

どうも！ポテチ君です！
コンソメではなくうすしお～～～～～
と見せかけての！コンポタ味です(???)
いやガリガリ君か！

はいということでガリガリ君です(違います)
今日は久しぶりに数学の記事でも書こうかなって思って
確率統計の方から気になる話をします！
まあ僕確率苦手なんだけど…

まず例として考える状況は

  1. N枚のくじがある(2<=N)

  2. くじを引く人がx人いる(0<x<=N/2)

  3. 当たりがN枚中a枚ある(x<=a<=N/2)

  4. くじを引いた後戻さない

  5. 引く順番を事前に固定し、引き始めたら変更しない

当たりの枚数は人数分以上あって半数を超えないという条件はここでは大事(全員が当たることも全員がハズれることもあり得るということ)

例として
N=10,a=4,x=3のときを考えてみます
1人目が当たる確率は10枚から4枚のどれかを当てればいいから
4/10=0.4
まあこうなるでしょう
2人目は、1人目が当たったかどうかで確率が変わりますね(なので引く順番は固定しておきました)
1人目が当たっていたとしたら
3/9=1/3
1人目がハズれていたとしたら
4/9
の確率で2人目は当たるわけです
これを1人目が当たる確率と合わせれば2人目が当たる確率は
1人目が当たる確率*(1人目が当たった上で2人目が当たる確率)
 +1人目がハズれる確率*(1人目がハズれた上で2人目が当たる確率)
=0.4*1/3+(1-0.4)*4/9=0.4
おっと?
じゃあ同じように考えると3人目は…
0.4*(1/3 * 2/8 + (1-1/3)*3/8)+(1-0.4)*(4/9 * 3/8+(1-4/9)*4/8)=0.4
お???
順番に関わらず全員同じ確率になりました

もしやこれ条件を満たしていれば全部同じ確率になるのでは?
ということで検証を始めていきます

まずちょっと曖昧だけど結果が合うものから
1人目はどんな状況にしても
a/N
なのは疑いようがないので
2人目のときに
「a/N枚の当たりが減っている」
と考えます
そうすると2人目は
$${\frac{a-\frac{a}{N}}{N-1}=\frac{Na-a}{N(N-1)}=\frac{a(N-1)}{N(N-1)}=\frac{a}{N}}$$
の確率で当たることになります
こんな感じでn人目(1<=n<=x)は
「n*a/N枚の当たりが減っている」ので
$${\frac{a-(n-1)*\frac{a}{N}}{N-(n-1)}=\frac{a}{N}}$$
ということで示されてそうな感じがします
問題は「a/N枚の当たりが減っている」というのが曖昧なことです

他に式に起こしてみたりしたんですけど
どうにも=a/Nにできなくて

$${f(N',a')=\frac{a'}{N'}f(N'-1,a'1)+\frac{N'-a'}{N'}f(N'-1,a') (N'>N-xのとき)}$$
$${\frac{a'}{N'} (N'=N-xのとき)}$$
て感じにしてみたり
これはプログラムにするとき書きやすいけど再帰だからちょっとメモリが心配で

なのでまあ結局は地道に全パターンを計算するのが確実

言ってもそうなると30人のときは1073741824回計算することになるんだけど

少なくとも落ちることはないかなって

ということで

N<-100
x<-25
a<-45
pattern<-logical(x+1) #各々が当たったかどうか
data<-numeric(x) #各々が当たる確率
for(i in 0:(x-1)){
  print(i)
  for(j in 0:i){#patternを初期化
    pattern[j+1]<-FALSE
  }
  doNext<-TRUE
  p<-0 #確率
  while(doNext){
    if(pattern[i+1]){ #当たるときだけ処理
      trueCount<-0 #今の時点での当たった数
      pAdd<-1.0 #確率の追加分
      if(i-1>=0){
        for(j in 0:(i-1)){
          if(pattern[j+1]){
            pAdd<-pAdd*(a-trueCount)/(N-j)
            trueCount<-trueCount+1
          }
          else{
            pAdd<-pAdd*(1-(a-trueCount)/(N-j))
          }
        }
      }
      p<-pAdd*(a-trueCount)/(N-i)+p
      doNext<-FALSE
      for(j in 0:i){ #全パターン終わったらループを終了
        if(!pattern[j+1]){doNext<-TRUE}
      }
    }
    j<-0
    while(j<x){ #1つ次の状態へ
      pattern[j+1]=!pattern[j+1]
      if(pattern[j+1]){break}
      j<-j+1
    }
  }
  data[i+1]<-p
}
barplot(data,names.arg=1:x,ylim=c(0,0.5),ylab="確率",xlab="順番")

Rで作ってみました
普通にPythonで作った方がやりやすかっt(

これは100枚中45枚が当たりで、25人が引くときですね

実行結果は

100枚中45枚が当たりのくじ

ということでちゃんと同じ確率になってくれましたね
普通に処理1分ぐらいかかった気がする

これは式の検証なので実際にシミュレーションもしていきますか

N<-100
x<-25
a<-45
count=integer(x) #当たった数
loop=1000000 #試行回数
for(i in 1:loop){
  a_tmp<-a
  for(j in 1:x){
    if(floor(runif(1,1,N-j+2))<a_tmp){
      a_tmp<-a_tmp-1
      count[j]<-count[j]+1
    }
  }
}
barplot(count/loop,names.arg=1:x,ylim=c(0,0.5),ylab="確率",xlab="順番")

この結果もちゃんと

シミュレーション結果

ここで気になる検証
当たりが少ないときはどうなるのかなって
a=1にしてみる

100枚中1枚が当たりのとき

まさかの全部大体一緒という結果

じゃあ逆にa=80してみるとってそれって当たりとハズレが逆転しただけだから

100枚中80枚が当たりのとき

まあこうなるか

なので
順番はほんとに関係ないらしい

じゃあ、その人が判断して次の人が引くか自分が引くか決めることができるようにしてみる(回された次の人は自分が引くかどうか決めることはできない)

N<-100
x<-25
a<-45
count=integer(x) #当たった数
loop=1000000 #試行回数
for(i in 1:loop){
  a_tmp<-a
  wait<-1
  nextWait<-2
  for(j in 1:x){
    index<-wait
    if(a_tmp/(N-j+1)<0.45 && j<x){
      index<-nextWait
      nextWait<-nextWait+1
    }else if(j==x){
      index<-wait
    }else{
      wait<-nextWait
      nextWait<-nextWait+1
    }
    if(floor(runif(1,1,N-j+2))<=a_tmp){
      a_tmp<-a_tmp-1
      count[index]<-count[index]+1
    }
  }
}
barplot(count/loop,names.arg=1:x,ylim=c(0,0.6),ylab="確率",xlab="順番")

判断を決めるのは自分が0.45(平均)以上で引けるときに引くという感じ

さて結果は…!

まじか…
逆にどうしたらばらつくんだ?
判断基準を0.5に引き上げて見ても変わらず
逆に不利なときに引くようにしても変わらず…

今度はパスされた次の人もパスできるようにしてみる

N<-100
x<-25
a<-45
count=integer(x) #当たった数
loop=100000 #試行回数
for(i in 1:loop){
  already<-logical(x)
  a_tmp<-a
  end<-25
  for(j in 1:x){
    index<-1
    nextDo<-TRUE
    while(index<end && nextDo){
      for(k in index:x){
        if(!already[k]){
          index<-k
          break
        }
      }
      if(a_tmp/(N-j+1)<0.45 && index<end){
        nextDo<-TRUE
        index<-index+1
      }else{
        nextDo<-FALSE
      }
    }
    already[index]=TRUE
    for(k in x:1){
      if(!already[k]){
        end<-k
        break
      }
    }
    if(floor(runif(1,1,N-j+2))<=a_tmp){
      a_tmp<-a_tmp-1
      count[index]<-count[index]+1
    }
  }
}
barplot(count/loop,names.arg=1:x,ylim=c(0,0.6),ylab="確率",xlab="順番")

これで今度こそ…!

う、うわあああああああああ！！！

この記事が気に入ったらサポートをしてみませんか？