【R言語】統計解析フリーソフトＲ 第6章【GNU R】 [無断転載禁止]©2ch.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
# >>597
患者来院数４０人程度では定常状態に達しないということみたいだな。
１０万人での待ち時間の分布（理論値５０に近い）
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 来院患者数と待ち時間シミュレーションの結果をグラフにしてみた。
個人医院で待ち行列の長さが一定になるほどの受診があるとは考えがたいからこっちの方が現実に即しているんじゃなかろうか
と思うが、解析解はどんなふうになるのか想像もつかん。
https://i.imgur.com/tCOoxF7.png 初歩的な質問ですいません
truehistで軸を対数軸にして表示する方法って分かりますか？ truehistだと以前の質問でお礼すらしなかった人か >>602
答えたの俺だがソース改造できたのかなぁ。 histで縦軸を対数表示させるならこんな感じかな。
with(hist(c(rnorm(1e6),rnorm(1e5,5,0.5))),plot(mids,log10(counts),type='h',col=4,lwd=10))
truehistだとソース改造すればいい。 >>604
plot(log='y')でもいいな。 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]) すみませんがアドバイスお願いします。
heatmap.2でDendrogramつきヒートマップを描きたいのですが、カラムの並びを任意に変えたいです。
Dendrogramの細かい並びは変えないように、大きなクラスタの並びを変えたいです。
例えば、（１，２，３）（４，５，６）（７，８，９）とあるのを
（４，５，６）（７，８，９）（１，２，３）とならび変えるのが目的です。
このとき、4,5,6と7，8，9は近いクラスタを形成します。なので、樹形図は崩れないように書き換えられると思っています。
https://www.biostars.org/p/237067/
上記サイトをみると、as.dendrogramのアウトプットをreorderで並び替えてColvに入れるようですが、
うまくいきませんでした。並びは全く変わっていません。
どなたか教えていただけますか。情報に漏れがありましたらご指摘ください。
環境は以下の通りです。
R 3.4.0
Rstudio 1.0.143
Win10になります。 頓珍漢な答かもしれない。
並べ替えだけなら
x=rbind(c(1,2,3),c(4,5,6),c(7,8,9))
x[c(2,3,1),] listなら
x=list(c(1,2,3),c(4,5,6),c(7,8,9))
x[c(2,3,1)] 初歩的な話なんだが、一様分布の分散は無限大かと思っていたら区間[a,b]で(a-b)^2/12とのこと。
Wolframで計算したら確かにそうなった。
https://www.wolframalpha.com/input/?i=integral(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 NULLのときってどうしてこういう仕様なんだろ？
プログラムしていたら、これに気づかないのがバグの原因だったｗ
> any(NULL)
[1] FALSE
> all(NULL)
[1] TRUE = と <-で微妙に動作が違うな。
> switch (3,
+ x =1,
+ x =2,
+ x =3
+ )
[1] 3
> switch (2,
+ x <- 1,
+ x <- 2,
+ x <- 3
+ ) >>612
そのふたつは意味が異なる
それぞれちゃんと命令どおりの挙動 どちらが見やすいかという問題かな。
> rm(x)
> x=switch(1,x =1)
> x
[1] 1
> switch(1,x<-1)
> x
[1] 1 >>614
その二つは異なるアルゴリズムでたまたま結果が同じになっているだけ
=と<-の違いはRやるなら理解しておいたほうが良い 0^x =0
x^0=1
0^0=1とした方が辻褄が合うことが多いけど
Rのこの仕様には何のメリットがあるんだろ？
> any(NULL)
[1] FALSE
> all(NULL)
[1] TRUE
> >>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,
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) pushってありますか?
例えばベクトルに値を追加すると
先頭が消えて後尾に新しい値をついかして要素数を一定に保つような
一定要素数以下の平均を求めたいのでそういうのが簡単に実現できる方法あればなおよいです ベクトルをrev()して先頭を[1:50]とかでとりだしてman()すればいいとわかりました
ほかにもっと簡単な方法アレばお湿気てください 組み合わせると
f=function(x,y,n=50) mean(tail(append(x,y),n)) 時系列データを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=?????)
として縦線を追加したいのですが v=as.POSIXct("2019-01-06 01:00")
みたいな感じにすれば解決しました >>278
かなり遅れすだけど
formals(cor) >>611
>>617
俺の予想
判定するときにNULLは強制的に論理値に変換される
変換されると logical(0) になる
logical(0) は空のベクトル
からのベクトルが渡されて中身をチェックしていく
アルゴリズムとして早く処理するためには
any()はTRUE探しにいって、一個でもTRUEがみつかればその時点でTRUEをかえす
all()はFALSEを探しに行って、一個でもFALSEがみつかればその時点でFALSEをかえす
もちろん最後までみない
空のベクトルが渡されたのでTRUEもFALSEもみつからない、となると
any()ではFALSEとなり
all()ではTRUEとなる
これで辻褄はあう
ちなみに
any(c())
all(c())
でも同じ結果が出る 1行のテキストデータを最終的に数値とテキストの混在するN行M列のデータフレームにしたいのですが、なかなかうまく出来ません。
1行データの構造の設計からデータフレームの変換までどうすればシンプルに実現できるか助言ください。
たとえば
1 a
2 b
というようなデータを
1-a,2-b
というような一行のテキストデータから初めてテーブル構造にするというような形です
x <- "1-a,2-b"
y <- str_split(x, ",")
と試しにやってみたのですがyが行単位のベクトルになるだけでここからどうデータフレームにすればよいかわかりません >>630
何をしたいのか、まったく理解できない。最低限、N, Mが何なのか説明したらどう？ >>630
文字列 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]) 実行結果
> 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
> >>631
NMは任意の数字です
したいことは
データを一行に記録してそれを
テーブル構造にすることです。
これだけです。
一旦ファイルに保存してread.csvなどにすればよいのでしょうが
いちおう直接テキストデータをコピペしてからということにしたいのです。
なぜこんなことをするかと言うと
TamperMonkeyという拡張機能でJavaScriptで使ってウェブ上のデータを収集しているのですが、
localStrageというクッキーの拡張版のような機能をつかってデータをクライアントに保存するときに基本的にテキストデータベースでの保存になので
一行のデータとして後ろにどんどんデータを追加していくのが一番単純な処理に成るからです。
そのlocalStorageに保存されたデータはコピペで取り出すしかないので
それをRで処理する時に一行のテキストデータから始めないといけないのです。 >>634
, を行のデリミタ、- を列のデリミタとして、一行の文字列をデータフレームにするというのであれば、
s に文字列が入っているとして、
r <- unlist(strsplit(s, ","))
d <- lapply(r, function(x) unlist(strsplit(x, "-")))
as.data.frame(do.call(rbind, d))
列数が行によって異なるときにどうなるかは知らん。 >>629
内部動作の考証ありがとうございます。
未だにこれは理解できません
> logical(NULL)
Error in logical(NULL) : invalid 'length' argument
> logical(0)
logical(0)
> logical(1)
[1] FALSE >>630
> 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
こんな感じか？ >>636
自分の解釈では、
NULL は「無」なのでエラー（引数はベクトルの要素数を必ず与えなければいけない）
0のとき、0個という指定なので、0個の要素をもつベクトル
1のとき、1個という指定なので、1個の要素を持つベクトル（規定値はFALSE）
ちなみに、
> logical(2)
[1] FALSE FALSE
2のとき、2個という指定なので、２個の要素を持つベクトル（規定値はFALSE） 皆さんありがとうございます。
>>637
ファイル入れなくてもできるんですね 最古の元号って大化なん?
一巡回って大化2でいいよ
あとは干支みたいに回せばいい >>641は誤爆ですw
レス付くと思わなかった。。。 >>639
⇣の結果をすべて即答できるようになれば合格だと思います…
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 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) すみません自己解決しました。
このコードではなくもっと後のthemeで凡例のスタイルを変えるときに余計なことをしていました。 mean.a <- function(x) "a"
mean(a)
#> [1] "a"
これ本に書いてたんだけど
君たちこれ実行して"a"がでてくる?
自分とこでじっこうしても
> mean(a)
[1] NA
警告メッセージ:
mean.default(a) で: 引数は数値でも論理値でもありません。NA 値を返します
となるんやけど?
いみわからん。 あ、上の方から順に入力していったらでたわ
でもJSでOOかじった程度なので全然意味分からん
なにやってんのこれ? >>653
オブジェクト指向をちゃんと勉強してくれ。 >>652
横からだが、勉強になった。
クラスの自作とか考えたことがなかったけど、今度、俺様クラスを作ってみよう >>649
知ってるかもだけどいちいちmappingは書かなくてもいい
scale_shapeも値が1:3ならいらない
ggplot(data, aes(shape=Conditions, col=Conditions))+
geom_line(alpha=0.6)+
geom_point(alpha=0.8) 習ってから思ったけど、Fortranとgnuplotでいいよね 普及度、パッケージ数、専門度、文献数の点でそれはない glmer関数にgamma分布を適用したモデルをggplotで回帰曲線を引く方法がわかりません......。summary(model)$dispersionがNULLと返されてしまいます。 コードの例示がないからよく分からないけど、もともと戻り値が無いんじゃないの dplyr 0.8.0てもう来た？
2月1日とか聞いたような >>656
可読性を保つために省略しないほうがいい
プログラミングの基本 逆じゃね。可読性上げるために省略する。
特にmappingなんかは書くまでもない。 継承を無視してコピペしまくるのが基本なわけない
マウント取りもいきすぎると滑稽だわ >>666
お前のほうが滑稽
マウントってるのお前
イキってるだけで中身もないし
滑稽すぎる データテーブルの1列目のある文字列（例えばstart）を検索してその行含めた上の行を全部削除するスマートな書き方ある？
df［-1:-grep("start",df[,1])］
今こんな感じだけどパイプで繋げない。
slice使えばいいのか、grep以外にいい関数があるのか？ >>668
startが必ず一つしかないとか、restartとかがあったときにそれを検出しても良いならそれでもよいが、普通は
df[-(1:which(df[, 1] == "start")[1]), ]
すると思う。
経験的にはgrepはできるだけ避けた方がいい。 お邪魔します。
スレ違いが明らかなんですが、ここで聞かせて下さい。
numpyスレ立てていいですか? >>669
なるほど、ありがとうございます。
そしてカッコでくるんでからマイナス付ける書き方もあるのか。勉強になります。 データフレームのある列が全てNAの時、その列を削除するよい方法ある？
現状col<-apply(2,function(x)｛all(is.na(x))｝
で要らない列定義してからdf[,!col]としてるけどほんとはパイプの中に入れて処理したい。 >>674
自己解決しました。select_ifでいけそうです。 >>670
板違いだと思うので、やめておいた方が良いのでは tidymodelsどう？あまり日本語の資料ないからなあ RStudioをvi風のキーバインドにすると、
ノーマルモードのときに日本語入力してしまうとバグみたいになるんだが
あれどうにかならんの? RStudio と日本語入力ってほんと相性わるい
なんとかしてくれんかなぁ windowsのpreview版(1.2.1303)RStudio使ってるけど、
IMEが無効になるバグが直ってる感じがする。
まだ十分使い込んでないからだけかも。 https://wired.jp/2019/01/18/get-wired-kevin-kelly-5-videos/
シロンボヒトモドキゴキブリニホンザルの自由は偽自由と詐欺広告人を殺す自由ヒトモドキゴキブリシロンボアメ公はニホンザルゴキブリと自殺せよ？ ｢いきる｣とか最近ネットで使ってるやつ増えてきたがなんなんあれ?
あほな反抗期の中学生がつかってるイメージしかわかんのだが。 「生きる」と「熱る」のどっちだ？
いずれにしろ広辞苑に載ってるから調べればいい >>684
方言に対して広辞苑を持ち出すのは的外れ。
大阪弁辞典によると「意気る」または「粋る」なので、漢字も外れ ＞「生きる」と「熱る」
「熱る」って初めて聞いた
これは普通に使われているの？
方言？
生まれてこのかた、東京なんだけど聞いたこと見たことないんだが
単なる私の不勉強かな？ https://dictionary.goo.ne.jp/jn/10539/meaning/m0u/
いき・る【熱る】
１ あつくなる。ほてる。むしむしする。
２ 激しく怒る。
https://kotobank.jp/word/%E7%86%B1%E3%83%BB%E7%86%85-2005439
いき・る【熱る】
① あつくなる。ほてる。むしむしする。
② 息づかいを荒くして怒る。相手と争おうとしていきまく。言いたてる。
③ 調子に乗って勢いこむ。元気づく。 https://www.youtube.com/channel/UCKyNitDuOcc7Gi1cjQlNCwA
ヒトモドキ奇形なりすましシロンボアメ公はイスラム国に首切られろヒトモドキフィフィラクダイスラム国の売女とともに自殺しろテロゴキブリ民族アメ公 豚肉屋の豚肉自民のキモオタ奴隷豚障害者犯罪者窃盗犯山田太郎と犯罪者キモオタは今すぐ民族レベルで自決自殺しろw
キモ豚山田太郎は自民入党のマイノリティ下僕化しキモオタヒトモドキ障害者の存在価値は0に達したんだよ表現戦士の性犯罪者害虫
このキモオタヒトモドキネトウヨをガス室で抹殺してえよな？w 腐れシロンボテロアメ公ヒトモドキは核で根絶やしになれヒトモドキニホンザルゴキブリの親玉障害者欠陥遺伝子白塵ゴキブリ QAi2Acx7tz
奇形ネトウヨヒトモドキゴミ藪部落の顔気持ち悪いから自殺しろ糞犬hk ww.sankei.com/column/ 180307/clm1803070005-a.html
ヒトモドキ産経便所ゴキブリは下着ドロの常習犯変態暴論ゴミ便所自殺しなさい今すぐ 欠損値のないデータで解析したモデルをstepAIC()で処理しようとしたところ、TRUE/FALSEが必要なところが欠損値ですとエラーが出たのですが、どう処理すればよいのでしょうか。MuMInのdredge()でも同様になってしまいました。 https://www3.nhk.or.jp/news/html/20190327/k10011863181000.html
インフルエンザの新しい治療薬「ゾフルーザ」を投与されたＡ香港型のインフルエンザ患者30人を調べたところ、70％余りに当たる22人から、この薬が効きにくい耐性ウイルスが検出されたことが国立感染症研究所の調査で分かりました。
調査件数は多くないものの、専門家は現在のような使用を続けると、耐性ウイルスが広がるおそれがあるとして使用基準を見直すべきだと指摘しています。
耐性化率が50%以上である確率は
pbeta(0.5,1+22,1+8,lower=F)
[1] 0.9946631でいいかな？ ■ このスレッドは過去ログ倉庫に格納されています
