- 【R言語】統計解析フリーソフトＲ 第6章【GNU R】 [無断転載禁止]©2ch.net
563 ：１３２人目の素数さん[sage]：2018/11/09(金) 09:44:07.41 ID:5AnUTlVm - >>559
全部が０のとき、エラーになるので修正
rle01 <- function(x){ # c(0,1,1,1,0,0) => return 3
if(sum(x)==0) return(0) #c(0,0,0,0,0,0) => return 0
else{
r=rle(x) # Run Length Encoding
max(r$lengths[which(r$values==1)]) # max length of value 1
}
}
動作確認
> rle01(x<-rbinom(100,1,0.5)) ; x
[1] 8
[1] 1 1 0 1 0 1 0 0 0 1 1 1 0 1 0 0 0 0 1 0 0 0 1 1 1 1 0 0 1 1 1 1 1 0 0
[36] 1 1 0 1 1 1 0 1 1 0 0 1 1 0 1 1 1 0 1 1 0 0 1 1 0 1 0 1 1 0 1 0 1 0 0
[71] 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 1 0 0 1 0 0 1 1 1 0 1 1
> rle01(x<-rbinom(100,1,0)) ; x
[1] 0
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[36] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[71] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
| - 【R言語】統計解析フリーソフトＲ 第6章【GNU R】 [無断転載禁止]©2ch.net
565 ：１３２人目の素数さん[sage]：2018/11/09(金) 10:40:31.84 ID:5AnUTlVm - >>564
whichは余分でした。
sumの方がrleより高速だと思ったらから、すべて０の場合はrleを呼ばないことにしただけ。
0連続でやるとこれだけ差がつく。
rle01 <- function(x){ # c(0,1,1,1,0,0) => return 3
if(sum(x)==0) return(0) #c(0,0,0,0,0,0) => return 0
else{
r=rle(x) # Run Length Encoding
max(r$lengths[which(r$values==1])] # max length of value 1
}
}
rle012 <- function(x) {
r <- rle(x)
one <- r$values == 1
if (any(one)) max(r$length[one]) else 0
}
> x=rep(0,1e8)
> system.time(rle01(x))
user system elapsed
0.3 0.0 0.3
> system.time(rle012(x))
user system elapsed
7.36 4.52 13.25
>
| - 【R言語】統計解析フリーソフトＲ 第6章【GNU R】 [無断転載禁止]©2ch.net
567 ：１３２人目の素数さん[sage]：2018/11/09(金) 11:01:40.74 ID:5AnUTlVm - 1000本に1本あたる宝クジを100本買って続けて2本あたる確率のシミュレーション解の算出時間比較。
確率の理論値は9.8897353347449091e-05
> system.time(mean(replicate(1e4,rle01(rbinom(N,1,p))>=n)))
user system elapsed
0.61 0.00 0.64
> system.time(mean(replicate(1e4,rle12(rbinom(N,1,p))>=n)))
user system elapsed
1.97 0.00 2.03
| - 【R言語】統計解析フリーソフトＲ 第6章【GNU R】 [無断転載禁止]©2ch.net
568 ：１３２人目の素数さん[sage]：2018/11/09(金) 15:22:38.24 ID:5AnUTlVm - >>567
分数で表すと
890788167367/9007199254740992
9.889735334744909e-05
| - 「安倍が死ねば、日本は幸せになる」という命題があるとき、安倍が死んでないならこの命題は必ず真 [無断転載禁止]©2ch.net
210 ：１３２人目の素数さん[sage]：2018/11/09(金) 20:08:31.89 ID:5AnUTlVm - pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95){
nxx=1001
xx=seq(xMIN,xMAX,length=nxx)
xx=xx[-nxx]
xx=xx[-1]
xmin=xx[1]
xmax=xx[nxx-2]
AUC=integrate(pdf,xmin,xmax)$value
PDF=function(x)pdf(x)/AUC
cdf <- function(x) integrate(PDF,xmin,x)$value
ICDF <- function(x) uniroot(function(y) cdf(y)-x,c(xmin,xmax))$root
hdi=HDInterval::hdi(ICDF,credMass=cred)
print(c(hdi[1],hdi[2]),digits=5)
invisible(ICDF)
}
# N個のクジでr個めで初めてあたった時のN個内の当たり数の推測
Atari <- function(N,r){
pmf <- function(x) ifelse(x>N-r+1,0,(1-x/N)^(r-1)*x/N) # dnbinom(r-1,1,x/N) ; dgeom(r-1,x/N)
# curve((1-x/N)^(r-1)*x/N,0,N)
AUC=integrate(pmf,0,N)$value
pdf <- function(x) pmf(x)/AUC
mode=optimise(pdf,c(0,N),maximum=TRUE)$maximum
mean=integrate(function(x)x*pdf(x),0,N)$value
cdf <- function(x) integrate(pdf,0,x)$value
median=uniroot(function(x)cdf(x)-0.5,c(0,N))$root
print(c(mode=mode,median=median,mean=mean))
pdf2hdi(pdf,0,N,cred=0.95)
}
Atari(100,3)
# 全体N個中当たりS個、1個ずつ籤を引いて当たったらやめる.
# r個めが初めて当たりであったときSの信頼区間を推定するシミュレーション。
atari <- function(N,r,k=1e3){ # k: simlation times
f <- function(S,n=N){which.max(sample(c(rep(1,S),rep(0,n-S))))}
vf=Vectorize(f)
sim=replicate(k,vf(1:(N-r)))
s=which(sim==r)%%(N-r)
s[which(s==0)]=N-r
hist(s,freq=T,col='skyblue')
print(quantile(s,c(.025,.05,.50,.95,.975)))
print(HDInterval::hdi(s))
}
atari(100,3,1e3)
| - 面白い問題おしえて〜な 28問目
144 ：１３２人目の素数さん[sage]：2018/11/09(金) 22:58:54.45 ID:5AnUTlVm - コインを１００回投げて表が連続した最大数が５のとき、表がでる確率の９５％信頼区間を求めよ。
近似解計算で
lower upper
0.2487456 0.6386493
になったけど、自信がない。
|
|