18/10/26 12:58:59.43 8LsI2NoU.net
昨日、PCを新しくして最新のR入れたんだけど、
作業ディレクトリの固定化できなくね？
前まではショートカットのプロパティで作業フォルダ指定してたんだけど、
新しいR、何度やっても作業ディレクトリがデフォルトのドキュメントに戻りやがる
どうすればいいか、誰がご教示を。
561:１３２人目の素数さん
18/10/26 13:04:50.23 I9GOHEF6.net
>>541
プロパティでコマンドラインオプションを消せばよい。
--cd-ほにゃらら
というのが付いてるでしょ。
562:１３２人目の素数さん
18/10/26 13:13:30.11 8LsI2NoU.net
>>542
あ、なるほど
ありがとうございます、多謝
563:１３２人目の素数さん
18/10/26 23:05:34.11 5yNEfvRx.net
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
564:１３２人目の素数さん
18/10/26 23:12:48.52 Xg1/ChR5.net
失礼、>>544の下のコードはこっちです
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
565:１３２人目の素数さん
18/10/26 23:38:29.60 oEWqNGH8.net
>>545
idx <- seq(10, length(time), by=24)
result <- data[idx]
じゃダメなのか？
566:１３２人目の素数さん
18/10/27 01:10:34.48 NQVW6zPX.net
>>546
ああ、これがいいですね！お恥ずかしい…
ありがとうございます
567:１３２人目の素数さん
18/11/01 22:44:27.08 xCdOvDq8.net
巨大数を扱えるというふれこみの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:
URLﾘﾝｸ(www.wolframalpha.com)
568: 30414093201713378043612608166064768844377641568960512000000000000
569:１３２人目の素数さん
18/11/02 23:29:01.53 fMme23mF.net
こまけぇこたぁいいんだよ!!
570:１３２人目の素数さん
18/11/03 09:32:49.16 NiwdVATo.net
R3.5.xから、Windows環境だと日本語パスまわりでエラーがでてどうにもならない
enc2native()しても駄目だし
とりあえずcairo_png()だけでも…
571:１３２人目の素数さん
18/11/03 10:31:12.42 pqZePX+I.net
一度つかって日本語周りの処理でエラー出てきたんで未だに3.4.4浸かってる
572:１３２人目の素数さん
18/11/04 21:55:16.95 lTCeMsqQ.net
統計の質問はこのスレでいいんだろうか？
ある医院に１時間あたり平均５人の患者が来院し、その人数の分布はポアソン分布にしたがうとする。
1時間あたりの平均診療人数は６人で、一人あたりの診療時間は指数分布に従うとする。
診察までの平均の待ち時間は何時間か？
このシミュレーション解はこれであってる？
N=1e6
λ=5
μ=6
sum(rpois(N,λ)*rexp(N,μ))/N
> sum(rpois(N,λ)*rexp(N,μ))/N
[1] 0.8331036
60かけて、分にすると
[1] 49.86565
573:１３２人目の素数さん
18/11/08 12:35:54.73 wKTjJ6Fa.net
>>307
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
逆に４倍時間がかかるようになった
574:１３２人目の素数さん
18/11/08 14:05:00.34 1h1GfW+d.net
今さらだけどそれってrle使っちゃだめなの？
575:１３２人目の素数さん
18/11/08 15:13:38.19 UNhl34kg.net
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
576:１３２人目の素数さん
18/11/08 22:10:55.78 I6htNxzg.net
あっ、コード中の5はnに置き換えてね
577:１３２人目の素数さん
18/11/08 23:32:48.68 /ZlhhGjJ.net
>>555
レスありがとうございました。
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
578:１３２人目の素数さん
18/11/08 23:35:10.38 /ZlhhGjJ.net
>>557
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
579:１３２人目の素数さん
18/11/09 00:20:12.33 BZRFS9lT.net
不勉強ながらrle関数って初めて知ったけど使いやすそうだな
580:１３２人目の素数さん
18/11/09 00:36:46.75 Ui6BpICy.net
>>557
rleは0も数えているから間違っているんじゃ？
581:１３２人目の素数さん
18/11/09 02:01:36.95 Qla0VxTD.net
>>560
ご指摘ありがとうございました。
修正しました。
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
>
582:１３２人目の素数さん
18/11/09 07:25:00.97 Qla0VxTD.net
無理矢理１行にして実行
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
583:１３２人目の素数さん
18/11/09 07:59:47.47 rijPDFSR.net
>>562
意味もなくforループ回してた上に毎回sum使って真偽値を数値に変換してたけど
replicate使って最後に一回だけmean取ると2.066→1.886で1割短くなるのね
他人のコード読むのは勉強になる
584:１３２人目の素数さん
18/11/09 09:06:56.99 Ui6BpICy.net
>>561
whichいらねーよ。
585:１３２人目の素数さん
18/11/09 09:44:07.41 5AnUTlVm.net
>>561
全部が０のとき、エラーになるので修正
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
586:１３２人目の素数さん
18/11/09 10:10:01.39 Ui6BpICy.net
>>565
センスないなあ。
rle01 <- function(x) {
r <- rle(x)
one <- r$values == 1
if (any(one)) max(r$length[one]) else 0
}
587:１３２人目の素数さん
18/11/09 10:40:31.84 5AnUTlVm.net
>>566
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
>
588:１３２人目の素数さん
18/11/09 10:52:23.39 Ui6BpICy.net
>>567
そんな特別な場合のことに対して高速化するのは愚の骨頂。
1が一つでもある場合にsumを呼ぶのは余計じゃないのか？
余計なwhichがあるくらいだから、まずは素直にやりたいこと、やるべきことを正しく書くようにしたら？
589:１３２人目の素数さん
18/11/09 11:01:40.74 5AnUTlVm.net
1000本に1本あたる宝クジを100本買って続けて2本あたる確率のシミュレーション解の算出時間比較。
590:確率の理論値は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
591:１３２人目の素数さん
18/11/09 15:22:38.24 5AnUTlVm.net
>>569
分数で表すと
890788167367/9007199254740992
9.889735334744909e-05
592:１３２人目の素数さん
18/11/09 19:43:44.34 0M0agNOf.net
JPXのデータ、ファイル形式csvを読み込もうとするとうまく行かないんですが
どんな引数をつければいいですか
593:１３２人目の素数さん
18/11/10 10:36:49.41 QJ6NJqU7.net
>>568
実はシミュレーションじゃなくて
漸化式からのプログラム解を分数表示するプログラムはpythonで作成済。
ここに置いた。
URLﾘﾝｸ(egg.2ch.net)
594:１３２人目の素数さん
18/11/10 11:20:41.80 QJ6NJqU7.net
コインを100回ふったときの表連続の最大数が5であったときの
このコインの表がでる確率の期待値、モード比、信頼区間を求めるのが次のネタ。
unirootで算出できたけど
シミュレーションはどうすればいいのかアイデアが浮かばない。
MCMCで解決できるかなぁ？
これをシミュレーションで検証したい。
$Rscript main.r
lower mean mode upper
0.2487456 0.4469764 0.4589692 0.6386493
URLﾘﾝｸ(tpcg.io)
595:１３２人目の素数さん
18/11/11 20:11:21.81 ODPKEEGK.net
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万回で暴走。
ｽﾚﾘﾝｸ(hosp板:132番)
596:１３２人目の素数さん
18/11/12 21:04:57.19 boi8bvdM.net
"マラソン大会の選手に１から順番に番号の書かれたゼッケンをつける。
用意されたゼッケンN(=100)枚以下の参加であった。
無作為に抽出したM(=5)人のゼッケン番号の最大値はMmax(=60)であった。
参加人数推定値の期待値とその95%
597:信頼区間を求めよ" 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 結果は確認できたけど、もっと高速なシミュレーションアルゴリズムはあるだろうか？
598:１３２人目の素数さん
18/11/16 13:44:59.80 U19cHKqd.net
重複があるか否かを返す、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でコードは短くなるが速さが犠牲になった。
599:１３２人目の素数さん
18/11/16 16:13:49.03 zcEc4+wx.net
カルマンフィールタ
600:１３２人目の素数さん
18/11/18 12:56:58.01 tDkQkz0a.net
初歩的な質問で恐縮ですが、
truehistで出したヒストグラムの情報をテキストファイルとして保存する方法、教えて頂けませんか？
601:１３２人目の素数さん
18/11/18 16:15:11.61 hHA2tyUi.net
>>578
truehist関数はhist関数と違って戻り値を返さないから無理なのでは
602:１３２人目の素数さん
18/11/18 17:12:37.03 HQweHJGH.net
>>578
MASS::truehist
でソースを覗いたら最後はinvisible()になっていた。
ここを
return(list(breaks=breaks,h=h,nbins=nbins,xlab=xlab))
にしたら
$`breaks`
[1] -0.001 0.999 1.999 2.999 3.999 4.999 5.999 6.999 7.999 8.999 9.999 10.999 11.999
[14] 12.999 13.999
$h
[1] 1
$nbins
[1] 17
$xlab
[1] "x"
で出力されるけど、
どのパラメータがあればヒストグラムが再現できるのかは不勉強に
603:てわからない。
604:１３２人目の素数さん
18/11/18 18:13:02.48 hHA2tyUi.net
関数いじるならbreaksとestを出力させればよい
605:１３２人目の素数さん
18/11/18 20:37:19.22 HQweHJGH.net
>>581
レスありがとうございます。
すると
truehist のソースの最後の
invisible()
を
return(list(breaks=breaks,est=est))
として、breakの中点を横軸、estを縦軸にするとヒストグラムが再現できるわけですね。
606:１３２人目の素数さん
18/11/18 20:40:50.06 HQweHJGH.net
invisible() → return(list(breaks=breaks,est=est))
に改造したソースをtruehist0として
midpoint <- function(x){ # c(1,2,3,4) -> c(1.5,2.5,3.5)
n=length(x)
mpt=numeric(n-1)
for(i in 1:(n-1)){
mpt[i]=mean(x[i],x[i+1])
}
return(mpt)
}
x=rnorm(10000)
(data=truehist0(x))
with(data, plot(midpoint(breaks),est,type='h',lwd=5,col='cyan'))
で元のヒストグラムが再現できた。
607:１３２人目の素数さん
18/11/18 20:52:08.49 HQweHJGH.net
graphics:::hist.default
でhistのソースを表示させてみた。
r <- structure(list(breaks = breaks, counts = counts, density = dens,
mids = mids, xname = xname, equidist = equidist), class = "histogram")
if (plot) {
plot(r, freq = freq1, col = col, border = border, angle = angle,
density = density, main = main, xlim = xlim, ylim = ylim,
xlab = xlab, ylab = ylab, axes = axes, labels = labels,
...)
invisible(r)
}
histでのcountsが
truehistではestと呼ばれているようだ。
estimateの略かな？
608:１３２人目の素数さん
18/11/20 18:37:20.95 kwMT/Yc9.net
どいつもこいつもナイル川で説明しやがって
データの操作が一番むずいんだよ！
609:１３２人目の素数さん
18/11/21 12:43:45.15 4Qrr7yQM.net
あるデータ群に対して、確率密度関数のパラメータをフィッティングさせる方法ってないですか？
ちなみに、フィッティングさせたいのはレブィフライト確率密度関数です。
610:１３２人目の素数さん
18/11/21 17:15:18.99 /i6dFnq8.net
普通に最尤推定できないっけ？
最小二乗法でもできた気がする
611:１３２人目の素数さん
18/11/21 21:23:46.28 kb0MtFao.net
>>586
MASSのfitdistとVGAMのlevyを使うとなんとかなるかも。
やったことないけど。
612:１３２人目の素数さん
18/11/21 21:59:06.34 kb0MtFao.net
>>586-587
とりあえず、最小二乗法でやってみた。
ガンマ分布の乱数を近似してみた。
dlevy <- function (x,m,c) sqrt(c/2/pi)*exp(-c/2/(x-m))/(x-m)^3/2
set.seed(123)
dat=rgamma(1e3,1) ; hist(dat,freq=F)
x=density(dat)$x ; y=density(dat)$y
lines(x,y)
f<-function(mc){
m=mc[1];c=mc[2]
sum((y-dlevy(x,m,c))^2)
}
(mc=optim(c(0,1),f, method='N')$par)
curve(dlevy(x,mc[1],mc[2]),add=T,col=2)
613:１３２人目の素数さん
18/11/21 23:39:04.38 8W1KB4Wk.net
>>588
fitdistは対応できる分布が限定されているから
のソースを改造しないと無理だな。
614:１３２人目の素数さん
18/11/22 19:54:16.79 bptUComJ.net
ソースが長かったので
sink('print.out')
print(MASS::fitdistr)
sink()
でprint.outに出力してみた。
これもoptimを使っているようだが、最小二乗法なのかどうなのかわからなかった。
615:１３２人目の素数さん
18/11/23 15:40:10.08 lNNKNPJr.net
>>552
統計というより待ち行列理論だね
50分が答えになるから合ってそう
616:１３２人目の素数さん
18/11/23 18:16:00.74 rwOyVQ2V.net
>>592
レスありがとうございます。
λ=5
μ=6
N=1e6
sum(rpois(N,λ)*rexp(N,μ))/N
> sum(rpois(N,λ)*rexp(N,μ))/N
[1] 0.833631
ρ=λ/μ
ρ/(1-�
617:ﾏ) ρ/(1-ρ)*1/μ > ρ/(1-ρ)*1/μ [1] 0.8333333 なのでいいのだろうとは思っていたのですが、時系列のシミュレーションは自信がありませんでした。
618:１３２人目の素数さん
18/11/23 21:47:20.33 rwOyVQ2V.net
>>593
>552の設定で12分毎に患者が来院したら、待ち時間は全員０だと思うのだが
ρ/(1-ρ) の公式って正しいんだろうかな？
619:１３２人目の素数さん
18/11/24 00:40:11.23 +sZNyHbC.net
50分待ちだと常時待合室で5人は待ってることになるな
620:１３２人目の素数さん
18/11/27 11:01:28.21 V3tvhpxu.net
>>594
自己レス
公式は定常状態に達したときという前提での計算なんだな。
MMS = function(n, lamda=5,mu=6,s=1){
rho=lamda/mu
sig=0
for(i in 0:s) sig=sig+rho^i/factorial(i)
p0=1/( sig + rho^(s+1)/factorial(s)/(s-rho) )
ifelse(n >= s, rho^n/factorial(s)/s^(n-s)*p0, rho^n/factorial(n)*p0)
}
E=0
for(i in 0:1000) E=E+i*MMS(i)
> E
[1] 5
621:１３２人目の素数さん
18/11/27 11:12:11.41 +K2fEP5F.net
みやぞん分布の話？
622:１３２人目の素数さん
18/11/27 16:44:04.82 V3tvhpxu.net
>>597
こういう類の待ち時間の話。
ある医院では、患者が平均１０分間隔でポアソン分布にしたがって訪ねてくることがわかった。
医者は１人であり、１人の患者の診療にかかる時間は平均８分の指数分布であった。
「平均待ち時間」を5分以下にするには同じ診察効率の医師が何人に必要か？
その最小人数で「平均待ち時間」を5分以下に保って診療するには１時間に何人まで受付可能か？
公式に当て嵌めれば解けるのだけど
どうやってシミュレーションすればいいのか思い浮かばない。
コイントスやサイコロだとシミュレーションは容易なんだが。
623:１３２人目の素数さん
18/11/28 14:59:50.66 r8zTzMor.net
# シミュレーションしたみたが、結果が合致しない（特定のseedでは合致したけど）
# ある医院に１時間あたり平均５人の患者が来院し、その人数の分布はポアソン分布にしたがうとする。
# 1時間あたりの平均診療人数は６人で、一人あたりの診療時間は指数分布に従うとする。
# 診察までの平均の待ち時間は何時間か？
MM1sim <- function(n=40,lambda=5/60,mu=6/60,seed=FALSE,Print=TRUE){
# service starc clock time(ssct) since 9:00
ssct=numeric(n)
# waiting time(w8)
w8=numeric(n)
# service end clock time(sect)
sect=numeric(n)
# arrival clock time(act)
if(seed) set.seed(1234) ;
act=round(cumsum(rexp(n,lambda)))
# duration of service(ds)
if(seed) set.seed(5678) ;
ds=round(rexp(n,mu))
# simulation assuming service starts at 9:00
head(act) # act : arrival clock time
head(ds) # ds : duration of service
# initial values
ssct[1]=act[1] # 9:15 service start clock time for 1st guest
sect[1]=act[1]+ds[1] # 9:25 sevice end clock time for 1st guest
w8[1]=0
for(i in 2:n){
w8[i]=max(sect[i-1]-act[i],0)
ssct[i]=max(sect[i-1],act[i])
sect[i]=ssct[i]+ds[i]
}
if(Print){
print(summary(w8))
hist(w8,freq=FALSE,col="lightblue",main="")
}
invisible(w8)
}
w8m=replicate(1e3,mean(MM1sim(P=F)))
summary(w8m)
624:１３２人目の素数さん
18/11/28 15:00:27.16 r8zTzMor.net
途中の計算サンプル
# simulation step by step
#
# act[2] # 9:16 arrival clock time of 2nd
# max(sect[1]-act[2],0) # 9:25-9:16 vs 0 = ?sevice for 1st ends b4 2nd arrival
# w8[2]=max(sect[1]-act[2],0) # 9 min : w8ing time of 2nd
# ssct[2]=max(sect[1],act[2]) # 9:25 vs 9:16 = service start clock time for 2nd
# sect[2]=ssct[2]+ds[2] # 9:25 + 8 = 9:33 service end clock time for 2nd
#
# act[3] # 9:17 arrival clock time of 3rd
# max(sect[2]-act[3],0) # 9:33 - 9:17 vs 0 = ?serivce for 2nd ends b4 3rd arrival?
# w8[3]=max(sect[2]-act[3],0) # 16 min : w8ting time of 3rd
# ssct[3]=max(sect[2],act[3]) # 9:33 vs 9:17 = service start clock time for 3rd
# sect[3]=ssct[3]+ds[3] # 9:33 + 11 = 9:44 service end clock time for 3rd
#
625:１３２人目の素数さん
18/11/28 19:28:45.39 r8zTzMor.net
>>59
626:9 患者来院数４０人程度では定常状態に達しないということみたいだな。 １０万人での待ち時間の分布（理論値５０に近い） Min. 1st Qu. Median Mean 3rd Qu. Max. 0.000 5.532 28.619 47.540 67.370 528.511 １００人での待ち時間の分布 （シミュレーションの度にばらつく） Min. 1st Qu. Median Mean 3rd Qu. Max. 0.000 6.958 16.564 22.200 31.652 109.204 ４０人来院を１万回繰り返した平均の分布 Min. 1st Qu. Median Mean 3rd Qu. Max. 1.304 11.764 19.357 25.449 32.536 177.913 ここに上げておきました。 http://tpcg.io/7psmrQ 結局のところ、待ち時間行列の理論を個人医院に適応しても 定常状態（待ち行列の長さが一定）に達していないと実用性がなさそうだな。 シミュレーションはなんとなく機能していた感じ。 ここに上げておきました。 https://www.tutorialspoint.com/tpcg.php?p=7psmrQ http://tpcg.io/7psmrQ
627:１３２人目の素数さん
18/11/29 18:50:23.24 GEvMs+K3.net
来院患者数と待ち時間シミュレーションの結果をグラフにしてみた。
個人医院で待ち行列の長さが一定になるほどの受診があるとは考えがたいからこっちの方が現実に即しているんじゃなかろうか
と思うが、解析解はどんなふうになるのか想像もつかん。
URLﾘﾝｸ(i.imgur.com)
628:１３２人目の素数さん
18/12/01 18:03:07.99 gm1XmpT3.net
初歩的な質問ですいません
truehistで軸を対数軸にして表示する方法って分かりますか？
629:１３２人目の素数さん
18/12/02 00:37:54.68 XyqxkFs6.net
truehistだと以前の質問でお礼すらしなかった人か
630:１３２人目の素数さん
18/12/02 07:13:03.31 ailhRsax.net
>>604
答えたの俺だがソース改造できたのかなぁ。
631:１３２人目の素数さん
18/12/02 08:28:26.92 gEGEypwn.net
histで縦軸を対数表示させるならこんな感じかな。
with(hist(c(rnorm(1e6),rnorm(1e5,5,0.5))),plot(mids,log10(counts),type='h',col=4,lwd=10))
truehistだとソース改造すればいい。
632:１３２人目の素数さん
18/12/02 09:43:09.85 bQp1mbWe.net
>>606
plot(log='y')でもいいな。
633:１３２人目の素数さん
18/12/02 13:24:03.04 8De20Fh0.net
histグラムぽく対数表示
例
dat=hist(c(rnorm(1e6),rnorm(1e5,5,0.5)))
attach(dat)
plot(breaks[-1],counts,type='s',log='y',ylim=range(counts))
segments(x0=breaks[-1],y=min(counts),y1=counts)
segments(x0=breaks[1],y=min(counts),y1=counts[1])
634:１３２人目の素数さん
18/12/04 00:31:54.32 4UCEIb49.net
すみませんがアドバイスお願いします。
heatmap.2でDendrogramつきヒートマップを描きたいのですが、カラムの並びを任意に変えたいです。
Dendrogramの細かい並びは変えないように、大きなクラスタの並びを変えたいです。
例えば、（１，２，３）（４，５，６）（７，８，９）とあるのを
（４，５，６）（７，８，９）（１，２，３）とならび変えるのが目的です。
このとき、4,5,6と7，8，9は近いクラスタを形成します。なので、樹形図は崩れないように書き換えられると思っています。
URLﾘﾝｸ(www.biostars.org)
上記サイトをみると、as.dendrogramのアウトプットをreorderで並び替えてColvに入れるようですが、
うまくいきませんでした。並びは全く変わっていません。
どなたか教えていただけますか。情報に漏れがありましたらご指摘ください。
環境は以下の通りです。
R 3.4.0
Rstudio 1.0.143
Win10になります。
635:１３２人目の素数さん
18/12/04 10:04:43.34 QKKYvADK.net
頓珍漢な答かもしれない
636:。 並べ替えだけなら x=rbind(c(1,2,3),c(4,5,6),c(7,8,9)) x[c(2,3,1),]
637:１３２人目の素数さん
18/12/04 10:06:24.27 QKKYvADK.net
listなら
x=list(c(1,2,3),c(4,5,6),c(7,8,9))
x[c(2,3,1)]
638:１３２人目の素数さん
18/12/04 18:58:07.88 7f8uMrnq.net
初歩的な話なんだが、一様分布の分散は無限大かと思っていたら区間[a,b]で(a-b)^2/12とのこと。
Wolframで計算したら確かにそうなった。
URLﾘﾝｸ(www.wolframalpha.com)(x-(a%2Bb)%2F2)%5E2%2F(b-a),from+a+to+b
バスの到着時間が平均１０分の指数分布に従うときにランダムにバス停に行ったときの平均待ち時間は１０分。
バスの到着時間が平均１０分の一様分布に従うときにランダムにバス停に行ったときの平均待ち時間は６分４０秒。
バスがきちんと１０分毎に到着するときはランダムにバス停に行ったときの平均待ち時間は５分。
乱数発生させて公式でのシミュレーション
> d2w8 <- function(x){# w8=E[X2]/2E[X]=(V[X]+E[X]^2)/2E[X]
+ c(mean=mean(x),var=var(x),w8=mean(x^2)/mean(x)/2)
+ }
> N=1e6
> d2w8(rexp(N,1/10)) # exp average:10
mean var w8
10.02477 100.67652 10.03377
> d2w8(runif(N,0,20)) # unif average:10
mean var w8
9.997470 33.325065 6.665408
> d2w8(rep(10,N)) # regular interval 10
mean var w8
10 0 5
639:１３２人目の素数さん
18/12/06 14:39:28.97 P/rPOK1I.net
NULLのときってどうしてこういう仕様なんだろ？
プログラムしていたら、これに気づかないのがバグの原因だったｗ
> any(NULL)
[1] FALSE
> all(NULL)
[1] TRUE
640:１３２人目の素数さん
18/12/06 20:13:51.63 P/rPOK1I.net
= と <-で微妙に動作が違うな。
> switch (3,
+ x =1,
+ x =2,
+ x =3
+ )
[1] 3
> switch (2,
+ x <- 1,
+ x <- 2,
+ x <- 3
+ )
641:１３２人目の素数さん
18/12/06 21:58:17.80 EMJ7DN40.net
>>614
そのふたつは意味が異なる
それぞれちゃんと命令どおりの挙動
642:１３２人目の素数さん
18/12/06 22:22:18.13 P/rPOK1I.net
どちらが見やすいかという問題かな。
> rm(x)
> x=switch(1,x =1)
> x
[1] 1
> switch(1,x<-1)
> x
[1] 1
643:１３２人目の素数さん
18/12/06 23:13:35.91 EMJ7DN40.net
>>616
その二つは異なるアルゴリズムでたまたま結果が同じになっているだけ
=と<-の違いはRやるなら理解しておいたほうが良い
644:１３２人目の素数さん
18/12/07 07:46:59.34 H6LI4wTx.net
>>617
scopeが違うってことですね。
645:１３２人目の素数さん
18/12/07 08:01:49.32 H6LI4wTx.net
0^x =0
x^0=1
0^0=1とした方が辻褄が合うことが多いけど
Rのこの仕様には何のメリットがあるんだろ？
> any(NULL)
[1] FALSE
> all(NULL)
[1] TRUE
>
646:１３２人目の素数さん
18/12/07 16:08:57.61 3jKAkgsb.net
奥村さんみたいなこと言わんといて
神経質すぎ
647:１３２人目の素数さん
18/12/23 11:38:19.16 CnP6hLfL.net
>>120
平成29年の簡易生命表から
f=c(179,28,19,13,9,7,6,5,5,4,4,4,5,7,9,10,11,12,14,16,18,19,19,20,21,22,23,24,25,27,28,30,
31,34,37,39,41,43,46,51,57,63,70,77,83,91,99,109,119,130,142,155,167,178,190,202,216,
233,249,265,282,302,326,354,387,420,457,502,549,598,648,700,761,836,927,1026,1142,1284,
1455,1651,1862,2089,2341,2625,2934,3264,3598,3923,4233,4508,4740,4893,4973,5007,4999,
4729,4314,3797,3222,2634,2071,1566,1135,788,523,761)
m=c(191,31,21,13,10,8,8,8,7,7,7,8,9,11,14,17,21,26,32,37,42,46,49,50,50,50,49,49,50,52,55,
648:58,60,63,65,67,71,76,82,89,97,104,112,122,134,148,165,183,203,224,246,268,294,324,357, 391,425,461,502,549,601,659,722,792,872,958,1052,1147,1239,1331,1433,1546,1663,1783, 1905,2026,2167,2333,2532,2750,2973,3195,3414,3630,3827,4000,4133,4200,4194,4104,3916, 3681,3388,3046,2669,2272,1875,1494,1145,841,589,392,246,144,79,71) LE <-function(ndx,Y,N0=10^5){ # life expectancy n=length(ndx) lx=numeric(n) lx[1]=N0 for(i in 1:(n-1)) lx[i+1] <- lx[i] - ndx[i] nqx=ndx/lx nLx=numeric(n) for(i in 1:n) nLx[i] <- mean(c(lx[i],lx[i+1])) nLx[n]=0 Tx=rev(cumsum(rev(nLx))) le=Tx/lx return(round(le[Y+1],1)) } LE(m,65) LE(m,61)
649:１３２人目の素数さん
18/12/28 00:54:31.27 fO3GkB5Y.net
pushってありますか?
例えばベクトルに値を追加すると
先頭が消えて後尾に新しい値をついかして要素数を一定に保つような
一定要素数以下の平均を求めたいのでそういうのが簡単に実現できる方法あればなおよいです
650:１３２人目の素数さん
18/12/28 01:17:50.69 fO3GkB5Y.net
ベクトルをrev()して先頭を[1:50]とかでとりだしてman()すればいいとわかりました
ほかにもっと簡単な方法アレばお湿気てください
651:１３２人目の素数さん
18/12/28 07:19:51.35 1iwNDVav.net
append(x[-1],y)
652:１３２人目の素数さん
18/12/28 08:10:10.49 4Fwuxgb+.net
>>623
mean(tail(x, 50))
653:１３２人目の素数さん
18/12/29 08:09:54.89 Pk3SjXcs.net
組み合わせると
f=function(x,y,n=50) mean(tail(append(x,y),n))
654:１３２人目の素数さん
18/12/30 10:52:59.30 3dqEgqR+.net
>>626
これはおかしい
655:１３２人目の素数さん
19/01/07 02:16:23.34 OfD+CJ7t.net
時系列データをplot()したときに縦線をabline()でいれたいのだがどうすればいいかよくわからないです
具体的には
timeStr <- "2018-01-07 01:00"
dateTime <- strptime(timeStr, format="%Y-%m-%d %H:%M") #"POSIXlt" クラスオブジェクトに変換
として時間データに変換したものをx軸としてプロットしたものです
たとえば
plot(x=dateTime, y=1)
abline(v=?????)
として縦線を追加したいのですが
656:１３２人目の素数さん
19/01/07 02:50:32.42 OfD+CJ7t.net
v=as.POSIXct("2019-01-06 01:00")
みたいな感じにすれば解決しました
657:１３２人目の素数さん
19/01/07 18:42:57.56 oKAYsRfx.net
>>278
かなり遅れすだけど
formals(cor)
658:１３２人目の素数さん
19/01/07 19:36:20.34 oKAYsRfx.net
>>613
>>619
俺の予想
判定するときにNULLは強制的に論理値に変換される
変換されると logical(0) になる
logical(0) は空のベクトル
からのベクトルが渡されて中身をチェックしていく
アルゴリズムとして早く処理するためには
any()はTRUE探しにいって、一個でもTRUEがみつかればその時点でTRUEをかえす
all()はFALSEを探しに行って、一個でもFALSEがみつかればその時点でFALSEをかえす
もちろん最後までみない
空のベクトルが渡されたのでTRUEもFALSEもみつからない、となると
any()ではFALSEとなり
all()ではTRUEとなる
これで辻褄はあう
ちなみに
any(c())
all(c())
でも同じ結果が出る
659:１３２人目の素数さん
19/01/08 12:14:42.13 LiVTLUAx.net
1行のテキストデータを最終的に数値とテキストの混在するN行M列のデータフレームにしたいのですが、なかなかうまく出来ません。
1行データの構造の設計からデータフレームの変換までどうすればシンプルに実現できるか助言ください。
たとえば
1 a
2 b
というようなデータを
1-a,2-b
というような一行のテキストデータから初めてテーブル構造にするというような形です
x <- "1-a,2-b"
y <- str_split(x, ",")
と試しにやってみたのですがyが行単位のベクトルになるだけでここからどうデータフレームにすればよいかわかりません
660:１３２人目の素数さん
19/01/08 12:50:44.29 qzUBBMdZ.net
>>632
何をしたいのか、まったく理解できない。最低限、N, Mが何なのか説明したらどう？
661:１３２人目の素数さん
19/01/08 13:31:30.92 By0HLtnN.net
>>632
文字列 1-a,2-b を 数値とテキストのデータフレームにしたいという意味と解した。
x="1-a,2-b"
y=strsplit(x,",")
z=unlist(y)
w=NULL
for(i in 1:length(z)){
w=rbind(w,unlist(strsplit(z[i],"-")))
}
data.frame(NUM=as.numeric(w[,1]),TEXT=w[,2])
662:１３２人目の素数さん
19/01/08 13:33:34.43 By0HLtnN.net
実行結果
> x="1-a,2-b"
> y=strsplit(x,",")
> z=unlist(y)
> w=NULL
> for(i in 1:length(z)){
+ w=rbind(w,unlist(strsplit(z[i],"-")))
+ }
> data.frame(NUM=as.numeric(w[,1]),TEXT=w[,2])
NUM TEXT
1 1 a
2 2 b
>
663:１３２人目の素数さん
19/01/08 13:34:00.52 LiVTLUAx.net
>>633
NMは任意の数字です
したいことは
データを一行に記録してそれを
テーブル構造にすることです。
これだけです。
一旦ファイルに保存してread.csvなどにすればよいのでしょうが
いちおう直接テキストデータをコピペしてからということにしたいのです。
なぜこんなことをするかと言うと
TamperMonkeyという拡張機能でJavaScriptで使ってウェブ上のデータを収集しているのですが、
localStrageというクッキーの拡張版のような機能をつかってデータをクライアントに保存するときに基本的にテキストデータベースでの保存になので
一行のデータとして後ろにどんどんデータを追加していくのが一番単純な処理に成るからです。
そのlocalStorageに保存されたデータはコピペで取り出すしかないので
それをRで処理する時に一行のテキストデータから始めないといけないのです。
664:１３２人目の素数さん
19/01/08 13:56:11.34 qzUBBMdZ.net
>>636
, を行のデリミタ、- を列のデリミタとして、一行の文字列をデータフレームにするというのであれば、
s に文字列が入っているとして、
r <- unlist(strsplit(s, ","))
d <- lapply(r, function(x) unlist(strsplit(x, "-")))
as.data.frame(do.call(rbind, d))
列数が行によって異なるときにどうなるかは知らん。
665:１３２人目の素数さん
19/01/08 14:57:46.44 By0HLtnN.net
>>631
内部動作の考証ありがとうございます。
未だにこれは理解できません
> logical(NULL)
Error in logical(NULL) : invalid 'length' argument
> logical(0)
logical(0)
> logical(1)
[1] FALSE
666:１３２人目の素数さん
19/01/08 19:53:14.03 CxwlQqo4.net
>>632
> txt <- "1-a,2-b,3-c"
> read.table(text = gsub(',', '¥n', txt), sep = '-')
V1 V2
1 1 a
2 2 b
3 3 c
こんな感じか？
667:１３２人目の素数さん
19/01/08 20:04:46.51 CxwlQqo4.net
>>638
自分の解釈では、
NULL は「無」なのでエラー（引数はベクトルの要素数を必ず与えなければいけない）
0のとき、0個という指定なので、0個の要素をもつベクトル
1のとき、1個という指定なので、1個の要素を持つベクトル（規定値はFALSE）
ちなみに、
> logical(2)
[1] FALSE FALSE
2のとき、2個という指定なので、２個の要素を持つベクトル（規定値はFALSE）
668:１３２人目の素数さん
19/01/08 20:41:30.30 gJU+fDZy.net
>>640
解説ありがとうございました。
669:１３２人目の素数さん
19/01/08 21:41:07.74 LiVTLUAx.net
皆さんありがとうございます。
>>639 ファイル入れなくてもできるんですね
671:１３２人目の素数さん
19/01/08 22:24:00.83 LiVTLUAx.net
最古の元号って大化なん?
一巡回って大化2でいいよ
あとは干支みたいに回せばいい
672:１３２人目の素数さん
19/01/08 22:43:03.83 9qwcPkmK.net
三文字やめれ
673:１３２人目の素数さん
19/01/08 23:33:49.53 /Jd5B1BR.net
じゃあ太蟹(たいかに)で
674:１３２人目の素数さん
19/01/08 23:44:56.11 LiVTLUAx.net
>>643は誤爆ですw
レス付くと思わなかった。。。
675:１３２人目の素数さん
19/01/08 23:47:19.60 LiVTLUAx.net
>>641
⇣の結果をすべて即答できるようになれば合格だと思います…
any(c(TRUE,NA))
any(c(FALSE,NA))
any(c(FALSE,NA,TRUE))
any(c(NA, TRUE))
any(c(NA,FALSE,TRUE))
all(c(TRUE,NA))
all(c(FALSE,NA))
all(c(FALSE,NA,TRUE))
all(c(NA, TRUE))
all(c(NA,FALSE,TRUE))
TRUE | NA
FALSE | NA
FALSE | NA | TRUE
NA | TRUE
NA | FALSE
NA | FALSE | TRUE
TRUE & NA
FALSE & NA
FALSE & NA & TRUE
NA & TRUE
NA & FALSE
NA & FALSE & TRUE
676:１３２人目の素数さん
19/01/09 00:23:21.60 aVxwJ5mP.net
鯛蟹で
677:１３２人目の素数さん
19/01/09 02:19:38.70 tVKwPfCD.net
改元は改源で
678:１３２人目の素数さん
19/01/09 06:18:55.77 DWUqFfaE.net
源義光とか
679:１３２人目の素数さん
19/01/10 19:23:43.99 a/KDy24T.net
ggplot2で凡例をまとめる事ってできないでしょうか。
例えば、下記のコードでは線と点で別々の凡例になります。
線と点のスタイルを合わせて１つの凡例にしたいのですがどうすればいいでしょうか。
geom_line(mapping=aes(colour=Conditions),alpha=0.6)
+ geom_point(mapping=aes(shape=Conditions,colour=Conditions),alpha=0.8)
+ scale_shape_manual(values = 1:3)
680:649
19/01/10 20:22:57.38 a/KDy24T.net
すみません自己解決しました。
このコードではなくもっと後のthemeで凡例のスタイルを変えるときに余計なことをしていました。
681:１３２人目の素数さん
19/01/10 23:01:43.10 HnjZz/GW.net
mean.a <- function(x) "a"
mean(a)
#> [1] "a"
これ本に書いてたんだけど
君たちこれ実行して"a"がでてくる?
自分とこでじっこうしても
> mean(a)
[1] NA
警告メッセージ:
mean.default(a) で: 引数は数値でも論理値でもありません。NA 値を返します
となるんやけど?
いみわからん。
682:１３２人目の素数さん
19/01/10 23:06:00.24 HnjZz/GW.net
あ、一応↓に電子版あるので興味ある人はページ内検索してみてください
URLﾘﾝｸ(adv-r.had.co.nz)
683:１３２人目の素数さん
19/01/10 23:11:23.83 HnjZz/GW.net
あ、上の方から順に入力していったらでたわ
でもJSでOOかじった程度なので全然意味分からん
なにやってんのこれ?
684:１３２人目の素数さん
19/01/10 23:48:17.35 rLaGFww4.net
>>655
オブジェクト指向をちゃんと勉強してくれ。
685:１３２人目の素数さん
19/01/11 19:05:06.45 Q3ISqt43.net
>>654
横からだが、勉強になった。
クラスの自作とか考えたことがなかったけど、今度、俺様クラスを作ってみよう
686:１３２人目の素数さん
19/01/13 19:54:34.41 c90jMv3t.net
>>651
知ってるかもだけどいちいちmappingは書かなくてもいい
scale_shapeも値が1:3ならいらない
ggplot(data, aes(shape=Conditions, col=Conditions))+
geom_line(alpha=0.6)+
geom_point(alpha=0.8)
687:１３２人目の素数さん
19/01/27 04:45:13.06 1PEFhfyS.net
習ってから思ったけど、Fortranとgnuplotでいいよね
688:１３２人目の素数さん
19/01/29 23:19:36.86 6ZawiHpQ.net
普及度、パッケージ数、専門度、文献数の点でそれはない
689:１３２人目の素数さん
19/01/30 01:27:02.50 6Sa5WD/n.net
ヒカキンの年収が10億超
690:え!?明石家さんま・坂上忍も驚愕の総資産とは?? https://logtube.jp/variety/28439 【衝撃】ヒカキンの年収・月収を暴露！広告収入が15億円超え！？ https://nicotubers.com/yutuber/hikakin-nensyu-gessyu/ HIKAKIN(ヒカキン)の年収が14億円！？トップYouTuberになるまでの道のりは？ https://youtuberhyouron.com/hikakinnensyu/ ヒカキンの月収は１億円！読唇術でダウンタウンなうの坂上忍を検証！ https://mitarashi-highland.com/blog/fun/hikakin なぜか観てしまう！！サバイバル系youtuberまとめ http://tokyohitori.hatenablog.com/entry/2016/10/01/102830 あのPewDiePieがついに、初心YouTuber向けに「視聴回数」「チャンネル登録者数」を増やすコツを公開！ http://naototube.com/2017/08/14/for-new-youtubers/ 27歳で年収8億円 女性ユーチューバー「リリー・シン」の生き方 https://headlines.yahoo.co.jp/article?a=20170802-00017174-forbes-bus_all 1年で何十億円も稼ぐ高収入ユーチューバー世界ランキングトップ10 https://gigazine.net/news/20151016-highest-paid-youtuber-2015/ おもちゃのレビューで年間12億円！ 今、話題のYouTuberは6歳の男の子 https://www.businessinsider.jp/post-108355 彼女はいかにして750万人のファンがいるYouTubeスターとなったのか？ https://www.businessinsider.jp/post-242 1億円稼ぐ9歳のYouTuberがすごすぎる……アメリカで話題のEvanTubeHD https://weekly.ascii.jp/elem/000/000/305/305548/ 世界で最も稼ぐユーチューバー、2連覇の首位は年収17億円 https://forbesjapan.com/articles/detail/14474 ヒカルの収入が日収80万、月収2400万、年収3億と判明www https://matomenewsxx.com/hikaru-income-8181.html はじめしゃちょーの年収は6億？2017年は30億突破か？ https://2xmlabs.com/archives/1873
691:１３２人目の素数さん
19/01/30 11:36:05.44 ZoxNNwN2.net
glmer関数にgamma分布を適用したモデルをggplotで回帰曲線を引く方法がわかりません......。summary(model)$dispersionがNULLと返されてしまいます。
692:１３２人目の素数さん
19/01/30 22:41:24.83 eyCKZk6T.net
コードの例示がないからよく分からないけど、もともと戻り値が無いんじゃないの
693:１３２人目の素数さん
19/01/31 17:27:54.13 m67Rusd/.net
?ggplotを見たかい
694:１３２人目の素数さん
19/02/02 11:02:49.05 aYmPih5s.net
dplyr 0.8.0てもう来た？
2月1日とか聞いたような
695:１３２人目の素数さん
19/02/02 21:51:39.00 KjQdurrZ.net
>>658
可読性を保つために省略しないほうがいい
プログラミングの基本
696:１３２人目の素数さん
19/02/02 22:39:28.49 43EjN4Kg.net
逆じゃね。可読性上げるために省略する。
特にmappingなんかは書くまでもない。
697:１３２人目の素数さん
19/02/03 00:29:26.31 Ifbe1wQ3.net
継承を無視してコピペしまくるのが基本なわけない
マウント取りもいきすぎると滑稽だわ
698:１３２人目の素数さん
2019/02/14
699:(木) 03:33:35.52 ID:ilawOvx/.net
700:１３２人目の素数さん
19/02/14 07:38:25.36 0JSgR2gZ.net
データテーブルの1列目のある文字列（例えばstart）を検索してその行含めた上の行を全部削除するスマートな書き方ある？
df［-1:-grep("start",df[,1])］
今こんな感じだけどパイプで繋げない。
slice使えばいいのか、grep以外にいい関数があるのか？
701:１３２人目の素数さん
19/02/14 09:16:23.85 5Ge1SQSk.net
>>670
startが必ず一つしかないとか、restartとかがあったときにそれを検出しても良いならそれでもよいが、普通は
df[-(1:which(df[, 1] == "start")[1]), ]
すると思う。
経験的にはgrepはできるだけ避けた方がいい。
702:１３２人目の素数さん
19/02/14 11:48:21.71 x7agr8k6.net
お邪魔します。
スレ違いが明らかなんですが、ここで聞かせて下さい。
numpyスレ立てていいですか?
703:１３２人目の素数さん
19/02/14 12:08:18.97 9Y2hAxsG.net
もちろん
とんでもなくスレ違い
さっさと立てるね
704:１３２人目の素数さん
19/02/14 12:28:21.65 0JSgR2gZ.net
>>671
なるほど、ありがとうございます。
そしてカッコでくるんでからマイナス付ける書き方もあるのか。勉強になります。
705:１３２人目の素数さん
19/02/14 12:45:38.37 1Nqxk3XQ.net
>>669
可読性について勉強できてよかったな
706:１３２人目の素数さん
19/02/14 17:44:07.85 0JSgR2gZ.net
データフレームのある列が全てNAの時、その列を削除するよい方法ある？
現状col<-apply(2,function(x)｛all(is.na(x))｝
で要らない列定義してからdf[,!col]としてるけどほんとはパイプの中に入れて処理したい。
707:１３２人目の素数さん
19/02/14 23:27:05.62 ATjmovrF.net
>>676
自己解決しました。select_ifでいけそうです。
708:１３２人目の素数さん
19/02/16 12:45:31.76 WMhX0kGV.net
>>672
板違いだと思うので、やめておいた方が良いのでは
709:１３２人目の素数さん
19/02/19 23:32:31.86 TNXJWgVE.net
tidymodelsどう？あまり日本語の資料ないからなあ
710:１３２人目の素数さん
19/02/27 15:36:28.78 3fnUBEYQ.net
RStudioをvi風のキーバインドにすると、
ノーマルモードのときに日本語入力してしまうとバグみたいになるんだが
あれどうにかならんの?
711:１３２人目の素数さん
19/02/28 13:53:56.86 uGjkNcWR.net
どなかた↑よろしく
712:１３２人目の素数さん
19/03/01 22:53:08.03 GbS9kHJG.net
RStudio と日本語入力ってほんと相性わるい
なんとかしてくれんかなぁ
713:１３２人目の素数さん
19/03/02 11:44:51.96 xWBxsoOC.net
windowsのpreview版(1.2.1303)RStudio使ってるけど、
IMEが無効になるバグが直ってる感じがする。
まだ十分使い込んでないからだけかも。
714:１３２人目の素数さん
19/03/03 01:15:23.41 Hp71k0It.net
URLﾘﾝｸ(wired.jp)
シロンボヒトモドキゴキブリニホンザルの自由は偽自由と詐欺広告人を殺す自由ヒトモドキゴキブリシロンボアメ公はニホンザルゴキブリと自殺せよ？
715:１３２人目の素数さん
19/03/04 09:44:29.47 752mbxzc.net
｢いきる｣とか最近ネットで使ってるやつ増えてきたがなんなんあれ?
あほな反抗期の中学生がつかってるイメージしかわかんのだが。
716:１３２人目の素数さん
19/03/04 12:02:26.
717:52 ID:DLycryqB.net
718:１３２人目の素数さん
19/03/04 12:42:02.72 YzqCP1Bw.net
【人類を２つに分ける】 世界教師マⅰトレーヤ到来
ｽﾚﾘﾝｸ(seiji板)
まもなく世界経済破綻、みなさんの財産は紙くずになります。
719:１３２人目の素数さん
19/03/04 19:13:04.98 TnBdKqi2.net
>>686
方言に対して広辞苑を持ち出すのは的外れ。
大阪弁辞典によると「意気る」または「粋る」なので、漢字も外れ
720:１３２人目の素数さん
19/03/04 19:19:22.60 Z/nuwIrB.net
>>688
どっちも「熱る」の当て字やん
721:saga
19/03/05 02:07:35.96 myQXAbVL.net
活きる
722:１３２人目の素数さん
19/03/05 08:37:04.71 agNxkP9Y.net
＞「生きる」と「熱る」
「熱る」って初めて聞いた
これは普通に使われているの？
方言？
生まれてこのかた、東京なんだけど聞いたこと見たことないんだが
単なる私の不勉強かな？
723:１３２人目の素数さん
19/03/05 08:59:52.15 rIYlDwl5.net
URLﾘﾝｸ(dictionary.goo.ne.jp)
いき・る【熱る】
１ あつくなる。ほてる。むしむしする。
２ 激しく怒る。
URLﾘﾝｸ(kotobank.jp)
いき・る【熱る】
① あつくなる。ほてる。むしむしする。
② 息づかいを荒くして怒る。相手と争おうとしていきまく。言いたてる。
③ 調子に乗って勢いこむ。元気づく。
724:１３２人目の素数さん
19/03/06 20:09:12.61 93dcjOku.net
URLﾘﾝｸ(www.youtube.com)
ヒトモドキ奇形なりすましシロンボアメ公はイスラム国に首切られろヒトモドキフィフィラクダイスラム国の売女とともに自殺しろテロゴキブリ民族アメ公
725:１３２人目の素数さん
19/03/07 21:17:38.50 VjN8Z7nf.net
豚肉屋の豚肉自民のキモオタ奴隷豚障害者犯罪者窃盗犯山田太郎と犯罪者キモオタは今すぐ民族レベルで自決自殺しろw
キモ豚山田太郎は自民入党のマイノリティ下僕化しキモオタヒトモドキ障害者の存在価値は0に達したんだよ表現戦士の性犯罪者害虫
このキモオタヒトモドキネトウヨをガス室で抹殺してえよな？w
726:１３２人目の素数さん
19/03/07 22:45:56.15 VjN8Z7nf.net
腐れシロンボテロアメ公ヒトモドキは核で根絶やしになれヒトモドキニホンザルゴキブリの親玉障害者欠陥遺伝子白塵ゴキブリ
727:１３２人目の素数さん
19/03/10 10:47:48.29 FJeOC8+j.net
QAi2Acx7tz
奇形ネトウヨヒトモドキゴミ藪部落の顔気持ち悪いから自殺しろ糞犬hk
728:１３２人目の素数さん
19/03/10 10:49:53.12 QtQWWkXv.net
ww.sankei.com/column/ 180307/clm1803070005-a.html
ヒトモドキ産経便所ゴキブリは下着ドロの常習犯変態暴論ゴミ便所自殺しなさい今すぐ
729:１３２人目の素数さん
19/03/20 12:05:11.55 sWnpvDa6.net
欠損値のないデータで解析したモデルをstepAIC()で処理しようとしたところ、TRUE/FALSEが必要なところが欠損値ですとエラーが出たのですが、どう処理すればよいのでしょうか。MuMInのdredge()でも同様になってしまいました。
730:１３２人目の素数さん
19/03/28 15:23:23.97 ahK9
731:oO7y.net
732:１３２人目の素数さん
19/03/28 20:29:02.09 U0fUTvgM.net
> (p=1-binom.test(22,30)$p.v)
[1] 0.9838752
> binom.test(22,30,conf=p)
Exact binomial test
data: 22 and 30
number of successes = 22, number of trials = 30, p-value = 0.01612
alternative hypothesis: true probability of success is not equal to 0.5
98.38752 percent confidence interval:
0.5000000 0.8994036
sample estimates:
probability of success
0.7333333
733:１３２人目の素数さん
19/03/29 21:03:39.81 g6RZxVSs.net
「統計」は「疑似科学」な
734:１３２人目の素数さん
19/04/06 19:22:34.74 ZoRErNus.net
yahoo apiで距離取得するコード掲載したサイトないですかね
google料金高過ぎて鞍替えなんですが
735:１３２人目の素数さん
19/04/16 20:17:02.49 pCHBksve.net
>>682
使ってるQtが古いんや。1.2でだいぶましになった感じだけど入力途中の文字が残る妙な挙動が時々でる。
736:１３２人目の素数さん
19/04/16 21:58:14.98 aYkLyMdt.net
R書籍漁ったけどまともなのないな
737:１３２人目の素数さん
19/04/17 07:23:57.02 DjFMZd4L.net
>>704
書籍は出版される頃には情報が古くなってるから。自分はネットの情報とヘルプだけで十分だな。
738:１３２人目の素数さん
19/04/18 19:43:25.70 JlSlq7Ln.net
名前が一文字のせいで検索する時に苦労する
739:１３２人目の素数さん
19/04/18 22:13:27.80 vlWwPP8w.net
つ seekR
740:１３２人目の素数さん
19/04/29 14:15:34.67 EFJy96Ez.net
RStudioって32ビット版って存在しないの？？
ダウンロードしようとしても32ビット版見当たらないし、Windows7,10用の64ビット版セットアッププログラムを起動しようすると32ビットへの選択とかはなくて、単にエラーになっちゃうし
誰か教えてください
741:１３２人目の素数さん
19/04/30 14:31:58.91 dkpfnzJo.net
httpstwitter.com/shotkr16
低脳中卒万引きヒトモドキネトウヨ猿ヒトモドキをさしころせ
742:１３２人目の素数さん
19/04/30 17:29:05.11 N9kgT92U.net
>>708
古いのはここにあるけど今後のことをかんがえたらOSを64bit版に変えた方がいいよ。
URLﾘﾝｸ(support.rstudio.com)
743:い
19/04/30 20:32:55.24 DrrYw4j1.net
ここでRをつかった仕事をされてる人いますか？
学生ですか？
仕事あれば教えてほしい
744:１３２人目の素数さん
19/05/01 19:38:51.19 dNg1mxtc.net
>>711
Rを仕事�
745:ﾉしてはいないけど、仕事の一環でR使ってる 信頼性、入手しやすさ、扱いやすさのバランス的に自分の環境ではR一択なんだよな
746:１３２人目の素数さん
19/05/01 23:10:05.28 h9v4bd5W.net
>>712
ちなみにPythonはどう？
747:710
19/05/02 08:41:31.17 fIit4zkT.net
>>713
AI目当てでこれから使おうと思って勉強中
今の職場で求められるのは主に統計用途や作図なんでRは手放せないかな
748:１３２人目の素数さん
19/05/02 13:30:24.26 74c/hxJZ.net
twitter.com/tukuhae
ゴキブリネトウヨヒトモドキ奇形売春婦肉便器なつこババア滅多刺しにして解体しろ
749:１３２人目の素数さん
19/05/02 13:37:47.73 8Iminfbi.net
>>714
なるほど。
自分も機械学習やってるけど、R の方がやりやすいよ。ディープラーニングならPythonかもしれないけど
750:１３２人目の素数さん
19/05/02 18:56:28.17 zw8jEgb5.net
データ分析とか統計解析とかだとRかPythonかだね
機械学習よりは統計寄りだとRは仕事でもスタンダード、特にマニアックな統計解析だとRにしかパッケージがない
まあ使い慣れてるならRで機械学習やっても全然構わないけど
751:１３２人目の素数さん
19/05/03 01:30:52.83 z3krlkgk.net
>>710
サンクス
まあメモリが2Gなので
752:１３２人目の素数さん
19/05/03 01:33:11.23 z3krlkgk.net
>>717
少し前はRなんてだんだん廃れるみたいな雰囲気あったけど今はむしろ盛り返してデファクトな何時になってるからこの先10年は行きそうだなあ
機械学習特需でpythonの躍進はすごかったけど
753:１３２人目の素数さん
19/05/03 18:06:51.34 74TBOTfG.net
2Gって四則演算もままならないだろ
754:１３２人目の素数さん
19/05/11 16:40:31.81 d+R/7VHb.net
Rstudio初心者です。
データファイル名を変えてフォルダには変更名で保存されてるのに、読み込むともとの名前のままなんです。どうすれば変更後の名前になるのでしょう？
755:１３２人目の素数さん
19/05/12 20:48:36.39 R3SzgyK7.net
>>710
サンクス
当面は32ビット版で頑張る
756:１３２人目の素数さん
19/05/17 22:45:05.59 BQfw+Y+R.net
>>721
全く状況がわからん。
主語や目的語、必要な修飾語を省略せずに、分かるように説明しないと誰も助言できない。
757:１３２人目の素数さん
19/05/19 05:50:12.72 i0Ntz+gH.net
いつのまにかどうでもいい国の資格試験の選択科目にまでなったしな
なんなんだろう
758:１３２人目の素数さん
19/05/19 12:44:13.70 5B49QX4v.net
資格を馬鹿にする者は資格に泣く
759:１３２人目の素数さん
19/05/20 08:40:12.19 ys3TVMX0.net
資格にしがみついたらいつのまにか技術の進展から取り残された国がありましたとさw
760:１３２人目の素数さん
19/05/22 10:42:04.74 d0kaOjCb.net
>>721
それはsave, loadよ話か？
761:１３２人目の素数さん
19/07/11 19:34:28.29 A4GNZ5B3.net
Pivot_longer & wider
762:１３２人目の素数さん
19/07/20 11:06:35.86 bSAoQnjE.net
0645
ふうL@Fu_L12345654321
学コン1傑いただきました！
とても嬉しいです！
ｈｔｔｐｓ://pbs.twimg.com/media/D-IuUuqVUAALnAB.jpg
ｈｔｔｐｓ://twitter.com/Fu_L12345654321/status/1144528199654633477
(deleted an unsolicited ad)
763:１３２人目の素数さん
19/08/06 14:27:52.67 Y8L07+xi.net
>
764:臨時地震板 > 【備えあれば】防災用品 非常食スレ121【憂いなし】 > https://mao.5ch.net/test/read.cgi/eq/1562388281/
765:１３２人目の素数さん
19/08/10 11:28:08.57 YonpeJMb.net
ＲってMac版とWindows版で優劣ないの？
今から始めるならどっちがオススメ？
766:１３２人目の素数さん
19/08/10 13:44:59.32 /UVOwF70.net
R自体に優劣はないが文字コード絡みの問題が多いWindows環境はおすゝめしない。
767:１３２人目の素数さん
19/08/10 17:55:47.97 OOGR/e7Y.net
Macはフォントで苦労するけどな
768:１３２人目の素数さん
19/08/10 21:06:01.89 7gPmJrDz.net
じゃ、Linuxだ。
769:１３２人目の素数さん
19/08/10 21:23:29.54 YonpeJMb.net
フォントで苦労ってどういう種類の？
文字化け？
770:１３２人目の素数さん
19/08/12 04:58:06.03 TcXjNI8s.net
>>731
R本体だとグラフィックディバイスに違いはあるけど、実質的に差はない。
でもRは本体だけで使うことは稀で、多種多様なGUIがあるから、
そこでOSの違いによる差が生じる。
771:１３２人目の素数さん
19/12/09 08:29:17 g2fJs3Gj.net
面白い問題スレにあったのでシミュレーションしてみた。
# サイコロ
# 正６面体のサイコロがある．４面は青色、２面は赤色である．
# このサイコロを合計20回振るとき、最も起こりそうな順番はどれか？
# １．赤 青 赤 赤 赤
# ２．青 赤 青 赤 赤 赤
# ３．青 赤 赤 赤 赤 赤
sim <- function(){
a=sample(0:1,20, replace=TRUE, prob=c(4,2))
b=as.character(a)
c=paste(b,collapse="")
s1=paste(c(1,0,1,1,1),collapse="")
s2=paste(c(0,1,0,1,1,1),collapse="")
s3=paste(c(0,1,1,1,1,1),collapse="")
res=c(grepl(s1,c),grepl(s2,c), grepl(s3,c))
return(res)
}
k=1e6
re=replicate(k,sim())
mean(re[1,])
mean(re[2,])
mean(re[3,])
結果は、直感とおり、１が再頻
> mean(re[1,])
[1] 0.124672
> mean(re[2,])
[1] 0.080873
> mean(re[3,])
[1] 0.040564
>
grep使わない方法ってあるかな？
772:１３２人目の素数さん
19/12/10 18:31:14.56 KL6XgkNo.net
これって確率分布的には何になるんですっけ？
773:１３２人目の素数さん
19/12/11 10:39:16.38 xruQypao.net
歪んだコインやんか
774:１３２人目の素数さん
19/12/12 16:40:23 8G7On9nx.net
>>737
ワイの直感的解法
# １．赤 青 赤 赤 赤
# ２．青 赤 青 赤 赤 赤
# ３．青 赤 赤 赤 赤 赤
だが、以下でも確率同じ。何となく
# １．赤 赤 赤 赤 青 ★
# ２．赤 赤 赤 赤 青 青
# ３．赤 赤 赤 赤 青 赤
★は赤でも青でもどっちでもOK
P(# １) >P(# ２) >P(# ３) ∵青出やすい
R言語等でシミュレーションされ、
自身の確率直感が正しいのを
確認できるとは、素晴らしい。
775:>>738
19/12/12 18:31:03 8G7On9nx.net
上記の件、若干の訂正とする
# １．赤 青 赤 赤 赤
# １'．赤 赤 赤 赤 青 とすると、
# １と# １'は、直感で同じ確率と
思ってたが間違えのようだ。
当方のシミュレーションで、
# １は、0.1248
# １'は、0.1271
となった。微妙だけど、多分だ。
やっぱり確率計算をコンピュータで
モンテカルロシミュレーションのは
素晴らしい。
776:１３２人目の素数さん
19/12/12 19:21:51.88 8G7On9nx.net
>>737
【grepは未使用の糞真面目な方法】
# １．赤 青 赤 赤 赤
についての確率、ほぼ厳密解を得た
# １は、0.124774… だと思う
計算は、モンテカルロ法でない方法
でプログラム、計算した。
で、grepは使用してない。
ちなみに計算誤
777:差は、ほぼ皆無なハズ ソースコードイメージ p01 = (1/3)^4*(2/3) p02 = p01 p03 = p01 * (1 - p01) p04 = p01 * (1 - p01 - p02) p05 = p01 * (1 - p01 - p02 - p03) … p16 = p01 * (1 - p01 - p02 - p03 … - p14) とし、 p01~p16 の合計を算出したところ、 0.124774… となった
778:１３２人目の素数さん
19/12/13 07:46:08 lh0xGphL.net
>>740
レスありがとうございます。
私の直感
# １．赤 青 赤 赤 赤
# ２．青 赤 青 赤 赤 赤
# ３．青 赤 赤 赤 赤 赤
# １．を６個に書き換えて #２．と並べると
# １．★ 赤 青 赤 赤 赤
# ２．青 赤 青 赤 赤 赤
★は赤でも青でもどっちでもOKだから#１．の方が起こりやすい
# ２と# ３を比べると
# ２．青 赤 青 赤 赤 赤
# ３．青 赤 赤 赤 赤 赤
３個めでは青の方がでやすいので
ら#２．の方が起こりやすい
よって、P(# １) >P(# ２) >P(# ３)
779:１３２人目の素数さん
19/12/14 03:40:27.74 L1XaepqW.net
分からない問題スレから、
>>
1回3.6%で激レアが出るガチャを10回回した確率って
36%なのでしょうか？
それとも0.964*0.964*0.964(略 0.964を10回電卓にかけた数なのでしょうか？
教えてください。
<<
百万回シミュレーション
p=3.6/100
N=10
sim <- function() any(rbinom(N,1,p)==1)
mean(replicate(1e6,sim()))
> mean(replicate(1e6,sim()))
[1] 0.306628
780:１３２人目の素数さん
19/12/14 03:51:08 L1XaepqW.net
シミュレーション その２
gacha=c(rep(1,36),rep(0,1000-36))
sim2 <-function() any(sample(gacha,10,replace=TRUE)==1)
mean(replicate(1e6,sim2()))
> mean(replicate(1e6,sim2()))
[1] 0.306904
何故か0.36にならない、どうしてだろ？
781:１３２人目の素数さん
19/12/14 04:10:22.87 L1XaepqW.net
シミュレーションその３ （処理速度の関係で１０万回の平均）
gacha=c(rep(1,36),rep(0,1000-36))
sim3 <- function() any(replicate(10,sample(gacha,1))==1)
mean(replicate(1e5,sim3()))
> mean(replicate(1e5,sim3()))
[1] 0.30691
これも0.307弱だな。
何が悪いんだろ？ 俺の頭かな？？
782:１３２人目の素数さん
19/12/14 05:45:07.83 qYB5MEbs.net
3.6%のガチャを10回回して全部外れる確率は
(1 - 0.036) ^ 10 ≒ 0.6930592
したがって、
3.6%のガチャを10回まわして1回以上当たる確率は
1 - 0.6930592 = 0.3069408
783:１３２人目の素数さん
19/12/14 08:18:13 L1XaepqW.net
>>747
ありがとう、
やっぱ、
これだったｗ
何が悪いんだろ？ 俺の頭かな？？
784:１３２人目の素数さん
19/12/14 08:19:41 L1XaepqW.net
>>748
シミュレーションにバグはなかったｗ
785:１３２人目の素数さん
19/12/14 12:29:26.79 cW38CJV+.net
小学生でも解けるぞ。がんばれ。
786:１３２人目の素数さん
19/12/14 16:29:58.85 QSWNkRNd.net
確率統計なんぞ無意味つまりインチキ
事実、結果こそすべて
787:１３２人目の素数さん
19/12/14 19:21:15.54 1VaGaO0p.net
その確率計算の激ナイーブな解法を示す
激レアは、レアだから1個とみなす。
故に、母集団の個数は
1÷0.036 = 27.777… きっと28個だ
手順1) 28枚のカードがある
手順2) 重複しない1~28の番号を振る
手順3) 28枚のカードをシャッフル
手順4) 1~10枚目のどれか1となる確率
と絶対同じハズ、だから、
P(1枚目で当) =27P27 ÷ 28P28 = 1/28
P(2枚目で当) も同様に1/28
P(3枚目で当) も同様に1/28
…
P(10枚目で当) も同様に1/28
で此等10個の事象は背反事象だから、
P = 10/28 = 0.357
∵有効数字3桁と勝手にしちゃう
ところてガチャってゲーム何か
よくわかんないけど、計算しちゃった
788:１３２人目の素数さん
19/12/14 20:44:12.02 1VaGaO0p.net
追記というか突然ですが、
そのガチャ3.6%の件、
超幾何分布なのか。二項分布なのか。
確率の小さいとか、母集団が小さい とかだと無視できないと思われる。
確率統計はギャンブル派生数学ぢゃ。
生半可な知識ではカモにされる。
現代の若者たちは、数学特に
確率統計はじめとするギャンブル
の能力が特段に欠けており、
R言語等のプログラミング教育で
ギャンブルゲームを学習すべきだ。
奇麗事の学問だけの今日の数学ぢゃ
カモにされるだけ。
健全な娯楽として賭博系確率統計学
をC R Java Pyson Javascript BASIC
の何れかを学校で学習すべきだ。
789:１３２人目の素数さん
19/12/16 09:18:31.40 t3hxiPBr.net
私のガチャのイメージは、
かつて、昭和の駄菓子屋によくある
ガチャガチャすなわちカプセルトイ
そのような健全的なギャンブルが
超幾何分布とか二項分布の理解に
役立つのだ。
限られたお小遣の10円玉数枚で、
如何にレアアイテムをGetするかを
子供らは、思考するからだ。
「残り物には福がある」は
確率統計的には正しいのか
子供同士で文学的に議論したものだ。
さて、今のガチャは恐らくは、
デジタルの媒体のスマホゲームだ。
二項分布でよいだろう。
R言語には、
二項分布の密度関数、それの累積関数
はモチロン、それに従う乱数生成を
提供してるようだ。
ゲームの仕組みが複雑化する今、
R言語等の乱数生成プログラミングで
これからのデジタル化ギャンブル社会
で、お金より大切な激レアをドンドン
無限にゲットできる人材の増加を
期待できる可能性を秘めている
790:１３２人目の素数さん
19/12/17 06:37:18.85 gnX7/8eN.net
カプセルトイ（ガチャ）１台に1000個カプセルが入っていて36個がアタリ(レアアイテム）とする。
同じカプセルトイが10台ある。カプセル取り出し後は補充されない。
アタリを１個でも手に入れる確率は
１台から10個取り出す場合(G10)と
１台から１個を１０台で取り出す場合(G01)
ではどちらが高いか？
そのシミュレーション
rm(list=ls())
N=1000 ; K=36 ;n=10 # アタリ3.6％
g=rep(c(1,0),c(K,N-K))
G10 <- function() any(sample(g,n,replace=FALSE)>0) # 非復元（超幾何分布）
G01 <- function() any(sample(g,n,replace=TRUE )>0) # 復元(二項分布)
mean(replicate(1e6,G10()))
mean(replicate(1e6,G01()))
791:１３２人目の素数さん
19/12/17 08:10:38.97 gnX7/8eN.net
シミュレーションと理論値
> mean(replicate(1e6,G10())) ; 1-choose(N-K,n)/choose(N,n)
[1] 0.307745
[1] 0.3081121
> mean(replicate(1e6,G01())) ; 1-(1-K/N)^n
[1] 0.307295
[1] 0.3069408
１台から１０個取り出した方がいいみたい。
792:１３２人目の素数さん
19/12/17 14:14:51.68 JPqkotiS.net
>>756 なるほど、
どの台も全部で、1000カプセルで、
どの台も当りが、36カプセルだと
G10(1台で10個)取り出す方が
僅かですが、確率良さそうですね。
シミュレーションも理論値も
同様な結論のようであり、
G10(1台で10個)取り出す方が僅かに
有利が分かり何か楽しかったです。
仮に、もしどの台も
250カプセル中当たり9カプセルなら、さらにG10戦略がG01戦略より有利な
感触を掴めました。
793:１３２人目の素数さん
19/12/17 15:44:13.73 gnX7/8eN.net
>>757
1000個だとシミュレーションの差が微妙で再現性が不安だったけど
250個にすると、差が明らかにつきますね。
> N=250 ; K=9 ;n=10 # アタリ3.6％
> g=rep(c(1,0),c(K,N-K))
>
> G10 <- function() any(sample(g,n,replace=FALSE)>0) # 非復元（超幾何分布）
> G01 <- function() any(sample(g,n,replace=TRUE )>0) # 復元(二項分布)
>
> mean(replicate(1e6,G10())) ; 1-choose(N-K,n)/choose(N,n)
[1] 0.31082
[1] 0.3117069
> mean(replicate(1e6,G01())) ; 1-(1-K/N)^n
[1] 0.306782
[1] 0.3069408
794:１３２人目の素数さん
19/12/18 14:18:29.16 GQllmBub.net
rhyper で超幾何分布の乱数を発生させることができたんだな。
これと rbinomを使うと sample関数を使わなくても同じことができた。
N=250 ; K=9 ;n=10 # アタリ3.6％
mean(replicate(1e6,rhyper(1,K,N-K,n)>0)) ; 1-choose(N-K,n)/choose(N,n) # 非復元（超幾何分布）
mean(replicate(1e6,any(rbinom(n,1,K/N)>0))) ; 1-(1-K/N)^n # 復元(二項分布)
結果
> mean(replicate(1e6,rhyper(1,K,N-K,n)>0)) ; 1-choose(N-K,n)/choose(N,n)
[1] 0.311581
[1] 0.3117069
> mean(replicate(1e6,any(rbinom(n,1,K/N)>0))) ; 1-(1-K/N)^n
[1] 0.307213
[1] 0.3069408
1台から順次１０個取り出すのをイメージすれば残り物には福があるとも言えなくもないな。
>
795:１３２人目の素数さん
19/12/18 21:18:37 GQllmBub.net
"フランスの数学者パスカル（1623～1662)が1654年にフェルマーにあてた手紙が、現在の確率
論の始まりだと言われている。当時の有名な賭博師メレがパスカルに以下のような問題を持ち
込み、その問題についてがその手紙のやりとりの中で論じられているそうである。
甲乙二人がおのおの32ピストル（当時のお金の単位）の金を賭けて勝負したとする。
そしてどちらかが先に３点を得たものを勝ちとし、勝った方がかけ金の総額64ピストルをもら
えるとする。ところが甲が２点、乙が１点を得たとき、勝負が中止になってしまった。
このとき、二人のかけ金の総額64ピストルを甲と乙にどのように分配すればよいだろうか。
ただし二人の力は互角で、勝つ確率はそれぞれ1/2ずつだとする。"
# 先にw(=3)点を得たものを勝ち、甲がA(=2)点、乙がB(=1)点を得たとき、勝負が中止
# 甲の確率を求める
gambling <- function(A=2,B=1,w=3,k=1e5){ # k : how many simulate
sim <- function(){
while(A < w & B < w){
g = rbinom(1,1,p=0.5)
if(g==1){
A=A+1
}else{
B=B+1
}
}
A > B
}
mean(replicate(k,sim())) # Pr[A wins]
}
> gambling(2,1,3)*64 #64ピストル配分
[1] 48.13056
> gambling(2,0,4) # 日本シリーズ
[1] 0.80945
796:１３２人目の素数さん
19/12/19 09:59:55.45 V+OT4hGF.net
>>760の計算結果は、興味深い結果だ
で、
自己流の直感でラフな計算をしてみた。
甲1点、乙1点 ⇒1:1按分 ∵自明 ★
甲3点、乙1点 ⇒1:0按分 ∵自明 ☆
安直に★と☆の平均を考えると
甲2点、乙1点 ⇒1:0.5按分 ∵安直
で、総額で64ピストルだから、
甲に、64*1/1.5 = 42.67 ピストル
乙に、64*0.5/1.5 = 21.33 ピストル
を配分すれば、思ってしまった。
ラフな計算だと、少々不公平なことが
起こるようだ。
R言語等のシミュレーションは、
超天才のワィの直感より、
さらに厳密な確率計算ができそうだ。
797:１３２人目の素数さん
19/12/19 16:21:08.12 kVFzAP55.net
>>760
最後はAが得点して勝負が決定するのを忘れずにAが勝者になるまでのBの点数で場合分け
NS <- function(w,A,B,p){ # 先にw点得点した方が勝者、A,B：現在の点数 ,p:甲の勝率
ans=0
for(k in 0:(w-B-1)){ # k: Aが勝者になるまでのBの点数
ans=ans+choose(w-A-1+k,w-A-1)*(1-p)^k*p^(w-A)
}
return(ans)
}
> NS(3,2,1,0.5)*64 #64ピストルの分配
[1] 48
> NS(4,1,0,0.5) #日本シリーズで先勝したほうが優勝する確率
[1] 0.65625
NSは日本シリーズの頭文字から命名、No Skinではないｗｗ
798:１３２人目の素数さん
19/12/19 18:00:34.70 V+OT4hGF.net
賭博等の確率計算をイメージすると
ワィの脳内は
799:活性化される だから、 ドンドン、妄想全開となるぅぅぅ。 という訳で、妄想内容を記載する 17世紀の地球に行き、胴元になって 賭博ビジネスでガッポリ儲けるゼェ。 17世紀の地球で、胴元になって、 コイントス賭博事業をするゼェ。 で、システムは、2回コイントスし、 表の回数が0回⇒参加者が32万円GET 表の回数が1回⇒胴元に、32万円PAY 表の回数が2回⇒参加者が32万円GET 参加料は、1万円 まんまと、参加者が沢山 1000人いた ある参加者は、以下の様に怪説をした 組み合わせは、3通りだ。 1通り目 表の回数が0回 2通り目 表の回数が1回 3通り目 表の回数が2回 だから、参加者の勝つ確率は2/3だ で、モチロン、胴元ワィの商売は成功 およそ、1000万円儲かっちゃた。
800:１３２人目の素数さん
19/12/20 02:07:38.50 yiLw1Jz8.net
0745
しろ@huwa_cororon 11月27日
苦節6ヶ月、初満点&一等賞です！
URLﾘﾝｸ(twitter.com)
(deleted an unsolicited ad)
801:１３２人目の素数さん
19/12/21 10:49:06.51 +RV2yvS7.net
日本シリーズで賭けをする。
日本シリーズは先に4勝したチームが優勝。
勝率はそれまでの（引き分けを除いた）通算勝率に従うとする。
シリーズ開始前の通算成績はA:2勝、B:4勝であった。
今シリーズでAが先勝（第一試合に勝利）した。
この時点でA,Bのどちらに賭ける方が有利か？
802:１３２人目の素数さん
19/12/21 22:03:14.56 vLicX+v2.net
分からない問題スレで３進法の小数とかが話題になっていた
ｽﾚﾘﾝｸ(math板:845番)
>>
845 名前：１３２人目の素数さん[] 投稿日：2019/12/21(土) 05:52:50.74 ID:wum1jR1j
...
10進数で表された5/9を3進数になおせという問題があって、答は0.12だったのですが、
...
<<
こんなのを作ってみた。
# 小数点付きの数numをN進法で表示する
rm(list=ls())
dec2basen <- function(num, N, kmin = 5){ # kmin：最小小数点後桁
int=floor(num)
r=int%%N
q=int%/%N
while(q > 0){
r=append(q%%N,r)
q=q%/%N
} # rに整数のN進法表示を格納
x=num-floor(num)
k=max(nchar(x)-2,kmin) # 同長もしくはkminの長さの小数表示
a=numeric(k)
for(i in 1:k){
y=x*N
a[i]=floor(y)
x=y-a[i] # r . a[1] a[2] a[3] ... a[k]
}
if(N<=10){ # Nが10以下は数値として表示
cat(r,paste(".",paste(a,collapse = ''),sep=''),'\n',sep='')
}
else{ # Nが1１以上は整数部分と小数部分を数列で表示
print(list(int=r,decimal=a))
}
invisible()
}
> dec2basen(5/9,3)
0.120000000000000
> dec2basen(3.14159265359,7)
3.06636514320
> dec2basen(3.14159265359,8)
3.11037552421
> dec2basen(3.14159265359,9)
3.12418812407
> dec2basen(3.14159265359,14)
$int
[1] 3
$decimal
[1] 1 13 10 7 5 12 13 10 8 1 4
803:１３２人目の素数さん
19/12/22 08:18:39.30 Ez3wWboE.net
>>766
バグの原因を追求したら、Rの仕様のせいみたいだな。
> (1.2-1)*5==1
[1] FALSE
> (1.2-1)*5>1
[1] FALSE
> (1.2-1)*5<1
[1] TRUE
これがfloor関数が誤値を返してくる原因だった。
> floor(1)
[1] 1
> floor((1.2-1)*5)
[1] 0
804:１３２人目の素数さん
19/12/22 08:52:46.43 UYCjcQMO.net
>>767
数値演算についてもっと勉強しな。
たいていの言語でそうなる。
805:１３２人目の素数さん
19/12/22 13:06:36.02 Ez3wWboE.net
すると、pythonに移植しても同じ結果ということか。
IPython 6.5.0 -- An enhanced Interactive Python.
C:\bin\Anaconda3\lib\site-packages\ipykernel\parentpoller.py:116: UserWarning: Parent poll failed. If the frontend dies,
the kernel may be left running. Please let us know
about your system (bitness, Python, etc.) at
ipython-dev@scipy.org
ipython-dev@scipy.org""")
(1.2-1)*5==1
Out[1]: False
(1.2-1)*5>1
Out[2]: False
(1.2-1)*5<1
Out[3]: True
(1.2-1)*5
Out[4]: 0.9999999999999998
806:１３２人目の素数さん
19/12/22 17:15:43.69 Ez3wWboE.net
round を使って回避することにした
> 0.728*5-3.64
[1] -4.440892e-16
> round(0.728*5-3.64,5)
[1] 0
んで、デバッグ版
# 小数点付きの数numをN進法で表示する(62進法まで対応 0-9,a-z,A-Z)
rm(list=ls())
dec62 <- function(num, N, kmin = 5){ # kmin：最小小数点後桁
int=floor(num)
r=int%%N
q=int%/%N
while(q > 0){
r=append(q%%N,r)
q=q%/%N
} # rに整数のN進法表示数列を格納
k=max(nchar(num)-nchar(floor(num))-1,kmin) # 同長もしくはkminの長さの小数表示
a=numeric(k)
x=round(num-floor(num),k) # e.g. 7.28-floor(7.28)-0.28 != 0に対応
for(i in 1:k){
y=round(x*N,k) # e.g. 0.728*5-3.64 !=0 に対応
a[i]=floor(y)
x=y-a[i] # r . a[1] a[2] a[3] ... a[k]
}
b=list(integer=r,decimal=a,num=sum(c(int,a)*(1/N)^(0:k)))
fig=c(0:9,letters,LETTERS)[1:N]
if(N<=62){ # Nが62以下は数値として表示
cat(paste(fig[b$integer+1],sep='',collapse=''),
'.',paste(fig[b$decimal+1],sep='',collapse=''),sep='')
cat('\n')
}
else{ # Nが63以上は整数部分と小数部分を数列で表示
print(b[1:2])
}
invisible(b) # b$num:検証用
}
807:１３２人目の素数さん
19/12/22 17:25:01.53 Ez3wWboE.net
>>769
Haskellでもやってみた
Prelude> (1.2-1)*5
0.9999999999999998
Prelude> 0.72*5-3.6
-4.440892098500626e-16
808:１３２人目の素数さん
19/12/25 07:12:30.44 oEKznZ6+.net
>>768
レスありがとうございました。
アドバイスに従ってちょっと勉強してみました。
他の言語に移植しても無駄とわかって助かりました。
こうなる理由が理解できました。
(1+1/10-1)*10==1
[1] FALSE
> (1-1+1/10)*10==1
[1] TRUE
809:１３２人目の素数さん
20/01/06 12:56:45.01 e9wyGMBv.net
からnまでの順列を列挙するスクリプト。
これをwhile使って高速化するにはどうすればいいだろう？
俺の頭では思いつかない。
perm <- function (n) {
v=1:n
sub <- function(n, v) {
if (n == 1)
matrix(v, 1, 1)
else {
x = NULL
for (i in 1:n) x =rbind(x, cbind(v[i], sub(n -1,v[-i])))
x
}
}
sub(n, v)
}
> perm(4)
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 1 2 4 3
[3,] 1 3 2 4
[4,] 1 3 4 2
[5,] 1 4 2 3
[6,] 1 4 3 2
[7,] 2 1 3 4
[8,] 2 1 4 3
[9,] 2 3 1 4
[10,] 2 3 4 1
[11,] 2 4 1 3
[12,] 2 4 3 1
[13,] 3 1 2 4
[14,] 3 1 4 2
[15,] 3 2 1 4
[16,] 3 2 4 1
[17,] 3 4 1 2
[18,] 3 4 2 1
[19,] 4 1 2 3
[20,] 4 1 3 2
[21,] 4 2 1 3
[22,] 4 2 3 1
[23,] 4 3 1 2
[24,] 4 3 2 1
>
810:１３２人目の素数さん
20/01/06 21:34:47.52 yym51Tg7.net
>>773
自己解決 e1071というパッケージに再帰呼び出しなしのpermutationがありました。
library(e1071)
permu <- function (n) {
z <- matrix(1)
for (i in 2:n) {
x <- cbind(z, i)
a <- c(1:i, 1:(i - 1))
z <- matrix(0, ncol = ncol(x), nrow = i * nrow(x))
z[1:nrow(x), ] <- x
for (j in 2:i - 1) {
z[j * nrow(x) + 1:nrow(x), ] <- x[, a[1:i + j]]
}
}
return(z)
}
811:１３２人目の素数さん
20/01/09 14:28:06.78 rRnLkyt7.net
可変引数の扱い
Cで
#include <stdio.h>
int main(int argc, char* argv[])
{
while(argc--)
printf("%s\n", *argv++);
return 0;
}
をRでやるには ... を使って、c(...)もしくはlist(...)で引き出せばいいんだな。
検索してもなかなかわからなかった。
main <- function(...){
argv=c(...)
for(i in 1:length(argv)) cat(argv[i],'\n')
}
812:１３２人目の素数さん
20/01/15 18:08:54 l229ykjv.net
>>548
library(Rmpfr)
mpfr(factorial(50),5000)
one = mpfr(1, 5000)
factorial(one*50)
> mpfr(factorial(50),5000)
1 'mpfr' number of precision 5000 bits
[1] 30414093201713018969967457666435945132957882063457991132016803840
> one = mpfr(1, 5000)
> factorial(one*50)
1 'mpfr' number of precision 5000 bits
[1] 30414093201713378043612608166064768844377641568960512000000000000
何故か、oneを高精度に設定しておくと、正しい計算をしてくれた。
813:１３２人目の素数さん
20/01/26 16:33:21 G7gVG9Ku.net
RのroundのIEEE仕様
> round(2.5)
[1] 2
> round(1.5)
[1] 2
は丸め誤差が減るという説明だったので実験してみた。
# 四捨五入 vs ＩＥＥＥ with FUN(=mean) for x(=c(-0.9,-0.8,...,0.8,0.9))
comp <- function(n=10,FUN=mean,x=sample((-9:9)/10,n,rep=T),print=T){
X=FUN(x)
dif = abs(X-FUN(round(x))) - abs(X-FUN(f45(x))) # round後の実行と四捨五入後の実行の差
if(dif!=0 & print) cat(dif<0,' : ',sort(x),'\n') # 差があれば表示 round後が小さければTRUE
return(dif) # 差を返す dif<0:round優位 dif>0:四捨五入優位
}
comp()
k=1e5
# mean
re=replicate(k,comp(print=F))
c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
# squared sum
comp(FUN=function(x) sum(x^2))
re=replicate(k,comp(FUN=function(x)sqrt(sum(x^2)),print=F))
c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
> # 平均 mean
> re=replicate(k,comp(print=F))
> c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
ＩＥＥＥ 四捨五入 引き分け
0.25235 0.16199 0.58566
> # 平方和 squared sum
> re=replicate(k,comp(FUN=function(x)sqrt(sum(x^2)),print=F))
> c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
ＩＥＥＥ 四捨五入 引き分け
0.40519 0.01413 0.58068
平均だと小差だが、統計計算の内部処理で使う平方和だと大差がついたので、体感できた。
814:１３２人目の素数さん
20/02/08 17:42:30.00 +JttXwIS.net
"隔離中のクルーズ船では船内の換気が共通らしいから運悪く13日後に発症した奴がいるとその近くの部屋のやつは
プラスで14日隔離しないといけない。それが今の船内の状況という。
両隣のどちらかが感染したら14日延長、どの部屋も1日で感染する確率pは1%
部屋の配置は長方形でどの部屋にも両隣があるとする。
発症するか、隔離期間が終われば下船できる。
全員定員１の個室として客と乗務員を合わせた人数nは3000人。
クルーズ船から全員下船できる日数の期待値は？"
815:１３２人目の素数さん
20/02/08 17:42:46.53 +JttXwIS.net
rm(list=ls())
p=0.001 # 1日の感染確率
n=3000 # 収容者人数
d14=14 # 隔離日数
Q0=rep(d14,n) # 必要隔離日数の配列の初期値
#
ODA <- function(X){ # One Day After 必要隔離日数の配列 x から１日後の配列を返す
Q=X #
816:Xを作業用Qに代入 S=which(Q!=0) # 未感染者susceptibleのindex s=length(S) # 未感染者人数 m=rbinom(1,s,p) # その日の感染数 if(m==0){ # 感染者0なら未感染者の Q[S]=Q[S]-1 # 必要隔離日数を１日減らして返す invisible(Q) }else{ I.idx=sample(S,m) # 未感染者Sからm人が感染、感染者のindex (Infected.idx) id2ext <- function(id){ # 感染者idからの両隣のidを返す plus =ifelse(id==n,1,id+1) # 最大番号なら1を返す,それ以外は1増やす minus=ifelse(id==1,n,id-1) # １なら最大番号を返す,それ以外は１減らす unique(c(plus,minus)) # 重複を除く } E.idx=as.vector(sapply(I.idx,id2ext)) # 感染者の両隣のidを列挙し E.idx=unique(E.idx) # 重複を除く Extended quarantine.idx Q[E.idx]=d14 # 両隣は必要隔離日数をd14にリセット Q[I.idx]=0 # 感染者は必要隔離日数=0 (リセット者に重複していたら０で上書き) IE.idx=unique(c(I.idx,E.idx)) # 感染もしくは感染者の隣のindex minus1ifnotzero <- function(x) ifelse(x==0,0,x-1) # ゼロでなければ１を減じて返す Q[-IE.idx]=minus1ifnotzero(Q[-IE.idx]) # 感染もなく感染者の隣りでもない人の隔離日数を１日減らす invisible(Q) } } sim <- function(){ AZ=FALSE # All Zero 隔離日数がすべて０か？ i=0 # カウンター Q=Q0 # 隔離日数配列初期値 while(!AZ){ # 隔離日数がすべて0になるまで Q=ODA(Q) # One Day After関数を実行 AZ=all(Q==0) # All zero？ i=i+1 # カウンター } return(i) # 何日かかったかを返す } k=1e4 # k回のシミュレーションの RE=replicate(k,sim()) mean(RE) # 平均値＝期待値
817:１３２人目の素数さん
20/02/08 18:25:56.61 Xaabj2Ml.net
感染確率0.1％やん
低すぎぃ
818:１３２人目の素数さん
20/02/08 19:50:08.59 +JttXwIS.net
>>780
一日の感染確率。
3000人いて1日3人くらい報告されたから0.1%/日くらいかな
819:１３２人目の素数さん
20/02/08 19:51:16.44 +JttXwIS.net
実行したら全員が下船できるまで約1ヶ月になった。
コードにまだバグがあるかも。
820:１３２人目の素数さん
20/02/09 21:31:53 IDQ3+Iu7.net
Rのcsvの文字化け治らねえ
csvそのものをutf-8にしたら今度は不正なマルチバイトがあるってエラー
どうすればいいか教えて
821:１３２人目の素数さん
20/02/09 22:07:23 vcjg+5qP.net
>>783
どんなデータ？
ヘッダーだけじゃなくてセル内も日本語？
822:１３２人目の素数さん
20/02/09 23:10:12 KodPd2ez.net
>>783
エスパーじゃないから、まず確認だけど
・Rを使ってる環境(OS)とRのバージョンは？
・CSVのコード変換に使ってるソフトは？
・CSVを読み込んでる関数と指定している引数は？
823:１３２人目の素数さん
20/03/06 19:53:26 WwE0ymvi.net
nCov2019 for studying COVID-19 coronavirus outbreak URLﾘﾝｸ(guangchuangyu.github.io)
後で遊んでみよう
824:１３２人目の素数さん
20/03/07 09:46:33.43 3M7vmA2q.net
>>786
エラーがでた。R 3.6.3にアップデートしたけど
> library('remotes')
> remotes::install_github("GuangchuangYu/nCov2019", dependencies = TRUE)
Downloading GitHub repo GuangchuangYu/nCov2019@master
Skipping 1 packages not available: BiocStyle
Installing 5 packages: downloader, cowplot, maps, magick, BiocStyle
Error: Failed to install 'nCov2019' from GitHub:
(converted from warning) package ‘BiocStyle’