- 臨床統計もおもしろいですよ、その2
71 ：卵の名無しさん[sage]：2018/11/09(金) 07:26:09.53 ID:guOlTCed - (mean(replicate(1e4,any(diff(cumsum(rbinom(100,1,0.5)),5)==5))))
[1] 0.8089
>
> (mean(replicate(1e4,with(rle(rbinom(100,1,0.5)), max(lengths[which(values=
<rle(rbinom(100,1,0.5)), max(lengths[which(values== 1)])>=5))))
[1] 0.8117
>
| - 臨床統計もおもしろいですよ、その2
72 ：卵の名無しさん[sage]：2018/11/09(金) 07:26:39.67 ID:guOlTCed - (mean(replicate(1e4,any(diff(cumsum(rbinom(100,1,0.5)),5)==5))))
[1] 0.8089
>
> (mean(replicate(1e4,with(rle(rbinom(100,1,0.5)), max(lengths[which(values=
<rle(rbinom(100,1,0.5)), max(lengths[which(values==1)])>=5))))
[1] 0.8117
| - 臨床統計もおもしろいですよ、その2
73 ：卵の名無しさん[sage]：2018/11/09(金) 07:27:42.12 ID:guOlTCed - 実行速度
system.time(mean(replicate(1e4,any(diff(cumsum(rbinom(100,1,0.5)),5)==5))))
user system elapsed
1.840 0.000 1.875
>
> system.time(mean(replicate(1e4,with(rle(rbinom(100,1,0.5)), max(lengths[wh
<e(1e4,with(rle(rbinom(100,1,0.5)), max(lengths[whi ch(values==1)])>=5))))
user system elapsed
4.440 0.000 4.631
| - 底辺私立医大を卒業した医者って頭悪いよね？ Part10
426 ：卵の名無しさん[sage]：2018/11/09(金) 07:47:17.21 ID:guOlTCed - coin5 = function(x,n=5) any(diff(cumsum(x),n)==n)
mean(replicate(1e4,coin5(rbinom(100,1,0.5))))
数学板のRのスレにgrep使ったスクリプトをあげたら
別アルゴリズムで高速化してくれた。
頻度が小さい事象のシミュレーションは実行速度が問題になるのでありがたい。
ド底辺シリツ医大の裏口入学調査委員会が
裏口入学は高々10%と報告したとする。
その結果の検証に100人を調査したら4人以上続けて裏口入学生であった、という。
この検証から裏口入学率が10%であるか否かを有意水準5%で検定せよ。
実行速度比較
system.time(mean(replicate(1e5,any(diff(cumsum(rbinom(100,1,0.1)),4)==4))))
user system elapsed
19.210 0.050 21.897
>
> system.time(mean(replicate(1e5,with(rle(rbinom(100,1,0.1)), max(lengths[wh
<e(1e5,with(rle(rbinom(100,1,0.1)), max(lengths[which(values==1)])>=4))))
user system elapsed
42.360 0.020 43.736
| - 底辺私立医大を卒業した医者って頭悪いよね？ Part10
429 ：卵の名無しさん[sage]：2018/11/09(金) 14:58:35.10 ID:guOlTCed - >>420
import numpy as np
from fractions import Fraction
def seqNp(N = 10,K = 5,p = 0.5):
if N == K: return p**K
q = 1-p
a = [0]*(N+1)
for i in range(1,K+1):
a[i] = q/(p**(i))
for i in range(K,N):
a[i+1] = 0
for j in range(0,K):
a[i+1] = a[i+1]+a[i-j]
a[i+1] = q/p*a[i+1]
P0=[0]*(N+1)
for i in range(1,N+1):
P0[i] = a[i]*p**i
del P0[0]
MP = np.zeros([K,N])
MP[0] = P0
MP[1][0] = p
for i in range(K-3,K):
MP[i][0] = 0
for k in range(1,K):
for i in range(0,N-1):
MP[k][i+1] = p*MP[k-1][i]
re = 1-np.sum(MP,axis=0)
print(Fraction(re[N-1]))
print(re[N-1])
| - 底辺私立医大を卒業した医者って頭悪いよね？ Part10
430 ：卵の名無しさん[sage]：2018/11/09(金) 14:59:21.85 ID:guOlTCed - http://tpcg.io/5xvZhU
| - 【セブ島】東医佐野タカノリ先生【橋本環奈】
497 ：卵の名無しさん[sage]：2018/11/09(金) 20:25:11.58 ID:guOlTCed - >>495
裏口入学を除籍にすべき！
| - 臨床統計もおもしろいですよ、その2
84 ：卵の名無しさん[sage]：2018/11/09(金) 21:35:35.51 ID:guOlTCed - 次の課題はこれだな。
コインを100回投げて表が連続した最大数が10であったとき
このコインの表がでる確率の95%信頼区間はいくらか？
|
|