18/08/04 20:37:20.39 kYh4t8ZN.net
>>470
Stanとかに対応してんの？
488:１３２人目の素数さん
18/08/04 22:03:56.61 kTWj7ufe.net
>>470
化石みたいなバージョンやな
489:１３２人目の素数さん
18/08/05 00:56:32.45 biDT5wEJ.net
>>470
いまSPSSを使う理由って何なの？煽りじゃなくマジで
490:１３２人目の素数さん
18/08/05 14:27:54.32 dSJLM/da.net
maple使えよ
491:１３２人目の素数さん
18/08/05 19:10:34.43 mULz5b3r.net
>>473
周囲のspssユーザを見ると「指導教員がspssしか使えないので、
自分もspssしか使えない、
新しい（学術系）ソフトは誰かに丁寧に根気よく指導してみらえないと無理」、
これの再生産。
492:１３２人目の素数さん
18/08/05 19:11:30.48 zaJ/WJYb.net
>>470
研究室に予算がないから
493:１３２人目の素数さん
18/08/05 19:12:42.57 mULz5b3r.net
>>476
予算がなかったら、Rを使え
494:１３２人目の素数さん
18/08/05 19:13:59.36 mULz5b3r.net
>>470と>>476を間違えたorz
495:１３２人目の素数さん
18/08/05 19:33:14.43 mTZcaJ/n.net
予算がないならケーキを食べればいいのに
496:１３２人目の素数さん
18/08/08 13:11:20.17 7/3l/PxD.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
497:１３２人目の素数さん
18/08/09 17:47:22.65 l1M8GWWv.net
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
498:１３２人目の素数さん
18/08/09 18:31:24.68 l1M8GWWv.net
簡略化できた
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
>
499:１３２人目の素数さん
18/08/10 18:41:42.48 Hlm8Oe3x.net
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
500:１３２人目の素数さん
18/08/10 18:55:39.89 vdkgfABT.net
>>483
ifelse() はベクトル演算できる関数。
だから、戻り値の要素数は、test部分の要素数に一致する。
testが1つなので、x[2:1]の1つ目の要素が返される。
501:483
18/08/10 18:58:51.59 vdkgfABT.net
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
502:１３２人目の素数さん
18/08/10 23:06:17.77 Hlm8Oe3x.net
>>484
ありがとうございました。
バグではなくてそういう仕様だったのですね。
こういう使い方ができるということがわかって勉強になりました。
> x=2:1
> ifelse(c(x[1]!=x[2],x[1]==x[2]),1:2,3:4)
[1] 1 4
503:１３２人目の素数さん
18/08/10 23:17:14.10 VlOUWHrC.net
教わったので早速、
ifelse() はベクトル演算できると、再利用されるの動作確認。
> ifelse(c(TRUE,TRUE,FALSE,FALSE,TRUE,FALSE,FALSE),1:2,11:15)
[1] 1 2 13 14 1 11 12
504:１３２人目の素数さん
18/08/10 23:45:24.37 VlOUWHrC.net
２種類の文字の変換も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,] 真 真 偽
505:１３２人目の素数さん
18/08/11 00:28:09.90 OesSEWnz.net
行列も次元付きのベクトルだからもっと簡略化できた。
> (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
>
506:１３２人目の素数さん
18/08/11 18:04:43.46 05wxJPem.net
RやSPSSって名前の試験まであるのね
507:１３２人目の素数さん
18/08/16 23:26:21.81 3mKQ9SP5+
>>469
適合度検定の場合は、不適の基準は期待度数5未満ではなく、2.5未満では?
つまり一様性の検定なら、全ての期待度数が2.5以上なら良いのでは?
508:１３２人目の素数さん
18/08/22 14:02:15.21 rznk0lAS.net
規約分数にするパッケージが探せなかったので自分でスクリプトを組んでみた。
エラー処置は省略。
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
509:１３２人目の素数さん
18/08/23 01:12:39.77 0Wz4IoKN.net
rvest?でログイン必要なとこをスクレイピングする時って
login_page <- html_session("URLﾘﾝｸ(xxxxx"))
login_form <- html_form(login_page)[[1]] %>%
set_values(AAA="xxxx", BBB="xxxx")
というので行けるはずなんですが
最終行のAAAとBBBってそれぞれhtmlタグのname=""のとこからとってくるんですよね?
これにここがname="user[email]"みたいに[]という記号はいってたらどうすればいいでしょうか?
510:１３２人目の素数さん
18/08/25 17:48:59.21 MdsDkupV.net
>>493
よろしくおねがいいしまーーーーーーーす
511:１３２人目の素数さん
18/08/26 22:01:57.06 2bj/HrEX.net
>>494
質問の意味がいまいち分からないから、誰も助言できないのでは？
> x <- 'user[email]'
> x
[1] "user[email]"
rvestパッケージを使ったことがないけど、
記号が入っていたらなぜ問題があるのかよく分からない。
512:１３２人目の素数さん
18/08/26 22:53:43.10 WnqFRbi9.net
というかどこのページのスクレイピングしたいか晒して
仮でもいいから
513:１３２人目の素数さん
18/08/27 06:16:23.20 x8E4FG0O.net
>>495 ちがいます AAAがどれかよくみてみてください ””はありません
515:１３２人目の素数さん
18/08/27 06:19:43.00 x8E4FG0O.net
ちょっと直接はることはできないのでタグだけはります。
<input class="form-control" autofocus="autofocus" placeholder="Email address" required="required" type="text" value="" name="user[email]" id="user_email">>>496
516:１３２人目の素数さん
18/08/27 13:25:39.67 SdOxff6m.net
>>497
もっと分かるように説明しないと、追試できる情報も提供しないし。
もしかして変数名などに[]が入っている場合 にどうしたらよいかってこと？
> `user[email]` <- 1:5
> `user[email]`
[1] 1 2 3 4 5
517:１３２人目の素数さん
18/08/27 15:43:37.46 x8E4FG0O.net
>>499
多分行けました
超感謝!!
518:１３２人目の素数さん
18/08/27 20:20:32.13 5a+trmgU.net
ある問題のシミュレーションしようと思って
問題文の記号のまま
q <- function(x)
....
とやって
q(100)と入力すると
Rが終了することに気づいた。
519:１３２人目の素数さん
18/08/27 22:56:54.68 a+8R0Tu1.net
base::q()
を実行すると爆弾出るね
520:１３２人目の素数さん
18/08/29 22:45:10.53 BcwFyR33.net
ちょっとした疑問です。
空ベクトルの検出って長さ０以外に検出方法ってあるでしょうか？
> 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
521:１３２人目の素数さん
18/08/31 22:46:15.16 IWQvY6FL.net
４点の座標を入力するとそれらを結ぶ四面体の体積を求めるスクリプトを書いてみた。
高さはパッケージ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
多分、正常に動作していると思う。
522:１３２人目の素数さん
18/09/01 00:25:33.05 52Ub52jp.net
>>504
行列式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
523:１３２人目の素数さん
18/09/01 01:41:47.56 qG52f2Ee.net
ベクトルの三重積を教わったので、パッケージ 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)
}
４行で済んだ。
>>505
ありがとうございました。
パッケージに頼らずに計算できたのですね。
524:１３２人目の素数さん
18/09/03 18:55:17.77 S47YTHgP.net
データを解析する前にさらっと特徴を見たい時、皆さんはどんなコマンドを使っていますか？
私が思いつくのは
summary
boxplot
hist
pairs
です。こんなのも良いよってのがあったら教えてくださいm(_ _)m
※ライブラリの使用有無は問いません
525:１３２人目の素数さん
18/09/03 19:38:36.45 fowZfPON.net
>>507
BESTパッケージのplotPost
histに加えて95%CIを表示してくれる。
526:１３２人目の素数さん
18/09/03 20:05:14.66 OGj8hrn2.net
>>507
出力をデータフレーム型にしたいときはskimrパッケージ
527:１３２人目の素数さん
18/09/03 20:07:16.33 OGj8hrn2.net
あと、まだ試してないけsummarytoolsパッケージも面白そう
528:１３２人目の素数さん
18/09/09 22:35:22.23 9XY+z1xx.net
>>503
良い方法が見つからなかったので !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
529:１３２人目の素数さん
18/09/09 23:21:54.06 1udV8jVZ.net
>>511
なぜ if (length(x) == 0) と書かないのか？
可読性って知ってるか？
530:１３２人目の素数さん
18/09/10 16:54:23.64 89eIPezd.net
>>512
文脈読めば
!length(x)で空白ベクトル
の流れ
531:１３２人目の素数さん
18/09/10 18:22:05.74 gLTPxqtw.net
>>513
length は数値を返す関数なのに、なぜわざわざ論理演算をするんだ？
プログラムは、ここの議論を知らない人が読む可能性のほうが高いんだぞ。
532:１３２人目の素数さん
18/09/10 19:12:52.95 ZYY4OYkH.net
>>514
正規分布が一様分布より大きい期待値の算出に
mean(rnorm(100)>runif(100))
って書かない？
俺は
sum(ifelse(rnorm(100)>runif(100),1,0))/100
って書いたりしないけど。
533:１３２人目の素数さん
18/09/10 19:43:25.60 SKQ8XoAj.net
>>514
え、常套手段やん？！
534:１３２人目の素数さん
18/09/10 19:50:35.36 ZYY4OYkH.net
>>516
論理値を数値に置き換えて計算しているから
数値を論理値にしても別に違和感がない。
可読性は慣れの問題。
535:１３２人目の素数さん
18/09/10 19:54:17.32 SKQ8XoAj.net
>>517
アンカー間違えとる？
だから、そんなの常套手段やん？
536:１３２人目の素数さん
18/09/10 21:30:23.37 gLTPxqtw.net
>>515
それとこれとは話が別だ。
logical は TRUE か FALSE の二値しかとらないから、それから数値への変換は自明。
多種の値をとる数値に論理演算をするのはいただけない。
537:１３２人目の素数さん
18/09/10 21:36:20.92 gLTPxqtw.net
>>517
is.empty.vector のような関数があるならそれでよいが、length を使っているのだから、それを勝手にベクトルが空かどうかの判断に使っているのはあなたであって、他人はそのようには考えない、といっているのだ。
自分だけしかこーどを見ないならべつに構わないが、こういうところに晒すのはよくない。
ましてや常套手段などと言って他人に教え込むのはやめてもらいたい。
538:１３２人目の素数さん
18/09/10 21:51:19.40 5QS5/GHY.net
>>520
常套手段というのはmeanの話じゃね？
539:１３２人目の素数さん
18/09/10 21:55:50.32 5QS5/GHY.net
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
}
540:１３２人目の素数さん
18/09/10 22:00:57.10 5QS5/GHY.net
>511は
配列を逆順に並べる再帰呼出しのコード。
Rには
541:revという関数がある。 そのコードみてみ！ lengthを真偽判定に使ってますがな。 > base:::rev.default function (x) if (length(x)) x[length(x):1L] else x
542:１３２人目の素数さん
18/09/10 22:18:05.27 DJR6mPp+.net
>>520
>523みた？
if(length(x)!=0)になってないよ。
543:１３２人目の素数さん
18/09/11 00:19:39.45 WPUC3B84.net
>>520
>こういうところ
誰も鵜呑みにしない便所の落書き。。。
544:１３２人目の素数さん
18/09/11 08:48:42.34 QUqp/jpE.net
!を引数が数値のときは0か否かを返す関数と理解すれば
if(!0) print(!1)の結果もサクッとわかる。
545:１３２人目の素数さん
18/09/11 09:08:32.20 QUqp/jpE.net
数値を非零かどうか論理値に変換して処理するかは
関数によるな。
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)
>
546:１３２人目の素数さん
18/09/11 09:12:56.04 QUqp/jpE.net
!ってベクトル対応しているな
> which(!(-1:1))
[1] 2
> !(-2:2)
[1] FALSE FALSE TRUE FALSE FALSE
547:１３２人目の素数さん
18/09/15 08:05:49.52 Vl7XZ52q.net
!をnotではなくis.falseとかis.zero変換すれば( ・∀・)ｲｲ!!
548:１３２人目の素数さん
18/09/15 16:42:00.50 h0gUCZ3r.net
>>529
is.zero だと？
だったら最初から == 0 てよいじゃないか。
549:１３２人目の素数さん
18/09/15 16:45:38.59 h0gUCZ3r.net
>>527
だいたい、if や while は関数じゃないし。
while (1) などと書くのはCに毒されているんじゃねーの？
Rなら while (TRUE) のほうが良いし、repeat というのもある。
550:１３２人目の素数さん
18/09/15 16:45:41.31 Vl7XZ52q.net
>>530
スクリプトを読むときに脳内変換する癖をつけておけば
可読性が悪いと思わずにすむ。
551:１３２人目の素数さん
18/09/15 16:46:13.42 Vl7XZ52q.net
>>520
>523みた？
if(length(x)!=0)になってないよ。
552:１３２人目の素数さん
18/09/15 17:19:35.45 h0gUCZ3r.net
>>533
Rのライブラリ？の書き方なんてあてにならないよ。
そもそも、rev も 509 の reverse も演算の回数を必要最低限にするなら、条件は length(x)) >= 2 または > 1 だ。
まあ if (length(x)) のほうが速かったのかも知れないが、常にそうなるとは限らないのだから、早すぎる最適化は諸悪の根源だ。
553:１３２人目の素数さん
18/09/15 19:30:34.55 Vl7XZ52q.net
可読性の次は速度かよ。
自分の流儀と違う { の使い方でも文句いいそう。
554:１３２人目の素数さん
18/09/15 19:44:27.43 VibLIqgl.net
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
555:１３２人目の素数さん
18/09/15 19:57:45.72 VibLIqgl.net
!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
結果は再現された。
556:１３２人目の素数さん
18/09/15 20:06:37.73 VibLIqgl.net
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
御指摘の通り、!の負け
557:１３２人目の素数さん
18/10/04 17:02:48.04 247Ted9r.net
>>445
＞遅レスだが、Ubuntu 18.04LTSならデフォルトのIBus-mozcで日本語入力ができるそうだ。
＞詳しくはUbuntu We
558:ekly Recipeの第527回を読んでくだされ Ubuntu 18.04.1 LTSにアップグレードしたけどできない いままでスクレイピングにつかってrvestのログインコマンドもエラー出るようになったし最悪だ。。。。
559:１３２人目の素数さん
18/10/04 20:04:49.68 W6fWyA5T.net
アップグレードってことは文字入力がFcitxのままになってない？Fcitxならパッチ要だよ。
Voyager 18.04LTS（xubuntu 18.04LTS）ではパッチなしでRStudioに日本語入力できた。
560:１３２人目の素数さん
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弱だな。
何が悪いんだろ？ 俺の頭かな？？