【R言語】統計解析フリーソフトＲ 第6章【GNU R】 [無断転載禁止]©2ch.net
■ このスレッドは過去ログ倉庫に格納されています
再帰呼び出しに関数名の代わりにRecall使うと処理に要する時間が増えることに気づいた。
フィボナッチ数列の再帰呼び出し
# 関数名での呼び出し
fibo <- function(n){
if(n==1|n==2) return(1)
else fibo(n-1)+fibo(n-2)
}
# Recallを使っての呼び出し
fiboR <- function(n){
if(n==1|n==2) return(1)
else Recall(n-1)+Recall(n-2)
}
> system.time(fibo(30))
user system elapsed
4.27 0.02 4.36
> system.time(fiboR(30))
user system elapsed
6.65 0.00 6.70 cumsumを再帰関数で書いてみた。
何度も試行錯誤
# cumsum with recursive call
cumsumR <- function(v){
res=numeric(length(v))
cumsumR_sub <- function(v,res,i){
res[1]=v[1]
if(i==length(v)) return(res)
else{
res[i+1] = res[i] + v[i+1]
Recall(v,res,i+1)
}
}
cumsumR_sub(v,res,1)
}
cumsumR(1:10)
とりえあず、動作した。
> cumsumR(1:10)
[1] 1 3 6 10 15 21 28 36 45 55 簡略化できた
cumsumR <- function(v,res=NULL,i=1){
res[1]=v[1]
if(i==length(v)) return(res)
else{
res[i+1] = res[i] + v[i+1]
Recall(v,res,i+1)
}
}
> cumsumR(1:10)
[1] 1 3 6 10 15 21 28 36 45 55
> ifelseの動作がどうも納得できないでご教示いただきたいのですが。
これってifelseの仕様でしょうか？バグでしょうか？
> x=2:1
> if(x[1]<=x[2]) x else x[2:1]
[1] 1 2
> ifelse(x[1]<=x[2],x,x[2:1])
[1] 1 >>482
ifelse() はベクトル演算できる関数。
だから、戻り値の要素数は、test部分の要素数に一致する。
testが1つなので、x[2:1]の1つ目の要素が返される。 testの要素を2つにすれば、出力も2つになる。
> ifelse(c(x[1] <= x[2], x[1] <= x[2]), x, x[2:1])
[1] 1 2
さらに３つにすると、helpに書いている通り、yesやno部分は再利用される。
> ifelse(c(x[1] <= x[2], x[1] <= x[2], x[1] <= x[2]), x, x[2:1])
[1] 1 2 1 >>483
ありがとうございました。
バグではなくてそういう仕様だったのですね。
こういう使い方ができるということがわかって勉強になりました。
> x=2:1
> ifelse(c(x[1]!=x[2],x[1]==x[2]),1:2,3:4)
[1] 1 4 教わったので早速、
ifelse() はベクトル演算できると、再利用されるの動作確認。
> ifelse(c(TRUE,TRUE,FALSE,FALSE,TRUE,FALSE,FALSE),1:2,11:15)
[1] 1 2 13 14 1 11 12 ２種類の文字の変換もifelseでできるんだなぁ。
> mat
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 1 1 0
[3,] 1 1 0
[4,] 1 1 0
[5,] 1 1 0
[6,] 1 1 0
> print(apply(mat,2, function(x) ifelse(x==1,'真','偽')),quote=F)
[,1] [,2] [,3]
[1,] 真 真 真
[2,] 真 真 偽
[3,] 真 真 偽
[4,] 真 真 偽
[5,] 真 真 偽
[6,] 真 真 偽 行列も次元付きのベクトルだからもっと簡略化できた。
> (mat=matrix(sample(0:1,12,rep=T),3,4))
[,1] [,2] [,3] [,4]
[1,] 1 1 0 1
[2,] 0 1 0 0
[3,] 1 0 0 1
> print(ifelse(mat==1,'Right','Wrong'),quote=F)
[,1] [,2] [,3] [,4]
[1,] Right Right Wrong Right
[2,] Wrong Right Wrong Wrong
[3,] Right Wrong Wrong Right
> 規約分数にするパッケージが探せなかったので自分でスクリプトを組んでみた。
エラー処置は省略。
reduce_fraction <- function(x,y){
a=x
b=y
r = a%%b # a=bq+r -> a%%b=b%%r
while(r){
a = b
b = r
r = a%%b
}
gcd=b
cat(x/gcd,'/',y/gcd,'\n')
invisible(c(x/gcd,y/gcd))
}
> reduce_fraction(2860,1082900)
11 / 4165 rvest?でログイン必要なとこをスクレイピングする時って
login_page <- html_session("https://xxxxx";)
login_form <- html_form(login_page)[[1]] %>%
set_values(AAA="xxxx", BBB="xxxx")
というので行けるはずなんですが
最終行のAAAとBBBってそれぞれhtmlタグのname=""のとこからとってくるんですよね?
これにここがname="user[email]"みたいに[]という記号はいってたらどうすればいいでしょうか? >>491
よろしくおねがいいしまーーーーーーーす >>492
質問の意味がいまいち分からないから、誰も助言できないのでは？
> x <- 'user[email]'
> x
[1] "user[email]"
rvestパッケージを使ったことがないけど、
記号が入っていたらなぜ問題があるのかよく分からない。 というかどこのページのスクレイピングしたいか晒して
仮でもいいから >>493
ちがいます
AAAがどれかよくみてみてください
””はありません ちょっと直接はることはできないのでタグだけはります。
<input class="form-control" autofocus="autofocus" placeholder="Email address" required="required" type="text" value="" name="user[email]" id="user_email">>>494 >>495
もっと分かるように説明しないと、追試できる情報も提供しないし。
もしかして変数名などに[]が入っている場合 にどうしたらよいかってこと？
> `user[email]` <- 1:5
> `user[email]`
[1] 1 2 3 4 5 ある問題のシミュレーションしようと思って
問題文の記号のまま
q <- function(x)
....
とやって
q(100)と入力すると
Rが終了することに気づいた。 ちょっとした疑問です。
空ベクトルの検出って長さ０以外に検出方法ってあるでしょうか？
> x=c(1,2)
> x=x[-1]
> x=x[-1]
> length (x)
[1] 0
> x
numeric(0)
> x==numeric (0)
logical(0)
> is.null(x)
[1] FALSE
> is.na(x)
logical(0)
> x==NULL
logical(0)
> x==NA
logical(0)
> length(x)==0
[1] TRUE ４点の座標を入力するとそれらを結ぶ四面体の体積を求めるスクリプトを書いてみた。
高さはパッケージnleqslvを使った近似計算。
# Calculate tetrahedron volume from cordinates
library(nleqslv)
Tetra <- function(O=c(1/2,sqrt(3)/6,sqrt(2/3)),A=c(0,0,0),B=c(1,0,0),C=c(cos(pi/3),sin(pi/3),0)){
fn <- function(x,O,A,B,C){
AO=A-O
BO=B-O
CO=C-O
HO=x[1]*AO+x[2]*BO+(1-x[1]-x[2])*CO # H on triangle ABC
AB=B-A
AC=C-A
c(HO%*%AB,HO%*%AC) # HO vertial to AB and AC
}
fn1 <- function(x) fn(x,O,A,B,C)
x=nleqslv::nleqslv(c(1/3,1/3),fn1)$'x'
AO=A-O
BO=B-O
CO=C-O
HO=x[1]*AO+x[2]*BO+(1-x[1]-x[2])*CO
h=sqrt(sum(HO^2))
a=sqrt(sum((B-C)^2))
b=sqrt(sum((C-A)^2))
c=sqrt(sum((A-B)^2))
s=(a+b+c)/2
base=sqrt(s*(s-a)*(s-b)*(s-c))
V=1/3*base*h
return(V)
}
初期値は辺の長さ１の正四面体
> options(digits = 16)
> Tetra()
[1] 0.1178511301977579
> sqrt(2)/12
[1] 0.1178511301977579
多分、正常に動作していると思う。 >>502
行列式det使うと簡単
> po <- c(1/2, sqrt(3)/6, sqrt(2/3))
> pa <- c(0,0,0)
> pb <- c(1,0,0)
> pc <- c(cos(pi/3), sin(pi/3), 0)
> det(cbind(pa-po,pb-po,pc-po))/6
[1] -0.1178511 ベクトルの三重積を教わったので、パッケージ pracma の外積crossを使った
tetrahedron <- function(O=c(1/2,sqrt(3)/6,sqrt(2/3)),A=c(0,0,0),B=c(1,0,0),C=c(cos(pi/3),sin(pi/3),0)){
AO=A-O
BO=B-O
CO=C-O
as.numeric(abs(pracma::cross(AO,BO) %*% CO)/6)
}
４行で済んだ。
>>503
ありがとうございました。
パッケージに頼らずに計算できたのですね。 データを解析する前にさらっと特徴を見たい時、皆さんはどんなコマンドを使っていますか？
私が思いつくのは
summary
boxplot
hist
pairs
です。こんなのも良いよってのがあったら教えてくださいm(_ _)m
※ライブラリの使用有無は問いません >>505
BESTパッケージのplotPost
histに加えて95%CIを表示してくれる。 >>505
出力をデータフレーム型にしたいときはskimrパッケージ あと、まだ試してないけsummarytoolsパッケージも面白そう >>501
良い方法が見つからなかったので !length(x)で空白ベクトル判定とした。
文字列を逆に並べる再帰呼び出しスクリプト
> reverse <- function(x){
+ if(!length(x)) return(NULL)
+ c(Recall(x[-1]),x[1])
+ }
> cat(reverse(LETTERS[1:26]))
Z Y X W V U T S R Q P O N M L K J I H G F E D C B A >>509
なぜ if (length(x) == 0) と書かないのか？
可読性って知ってるか？ >>510
文脈読めば
!length(x)で空白ベクトル
の流れ >>511
length は数値を返す関数なのに、なぜわざわざ論理演算をするんだ？
プログラムは、ここの議論を知らない人が読む可能性のほうが高いんだぞ。 >>512
正規分布が一様分布より大きい期待値の算出に
mean(rnorm(100)>runif(100))
って書かない？
俺は
sum(ifelse(rnorm(100)>runif(100),1,0))/100
って書いたりしないけど。 >>514
論理値を数値に置き換えて計算しているから
数値を論理値にしても別に違和感がない。
可読性は慣れの問題。 >>515
アンカー間違えとる？
だから、そんなの常套手段やん？ >>513
それとこれとは話が別だ。
logical は TRUE か FALSE の二値しかとらないから、それから数値への変換は自明。
多種の値をとる数値に論理演算をするのはいただけない。 >>515
is.empty.vector のような関数があるならそれでよいが、length を使っているのだから、それを勝手にベクトルが空かどうかの判断に使っているのはあなたであって、他人はそのようには考えない、といっているのだ。
自分だけしかこーどを見ないならべつに構わないが、こういうところに晒すのはよくない。
ましてや常套手段などと言って他人に教え込むのはやめてもらいたい。 >>518
常套手段というのはmeanの話じゃね？ whileの中は１でも２でもよくね？
n=0
while(1){
if(n>10) return(10)
n=n+1
}
n=0
while(2){
if(n>10) return(10)
n=n+1
} >509は
配列を逆順に並べる再帰呼出しのコード。
Rにはrevという関数がある。
そのコードみてみ！
lengthを真偽判定に使ってますがな。
> base:::rev.default
function (x)
if (length(x)) x[length(x):1L] else x >>518
>521みた？
if(length(x)!=0)になってないよ。 >>518
>こういうところ
誰も鵜呑みにしない便所の落書き。。。 !を引数が数値のときは0か否かを返す関数と理解すれば
if(!0) print(!1)の結果もサクッとわかる。 数値を非零かどうか論理値に変換して処理するかは
関数によるな。
if やwhileは変換処理しているけど
whichだとエラーになるな。
> which(!0)
[1] 1
> which(0)
Error in which(0) : argument to 'which' is not logical
> which(2)
Error in which(2) : argument to 'which' is not logical
> which(!3)
integer(0)
> !ってベクトル対応しているな
> which(!(-1:1))
[1] 2
> !(-2:2)
[1] FALSE FALSE TRUE FALSE FALSE !をnotではなくis.falseとかis.zero変換すれば( ・∀・)ｲｲ!! >>527
is.zero だと？
だったら最初から == 0 てよいじゃないか。 >>525
だいたい、if や while は関数じゃないし。
while (1) などと書くのはCに毒されているんじゃねーの？
Rなら while (TRUE) のほうが良いし、repeat というのもある。 >>528
スクリプトを読むときに脳内変換する癖をつけておけば
可読性が悪いと思わずにすむ。 >>518
>521みた？
if(length(x)!=0)になってないよ。 >>531
Rのライブラリ？の書き方なんてあてにならないよ。
そもそも、rev も 509 の reverse も演算の回数を必要最低限にするなら、条件は length(x)) >= 2 または > 1 だ。
まあ if (length(x)) のほうが速かったのかも知れないが、常にそうなるとは限らないのだから、早すぎる最適化は諸悪の根源だ。 可読性の次は速度かよ。
自分の流儀と違う { の使い方でも文句いいそう。 0,1を1000万個作って
!
!=
>
==
での真偽判定を各々１０回する時間を出してみた。
> x=rbinom(1e7,1,0.5)
> system.time(replicate(10,!x))
user system elapsed
1.27 0.33 1.59
> system.time(replicate(10,x!=0))
user system elapsed
2.39 0.36 2.75
> system.time(replicate(10,x>0))
user system elapsed
2.58 0.31 2.92
> system.time(replicate(10,x==1))
user system elapsed
2.60 0.36 2.99 !x と x!=0 で１０回やったが、
> f <- function(){
+ x=rbinom(1e7,1,0.5)
+ (a=system.time(!x)[3])
+ (b=system.time(x!=0)[3])
+ a<b
+ }
> (re=replicate(10,f()))
elapsed elapsed elapsed elapsed elapsed elapsed elapsed elapsed elapsed
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
elapsed
TRUE
結果は再現された。 lengthが絡むと
> x=1:1e8
> f1 = function(x) if(!length(x)) return(NULL) else x=x[-1]
> f2 = function(x) if(length(x)==0) return(NULL) else x=x[-1]
> system.time(f1(x))
user system elapsed
2.38 0.53 2.93
> system.time(f2(x))
user system elapsed
2.05 0.61 2.69
御指摘の通り、!の負け >>444
＞遅レスだが、Ubuntu 18.04LTSならデフォルトのIBus-mozcで日本語入力ができるそうだ。
＞詳しくはUbuntu Weekly Recipeの第527回を読んでくだされ
Ubuntu 18.04.1 LTSにアップグレードしたけどできない
いままでスクレイピングにつかってrvestのログインコマンドもエラー出るようになったし最悪だ。。。。 アップグレードってことは文字入力がFcitxのままになってない？Fcitxならパッチ要だよ。
Voyager 18.04LTS（xubuntu 18.04LTS）ではパッチなしでRStudioに日本語入力できた。 昨日、PCを新しくして最新のR入れたんだけど、
作業ディレクトリの固定化できなくね？
前まではショートカットのプロパティで作業フォルダ指定してたんだけど、
新しいR、何度やっても作業ディレクトリがデフォルトのドキュメントに戻りやがる
どうすればいいか、誰がご教示を。 >>539
プロパティでコマンドラインオプションを消せばよい。
--cd-ほにゃらら
というのが付いてるでしょ。 >>540
あ、なるほど
ありがとうございます、多謝 xtsの取り扱い方について質問失礼します。
以下のコードで、dataには2018/10/20 00:00〜2018/10/25 00:00の1時間刻みのaとbの値のxts型オブジェクトが入りますが、（aとbは時系列データのイメージです）
このdataの毎日09:00の値を添え字を使って取り出す方法はありますか？
（data["2018-10-21 00:00"]とすると10/21の00:00が取り出せるような感じで）
library(xts)
time <- seq(as.POSIXct("2018-10-20 00:00"), as.POSIXct("2018-10-25 00:00"), by="1 hour")
a <- rnorm(length(time))
b <- rnorm(length(time))
data <- xts(x = cbind(a,b), order.by = time)
思いついたのはfor文を使って以下のようにするのかなと思ったのですがもう少しいい方法があるだろうなと思いまして…
data_result <- data[10] # まずxts型の入れ物を作る(2018/10/20 09:00)
for (i in 1:4) {
data_result <- rbind(temp,data[10+24*i]) # 2日目以降はこれにくっつける
}
data_result 失礼、>>542の下のコードはこっちです
data_result <- data[10] # まずxts型の入れ物を作る(2018/10/20 09:00)
for (i in 1:4) {
data_result <- rbind(data_result, ,data[10+24*i]) # 2日目以降はこれにくっつける
}
data_result >>543
idx <- seq(10, length(time), by=24)
result <- data[idx]
じゃダメなのか？ >>544
ああ、これがいいですね！お恥ずかしい…
ありがとうございます 巨大数を扱えるというふれこみのRmpfrって正確じゃないな。
５０の階乗を計算させてみた
R with Rmpfr
> mpfr(factorial(50),1e5)
1 'mpfr' number of precision 100000 bits
30414093201713018969967457666435945132957882063457991132016803840
Haskell:
Prelude> product[1..50]
30414093201713378043612608166064768844377641568960512000000000000
Python:
import math
print(math.factorial(50))
30414093201713378043612608166064768844377641568960512000000000000
Wolfram:
https://www.wolframalpha.com/input/?i=50!
30414093201713378043612608166064768844377641568960512000000000000 R3.5.xから、Windows環境だと日本語パスまわりでエラーがでてどうにもならない
enc2native()しても駄目だし
とりあえずcairo_png()だけでも… 一度つかって日本語周りの処理でエラー出てきたんで未だに3.4.4浸かってる 統計の質問はこのスレでいいんだろうか？
ある医院に１時間あたり平均５人の患者が来院し、その人数の分布はポアソン分布にしたがうとする。
1時間あたりの平均診療人数は６人で、一人あたりの診療時間は指数分布に従うとする。
診察までの平均の待ち時間は何時間か？
このシミュレーション解はこれであってる？
N=1e6
λ=5
μ=6
sum(rpois(N,λ)*rexp(N,μ))/N
> sum(rpois(N,λ)*rexp(N,μ))/N
[1] 0.8331036
60かけて、分にすると
[1] 49.86565 >>306
grep 使って書いてみた。
# mhs(c(1,0,1,1,0,1,1,1)) : return 3
mhs = function(x){ # maximum head sequence
y=paste(x,collapse='')
str="1"
if(!grepl(str,y)) return(0)
else{
while(grepl(str,y)){
str=paste0(str,"1")
}
return(nchar(str)-1)
}
}
system.time(mean(replicate(1e4,mhs(sample(0:1,100,rep=T))>=5)))
# N(=100)回コインをｎ(=5回)以上続けて表がでるか？TRUE or FALSE
seqn<-function(n=5,N=100,p=0.5){
rn=rbinom(N,1,p)
count=0
for(i in 1:N){
if(rn[i] & count<n){
count=count+1
}
else{
if(count==n) {return(TRUE)}
else{
count=0
}
}
}
return(count==n)
}
system.time(mean(replicate(1e4,seqn())))
結果は、
> system.time(mean(replicate(1e4,mhs(sample(0:1,100,rep=T))>=5)))
user system elapsed
4.56 0.01 4.61
> system.time(mean(replicate(1e4,seqn())))
user system elapsed
1.05 0.00 1.07
逆に４倍時間がかかるようになった coinの表裏を1と0で表して、その累積和を取ったベクトルを用意して
そのベクトルの5個前とのdiffの要素の中に5が一回でも現れることが、一回でも5回連続で表が出たことがあることの必要十分条件
grepはどうしても時間がかかるしif文もできれば使いたくない
searchseq <- function(n=5,N=100,p=0.5,trial=10000){
result <- 0 # 表が5回以上出た回数を数える入れ物
for (i in 1:trial){
coin <- rbinom(N,1,p)
coincumsum <- cumsum(coin)
coindiff <- diff(coincumsum,5)
#diff(coincumsum,5)の要素に一個でも5があれば表が5回以上出たことがあったということ
#anyでT/Fにして、sumで0/1にする
result <- result+sum(any(coindiff==5))
}
return(result/trial)
}
結果もよさそう
> searchseq()
[1] 0.80793
> system.time(mean(replicate(10000,mhs(sample(0:1,100,rep=T))>=5)))
ユーザ システム 経過
0.92 0.00 0.93
> system.time(mean(replicate(10000,seqn())))
ユーザ システム 経過
0.29 0.00 0.30
> system.time(searchseq())
ユーザ システム 経過
0.21 0.00 0.20 >>553
レスありがとうございました。
rleとも比べてみました。
> system.time(mean(replicate(1e4,seqn())))
user system elapsed
4.290 0.000 4.377
>system.time(mean(replicate(1e4,max(rle(rbinom(100,1,0.5))$len)>=5)))
user system elapsed
3.620 0.000 3.742
> system.time(mean(replicate(1e4,mhs(rbinom(100,1,0.5)>=5))))
user system elapsed
2.390 0.000 2.435
> system.time(searchseq())
user system elapsed
1.880 0.000 1.988 >>555
greplの逆転はsampleをrbinomに変えたためでしょう。
> system.time(replicate(1e5,sample(0:1,100,rep=T)))
user system elapsed
7.200 0.180 7.591
> system.time(replicate(1e5,rbinom(100,1,0.5)))
user system elapsed
5.980 0.230 6.319 不勉強ながらrle関数って初めて知ったけど使いやすそうだな >>555
rleは0も数えているから間違っているんじゃ？ >>558
ご指摘ありがとうございました。
修正しました。
rle1 <- function(N=100,n=5,p=0.5){
r=rle(rbinom(N,1,p))
max(r$len[which(r$val==1)])
}
結果
> system.time(mean(replicate(1e4,max(rle1()>=5))))
user system elapsed
4.430 0.000 4.546
> system.time(mean(replicate(1e4,seqn())))
user system elapsed
4.490 0.010 4.609
> system.time(mean(replicate(1e4,mhs(rbinom(100,1,0.5))>=5)))
user system elapsed
7.440 0.000 7.656
> system.time(searchseq())
user system elapsed
1.950 0.000 2.066
> 無理矢理１行にして実行
system.time(mean(replicate(1e4,any(diff(cumsum(rbinom(100,1,0.5)),5)==5))))
user system elapsed
1.820 0.000 1.886
>
> 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.370 0.010 4.478 >>560
意味もなくforループ回してた上に毎回sum使って真偽値を数値に変換してたけど
replicate使って最後に一回だけmean取ると2.066→1.886で1割短くなるのね
他人のコード読むのは勉強になる >>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 >>563
センスないなあ。
rle01 <- function(x) {
r <- rle(x)
one <- r$values == 1
if (any(one)) max(r$length[one]) else 0
} >>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
> >>565
そんな特別な場合のことに対して高速化するのは愚の骨頂。
1が一つでもある場合にsumを呼ぶのは余計じゃないのか？
余計なwhichがあるくらいだから、まずは素直にやりたいこと、やるべきことを正しく書くようにしたら？ 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 >>567
分数で表すと
890788167367/9007199254740992
9.889735334744909e-05 JPXのデータ、ファイル形式csvを読み込もうとするとうまく行かないんですが
どんな引数をつければいいですか >>566
実はシミュレーションじゃなくて
漸化式からのプログラム解を分数表示するプログラムはpythonで作成済。
ここに置いた。
https://egg.2ch.net/test/read.cgi/hosp/1540905566/77 コインを100回ふったときの表連続の最大数が5であったときの
このコインの表がでる確率の期待値、モード比、信頼区間を求めるのが次のネタ。
unirootで算出できたけど
シミュレーションはどうすればいいのかアイデアが浮かばない。
MCMCで解決できるかなぁ？
これをシミュレーションで検証したい。
$Rscript main.r
lower mean mode upper
0.2487456 0.4469764 0.4589692 0.6386493
http://tpcg.io/asKRE9 1億回のコイントスで何回連続して表がでる確率が高いかRでやってみた。
# maximal sequential head probability at 10^8 coin flip
> y
[1] 2.2204460492503131e-16 2.2204460492503131e-16
[3] 8.8817841970012523e-16 1.5543122344752192e-15
[5] 3.5527136788005009e-15 6.8833827526759706e-15
[7] 1.4210854715202004e-14 2.8199664825478976e-14
[9] 5.6843418860808015e-14 1.1346479311669100e-13
[11] 2.2737367544323206e-13 4.5452530628153909e-13
[13] 9.0949470177292824e-13 1.8187673589409314e-12
[15] 3.6379788070917130e-12 7.2757355695785009e-12
[17] 1.4551915228366852e-11 2.9103608412128779e-11
[19] 5.8207660913467407e-11 1.1641509978232989e-10
[21] 6.6493359939245877e-06 2.5720460687386204e-03
[23] 4.8202266324911647e-02 1.7456547460031679e-01
[25] 2.4936031630254019e-01 2.1428293501123058e-01
[27] 1.4106434838399229e-01 8.1018980443629832e-02
[29] 4.3428433624081136e-02 2.2484450838189007e-02
25回が続くのが4回に1回あることになる。
pythonで25回以上と25回ちょうどになるのを計算させてみた。
その結果、
Over 25
6977459029519597/9007199254740992
= 0.7746535667951356
Just 25
2246038053550679/9007199254740992
= 0.24936031612362342
高速化を狙ってCに移植したら
100万回で暴走。
https://egg.5ch.net/test/read.cgi/hosp/1540905566/132 "マラソン大会の選手に１から順番に番号の書かれたゼッケンをつける。
用意されたゼッケンN(=100)枚以下の参加であった。
無作為に抽出したM(=5)人のゼッケン番号の最大値はMmax(=60)であった。
参加人数推定値の期待値とその95%信頼区間を求めよ"
decken <- function(M, Mmax, N, conf.level=0.95){
if(Mmax < M) return(0)
n=Mmax:N
pmf=choose(Mmax-1,M-1)/choose(n,M)
pdf=pmf/sum (pmf)
mean=sum(n*pdf)
upr=n[which(cumsum(pdf)>conf.level)[1]]
lwr=Mmax
c(lower=lwr,mean=mean,upper=upr)
}
decken(M=5,Mmax=60,N=100)
> decken(M=5,Mmax=60,N=100)
lower mean upper
60.0000 71.4885 93.0000
これをシミュレーションで確認したい。
# simulation
M=5 ; Mmax=60 ; N=100
sub <- function(M,Mmax,N){
n=sample(Mmax:N,1) # n : 参加人数n
m.max=max(sample(n,M)) # m.max : n人からM人選んだ最大番号
if(m.max==Mmax) return(n)
}
sim <- function(){
n=sub(M,Mmax,N)
while(is.null(n)){
n=sub(M,Mmax,N) # 最大番号が一致するまで繰り返す
}
return(n)
}
runner=replicate(1e4,sim())
summary(runner) ; hist(runner,freq=F,col="lightblue")
quantile(runner,prob=c(0,0.95))
cat(names(table(runner)[which.max(table(runner))]))
> summary(runner) ; hist(runner,freq=F,col="lightblue")
Min. 1st Qu. Median Mean 3rd Qu. Max.
60.00 63.00 68.00 71.43 77.00 100.00
> quantile(runner,prob=c(0,0.95))
0% 95%
60 93
> cat(names(table(runner)[which.max(table(runner))])) # 最頻値
60
結果は確認できたけど、もっと高速なシミュレーションアルゴリズムはあるだろうか？ 重複があるか否かを返す、anyDuplicatedという関数を知ったので総当たり比較と早いかどうか比べてみた。
覆面算を □/□ ＊ □/□ ＝ □□/□を解くの使ってみた。
a/b * c/d == ef/g (c>a)として
F = function(fun){
n=1:9
ans=NULL
for(a in n){
for(b in n){
for(c in a:9){
for(d in n){
for(e in n){
for(f in n){
for(g in n){
if(fun(a,b,c,d,e,f,g)){
ans=rbind(ans,c(a,b,c,d,e,f,g))}}}}}}}}
return(ans)
}
で虱潰しに判定
F1=function(a,b,c,d,e,f,g){ #全部の組み合わせが等しくないのを確認
(a/b)*(c/d)==(10*e+f)/g &
a!=b & a!=c & a!=d & a!=e & a!=f & a!=g &
b!=c & b!=d & b!=e & b!=f & b!=g & c!=d &
c!=e & c!=f & c!=g & d!=e & d!=f & d!=g & e!=f & e!=g & f!=g
}
F2=function(a,b,c,d,e,f,g){ # anyDuplicatedで重複なしを判定
(a/b)*(c/d)==(10*e+f)/g & !anyDuplicated(c(a,b,c,d,e,f,g))
}
> system.time(F(F1))
user system elapsed
52.56 0.25 53.38
> system.time(F(F2))
user system elapsed
113.78 0.11 115.81
anyDuplicatedでコードは短くなるが速さが犠牲になった。 ■ このスレッドは過去ログ倉庫に格納されています
