18/06/04 13:05:22.57 3QOAvZa7.net
::なら名前空間カブールよけよけの術だよね
421:１３２人目の素数さん
18/06/04 14:24:00.43 d7l4LqX9.net
>>408
とりあえずこれ
URLﾘﾝｸ(github.com)
422:１３２人目の素数さん
18/06/05 12:44:14.88 LfwL/Kok.net
>>409
:::と3個なんだよなぁ
423:１３２人目の素数さん
18/06/05 13:06:03.27 CHJNgrVq.net
>>411
グダグダ言っていないで、さっさとヘルプを参照すればいいじゃないか？
?':::'
424:１３２人目の素数さん
18/06/05 19:45:28.11 LfwL/Kok.net
>>412
ありがとう。
stats:::t.test.default
でｔ検定のソースが見られた
425:１３２人目の素数さん
18/06/11 00:14:18.19 K1zGYlRd.net
関数を値で上書きしてしまった場合に元に戻す方法を教えて頂きたく。
例)
rm <- 1
としてしまった場合にrm関数を元に戻したいです。
よろしくお願いします。
426:１３２人目の素数さん
18/06/11 07:21:35.69 +ayC4RGv.net
>>414
直後ならRstudioでctrl+zで元に戻す。
再起動してたらこの方法は無理だろな
427:１３２人目の素数さん
18/06/11 08:24:34.24 2ZHCUqeW.net
>>414
自分で作ったのではない関数の rm だったら rm(rm) で数値ベクトルオブジェクトが消えて、関数の rm が見えるようになる。
428:１３２人目の素数さん
18/06/11 15:01:37.34 mlFyU0v4.net
>>414
> rm <- 1
> ls()
[1] "rm"
> rm(rm)
> ls()
character(0)
> rm
function (..., list = character(), pos = -1, envir = as.environment(pos),
[以下略]
深く考えなくてもrm(rm) で元に戻るけど。
rmが数値ではなくて自作関数だったら、少々ややこしいけど、
base::rm()で大丈夫だろう
429:１３２人目の素数さん
18/06/11 17:04:35.54 Z+okZT62.net
>>414
既に説明されてるけど上書きじゃなくて別オブジェクト扱いになるからrmで消すか明示的にパッケージを::演算子で関数を指定してやればおけ
自作関数だと上書きされちゃうけど
430:１３２人目の素数さん
18/06/11 19:33:14.14 rfEFkwvW.net
rmが自作関数だと戻せないのでは？
> rm = function(x) asin(sqrt(x))
> rm(1)
[1] 1.570796
> rm
function(x) asin(sqrt(x))
> rm = 1
> rm
[1] 1
> rm(rm)
> rm
function (..., list = character(), pos = -1, envir = as.environment(pos),
inherits = FALSE)
431:１３２人目の素数さん
18/06/11 19:51:41.60 Z+okZT62.net
同じ名前空間にあるなら上書きされちゃうよ
.Gloalenv（だったかな？）に自作関数があり、そこに上書きして変数にしてるから消したら消えるだけで戻せない
432:１３２人目の素数さん
18/06/11 19:55:10.50 Z+okZT62.net
ちょっと古いけど、この辺りを読むとなぜそうなるかが分かると思う
（Rの）環境問題について その１。
URLﾘﾝｸ(qiita.com)
433:１３２人目の素数さん
18/06/11 20:04:41.37 mlFyU0v4.net
>>419
戻せます。
> rm <- function(x) asin(sqrt(x))
> rm
function(x) asin(sqrt(x))
> base::rm(rm)
> rm
function (..., list = character(), pos = -1, envir = as.environment(pos),
[以下略]
434:１３２人目の素数さん
18/06/11 20:06:37.76 mlFyU0v4.net
>>420
>.Gloalenv
case sensitiveなRでそれはないのでは？
> help(".GlobalEnv")
435:１３２人目の素数さん
18/06/11 20:19:14.99 Z+okZT62.net
フォローありがとう
436:413
18/06/11 21:00:49.78 K1zGYlRd.net
>>415-423
みなさまご返信ありがとうございました。
例示したものがbase::rmだったので自作関数のことまでは想像だにしませんでした…(がとても良い勉強になりました。)
自環境でも上書きした関数を元に戻すことが出来ることを確認しました。
437:１３２人目の素数さん
18/06/12 07:42:40.88 9RfEtlLW.net
>>422
戻したいのが自作関数
rm = function(x) asin(sqrt(x))
だったら上書きされて無理じゃないの？
ゴミ箱から消去したファイルを復活させるソフトもあるから
PCに詳しければ可能かもしれないけれど
Rの機能としては復活させるのは無理と思う。
438:１３２人目の素数さん
18/06/12 11:00:45.94 cD/mQB+6.net
自作のrm関数はbaseパッケージとは別に定義されるから
base::rm(rm)で自作の方のrmは消えてbase::rmが「復活」します。
439:１３２人目の素数さん
18/06/12 12:07:52.53 cD/mQB+6.net
ああすまん
自作の方を上書きした場合は、そっちの復活は無理ですね。
440:１３２人目の素数さん
18/06/12 14:36:07.63 VF7OcVBc.net
>>426
ちょっと意味がわからない。
自作の関数なら定義を再実行するだけなのに、
消える消えないって何が問題なんだ？
441:１３２人目の素数さん
18/06/13 14:39:27.58 ZXkz+qgi.net
>>425
関数だけでなく既定値でも同じ。
> pi
[1] 3.141593
> pi=2
> pi
[1] 2
> rm(pi)
> pi
[1] 3.14159
442:１３２人目の素数さん
18/06/14 09:27:55.98 mxBGyFKT.net
共同ツール 1
URLﾘﾝｸ(seleck.cc)
URLﾘﾝｸ(trello.com)
ボードのメニュー → Power-Upsから拡張可能 Slack DropBoxなど
Trello Chrome拡張機能 elegant
URLﾘﾝｸ(www.kikakulabo.com)
trelloのオープンソースあり
共同ツール 2
URLﾘﾝｸ(www.google.com)
共同ツール 3
URLﾘﾝｸ(slack.com)
URLﾘﾝｸ(www.dropbox.com)
URLﾘﾝｸ(bitbucket.org)
URLﾘﾝｸ(ja.atlassian.com)
URLﾘﾝｸ(www.sketchapp.com)
URLﾘﾝｸ(photoshopvip.net)
URLﾘﾝｸ(goodpatch.com)
Trello Chrome拡張機能プラグイン集
URLﾘﾝｸ(chrome.google.com)
Slackプラグイン集
URLﾘﾝｸ(slack.com)
Sketchプラグイン集
URLﾘﾝｸ(sketchapp.com)
URLﾘﾝｸ(supernova.studio)
443:１３２人目の素数さん
18/06/16 19:29:24.89 10mr4q2e.net
はっとデシャンｸﾞ
444:１３２人目の素数さん
18/06/17 19:06:33.58 OYjqtCQI.net
>>307
知恵袋に漸化式があったので計算スクリプトを書いてみた。
表がでる確率ｐのコインをN回投げてK回以上表が続く確率。
# N: flips
# K: least sequential head
# p: probability of head
seqNp <- function(N=100,K=5,p=0.5){
if(N==K) return(p^K)
q=1-p
a=numeric(N) # a(n)=P0(n)/p^n , P0(n)=a(n)*p^n
for(i in 1:K) a[i]=q/p^i # P0(i)=q for any i
for(i in K:(N-1)){ # recursive formula
a[i+1]=0 # a[i+1]=q/p*(a[i]+a[i-1]+a[i-2]+...+a[i-(K-1)])
for(j in 0:(K-1)){
a[i+1]=(a[i+1]+a[i-j])
}
a[i+1]=q/p*a[i+1]
}
P0=numeric(N) # P0[n] : probability of ending with 0 head when flipped n times
for(i in 1:N) P0[i]=a[i]*p^i # P0(n)=a(n)*p^n
MP=matrix(rep(NA,N*K),ncol=K)
colnames
445:(MP)=paste0('P',0:(K-1)) MP[,'P0']=P0 MP[1,'P1']=p for(i in (K-2):K) MP[1,i]=0 # MP[k,n] = Pk[n] : probability of ending with k head when flipped n times for(k in 2:K){ # Pk(n+1)=p*P(k-1)(n) for(i in 1:(N-1)) MP[i+1,k]=p*MP[i,k-1] } ret=1-apply(MP,1,sum) ret[N] }
446:１３２人目の素数さん
18/06/18 09:16:00.35 JdABdagZ.net
>>307
この例で、答えとして何が返ってきて欲しいのか？
2で合ってる？
447:１３２人目の素数さん
18/06/18 21:23:58.67 fHaLluDH.net
>>434
1が5回以上出現しているからTRUE（数字なら１）が返り値。
448:１３２人目の素数さん
18/06/19 10:20:04.23 fcA67BpN.net
>>435
run.check <- function (x, n=5) {
# x の中に 1 が n 回以上連続していれば TRUE を返す.
chg <- c(TRUE, diff(x) != 0) # 変化があった場所
chgidx <- c(which(chg), length(x)+1) # 変化があった場所の添え字
run.length <- diff(chgidx) # 0や1の連続している個数
true.length <- run.length[x[chg] == 1] # 1の連続している個数
any(true.length >= n) # 連続している個数が n 以上のrunがあるか?
}
449:１３２人目の素数さん
18/06/19 21:57:17.85 QuKrVrE8.net
>>436
whichとdiffを用いての解法のお返事ありがとうございます。
実行時間を>303のスクリプトと比べてみました。
run.check <- function (x, n=5) {
# x の中に 1 が n 回以上連続していれば TRUE を返す.
chg <- c(TRUE, diff(x) != 0) # 変化があった場所
chgidx <- c(which(chg), length(x)+1) # 変化があった場所の添え字
run.length <- diff(chgidx) # 0や1の連続している個数
true.length <- run.length[x[chg] == 1] # 1の連続している個数
any(true.length >= n) # 連続している個数が n 以上のrunがあるか?
}
# N(=100)回コインをなげてｎ(=5回)以上続けて表がでる確率。
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(10^5,run.check(rbinom(100,1,0.5)))))
user system elapsed
17.56 0.04 17.68
> system.time(mean(replicate(10^5,seqn())))
user system elapsed
9.74 0.07 9.78
後者の方が速いのは1がn個連続したら、TRUEを返して終了するので以後のチェックはしないためだろうと思います。
diff や anyの使い方が勉強になりました。
時間を割いてスクリプトを作成していただいてありがとうございます。<(_ _)>
450:１３２人目の素数さん
18/06/21 10:12:41.96 Ze2kEGUX.net
Rstudioで日本語の入力ができないのですが、解決する方法知ってる方いませんか?
具体的には、IME自体をオンにすることができません。
日本語の表示やコピペは問題なくできます。
環境
Linux - Utubntu
RStudio Version 1.0.143
451:１３２人目の素数さん
18/06/21 10:26:53.84 Fds/4hNR.net
RStudioのQtが古いのでパッチ当てないと日本語入力できないのが現状ですわ
パッチはUbuntu 16.04.LTS用にしかないので他のバージョンならRStudioサーバー使った方が早いかと
なお、パッチはRStudio 日本語とかでグクってみて
452:１３２人目の素数さん
18/06/23 17:38:43.39 aZCZP6wm.net
windows版のRStudioでもたまに日本語入力ができなくなる。
別のソフトでIMEのON/OFFをやってから戻ると直るけど。
453:１３２人目の素数さん
18/06/23 20:18:55.21 4XxgOFIA.net
>>440
うちも同じ。
454:１３２人目の素数さん
18/06/23 21:30:47.03 JcH4CpHc.net
Windowsはストアアプリでも同じ症状でるから持病なんじゃないかと思う
Ubuntuはfcitxのパッチ当てとけばWindowsのような症状は出ないので快適
455:１３２人目の素数さん
18/06/25 13:23:09.19 rmZ6BvQE.net
これをやってみた。
URLﾘﾝｸ(img-cdn.jg.jugem.jp)
スクリプトは
URLﾘﾝｸ(imagizer.imageshack.com)
456:0/t1VcqA.jpg 小数点以下500桁でみると19個素数があった。 > vip=Vectorize(is.prime) > (idx=which (vip(N))) [1] 100 124 150 172 183 202 215 219 255 296 310 314 323 346 400 417 439 441 474 > cat(substring (e,idx[1],idx[1]+9)) 7427466391
457:１３２人目の素数さん
18/06/26 07:08:45.35 Na/Ih9Bj.net
問題
９９人の囚人がいます。彼らの頭に１～１００までのナンバーカードが貼りつけられた帽子をランダムにかぶせます。
他人の帽子は見ることができても、自分の帽子は見ることができません。
帽子の数は全部で１００なので、一つ使われずに余ります。
そのナンバーは囚人達にはわからないようにしておきます。
この状況で、囚人たちに一斉に自分のナンバーを宣言させて、全員が正解だったら釈放するという賭けをします。
囚人たちには帽子をかぶせられる前に相談タイムが設けられています。
どういう戦略を取れば、助かる確率を最も高くできるでしょうか？
URLﾘﾝｸ(000013.blogspot.com)
解答を読んでも数理が理解できないが
確率が0.5になるのはシミュレーションできた。
# URLﾘﾝｸ(000013.blogspot.com)
inversion <- function(x){
n=length(x)
ret=numeric(n)
for(i in 1:(n-1)){
ret[i] = sum(x[i] > x[(i+1):n])
}
sum(ret) # inversion number
}
is.even= function(x) !inversion(x)%%2 # is inverion number even?
prisoner99 <- function(n=100){
indx=sample(1:n,1) # defective number
X=sample((1:n)[-indx])
Y=numeric(n-1)
for (i in 1:(n-1)){ # select as even permutation
x1=X[-i]
x2=(1:n)[!(1:n) %in% x1] # two numbers unseen for i-th prisoner
tmp=X
tmp[i]=x2[1] ; tmp[n]=x2[2]
Y[i]=ifelse(is.even(tmp), x2[1],x2[2])
}
all(X==Y)
}
mean(replicate(1e3,prisoner99()))
458:１３２人目の素数さん
18/07/11 16:20:29.88 zFHa28EV.net
>>438
遅レスだが、Ubuntu 18.04LTSならデフォルトのIBus-mozcで日本語入力ができるそうだ。
詳しくはUbuntu Weekly Recipeの第527回を読んでくだされ
459:１３２人目の素数さん
18/07/22 11:22:38.35 Ott8rTSz.net
覆面算 RED + WHITE = COLOR
どうもＣでやったらうまくいかないのでＲでやってみた。
# RED + WHITE = COLOR
x=c('R','E','D','W','H','I','T','E','C','O','L','O','R')
unique(x)
redwhite <- function(R,E,D,W,H,I,T,C,O,L){
red=10^(2:0)
white=10^(4:0)
color=10^(4:0)
sum(red*c(R,E,D))+sum(white*c(W,H,I,T,E)) - sum(color*c(C,O,L,O,R))
}
x=unique(x) ; x
REDWHITECOLOR <- function(x){
R=x[1]
E=x[2]
D=x[3]
W=x[4]
H=x[5]
I=x[6]
T=x[7]
C=x[8]
O=x[9]
L=x[10]
cat(paste('RED = ',R,E,D,
' WHITE = ',W,H,I,T,E,
' COLOR = ',C,O,L,O,R),'\n')
}
REDWHITECOLOR(unique(x))
library(gtools)
perm=permutations(n=10,r=10,v=0:9)
perm=perm[perm[,1]!=0&perm[,4]!=0&perm[,8]!=0,] # R!=0,W!=0,C!=0
head(perm) ; tail(perm)
n=nrow(perm)
re=numeric(n)
for(i in 1:n){
re[i]=redwhite(perm[i,1],perm[i,2],perm[i,3],perm[i,4],perm[i,5],perm[i,6],perm[i,7],perm[i,8],perm[i,9],perm[i,10])
}
hist(re)
(indx=which(re==0))
REDWHITECOLOR(perm[indx,])
> REDWHITECOLOR(perm[indx,])
RED = 5 8 7 WHITE = 3 9 6 1 8 COLOR = 4
460: 0 2 0 5
461:１３２人目の素数さん
18/07/22 12:50:02.38 Ott8rTSz.net
# AAB
# × CD
# ------
# CCE
# BFAA
# ------
# BEBDE
f <- function(A,B,C,D,E,F){
(A*100+A*10+B)*D==C*100+C*10+E &(A*100+A*10+B)*C==B*1000+F*100+A*10+A &
(C*100+C*10+E)+(B*1000+F*100+A*10+A)*10==B*10000+E*1000+B*100+D*10+E
}
library(gtools)
perm=permutations(n=10,r=6,v=0:9)
n=nrow(perm)
re=numeric(n)
for(i in 1:n){
re[i]=f(perm[i,1],perm[i,2],perm[i,3],perm[i,4],perm[i,5],perm[i,6])
}
indx=which(re==TRUE)
perm[indx,]
462:１３２人目の素数さん
18/07/23 11:09:01.32 SDW/E7TY.net
大変初歩的な質問で恐縮ですが、ご教示いただければ幸いです。
PERT分布を学んでいるのですが、ベータ分布のスケーリングの仕方が分かりません。
最小値a、最頻値ｂ、最大値ｃ として、a, b, c から求めたパラメータをa1、a2 とします。
このとき、
dbeta(x, a1, a2)*(c-a)+a
とすると、当然ながら横軸の値の範囲は0から1で、縦軸の値のみ大きくなってしまいます。
横軸の値の範囲がaからcで、縦軸の値は確率密度のままにするには、どうしたらよいのでしょうか？
463:１３２人目の素数さん
18/07/23 12:15:25.71 llIMc79p.net
え？aだけ右に動かしたいならこれでいいのでは？
dbeta(x-a, a1, a2)
464:１３２人目の素数さん
18/07/23 18:17:20.65 SDW/E7TY.net
>>449
447です。ご教示どうもありがとうございました。
a <- 40 # 最小値
b <- 45 # 最頻値
c <- 70 # 最大値
mu <- (a+4*b+c)/6 # PERT分布の平均値
a1 <- b*(mu-a)/(c-a) # パラメータ a1
a2 <- b*(c-mu)/(c-a) # パラメータ a2
pert1 <- function(x) dbeta(x,a1,a2) # ベータ分布布の計算
まではできたのですが、ご教示いただいた方法で
pert1 <- function(x) dbeta(x-a ,a1,a2)
とすると、ｘ の範囲が0から1で、ｙが0の一様分布のようなグラフになってしまいます。
465:１３２人目の素数さん
18/07/23 20:05:23.59 5MGL8f5L.net
4パラべーたへの変換は下記かな？
URLﾘﾝｸ(en.wikipedia.org)
a1とa2のパラメータ推定って合ってるんだっけ？僕も正直わかんない
やりたいことはこれだと思う
install.packages("mc2d")してからやってね
----------------------------------------------------
rm(list=ls())
library(mc2d)
a <- 40 # 最小値
b <- 45 # 最頻値
c <- 70 # 最大値
pert4<-function(x) dpert(x,min=a,mode=b,max=c,shape=4)
x1<-40:70
plot(x1,pert4(x1))
466:１３２人目の素数さん
18/07/23 21:05:56.52 +HoGXx9d.net
>>451
447です。
450様、まさにこれ、これです。どうもありがとうございました。
おかげさまで大変助かりました。
448様にも、ご助言下さいましたこｔ、厚く御礼申し上げます。
467:１３２人目の素数さん
18/07/27 15:50:12.19 E7dz8jqT.net
lmtestのパッケージをインストールしようとすると不正なマルチバイト文字がありますと出てくるのですがどうすればいいですか？
468:１３２人目の素数さん
18/07/27 17:33:01.36 ufA3ZDTu.net
>>453
Windows 10 （Rの日本語ヘルプはインストールしていない）+ R studioだけど、そのメッセージ出ないな。
469:１３２人目の素数さん
18/07/27 23:29:04.58 klPq25x/.net
ユーザー名が日本語かな
470:１３２人目の素数さん
18/07/28 12:56:05.38 lhV9Awbf.net
職場のSPSS信者たちにRを認めてもらうにはどうしたら良いだろう
「安かろう悪かろう」「俺が知らないものは三流」的な考え方の人が多くて肩身が狭い…
471:１３２人目の素数さん
18/07/28 17:53:14.05 XHIBT/Gx.net
>>456
>>453のようなエラーが出たら「無料のものはやっぱりダメね�
472:v 有償(SPSS) →信頼できる ボランティア→いい加減 と言う思考らしい
473:１３２人目の素数さん
18/07/28 18:33:45.13 l8Bhwsl9.net
>>455
ユーザー名はアルファベット1文字です
studioインストールしてないんですけど、しないとダメですか？
474:１３２人目の素数さん
18/07/28 21:48:13.01 hNHqu0EH.net
他のパッケージはインストールできるの？
475:１３２人目の素数さん
18/07/29 07:46:06.57 u2P99O/7.net
linuxなら環境変数LANGをCにすれば解決する問題だからwindowsでも似たような対処すれば行けそう。
476:１３２人目の素数さん
18/07/29 19:38:42.16 hjP1s1BE.net
これかな？
> Sys.setlocale(locale="C") # locale の変更（USの設定にしないと表示が変になる）
[1] "C"
元に戻すには
> Sys.setlocale(locale="Japanese_Japan.932")
477:１３２人目の素数さん
18/08/02 20:59:18.66 DvtHDake.net
すみません、教えて下さい。
フィッシャーの正確確率検定で、
適合度（一様性）の検定ってできますか？
できるのは独立性の検定のみ？
Rでいろいろ試したみたんですが、
2群ないとエラーになる。。
478:１３２人目の素数さん
18/08/02 21:16:50.41 GWs8Hfcd.net
適合度ならカイ二乗検定でやればいいのでは？
479:１３２人目の素数さん
18/08/02 21:52:11.65 DvtHDake.net
>>463
レス感謝です。
そうなんですが、
期待度数が5未満になってしまった場合は、
カイ二乗検定は向かないんですよね？
480:１３２人目の素数さん
18/08/02 21:54:12.25 82PNyNf0.net
それで他の方法を探しています
481:１３２人目の素数さん
18/08/02 23:02:57.38 clch7kxN.net
>>464
こういうので
r1=5;r2=4;n1=10;n2=12
prop.test(c(r1,r2),c(n1,n2),correct=TRUE)
chisq.test(matrix(c(r1,r2,n1-r1,n2-r2),2,byrow=TRUE),correct=TRUE)
警告がでるという話かな？
Warning message:
In chisq.test(matrix(c(r1, r2, n1 - r1, n2 - r2), 2, byrow = TRUE), :
Chi-squared approximation may be incorrect
482:１３２人目の素数さん
18/08/03 09:04:31.58 VGe+pklE.net
>>464
適合度を見るのに関係ないだろ。
> cnt
[1] 1 0 2 4 3 4 12 6 12 11 21 22 21 37 40 30 59 44 49 47 38 48 36 43 38
[26] 46 33 34 30 21 26 33 17 13 14 12 8 16 7 7 9 6 6 5 3 8 4 4 4 1
[51] 0 2 1 0 0 0 0 2
> prb
[1] 0.0002947255 0.0006464052 0.0012538873 0.0022037831 0.0035720843
[6] 0.0054114364 0.0077413658 0.0105430123 0.0137588816 0.0172972083
[11] 0.0210398698 0.0248524558 0.0285950600 0.0321325300 0.0353432214
[16] 0.0381256526 0.0404028045 0.0421240952 0.0432652773 0.0438266419
[21] 0.0438299769 0.0433147339 0.0423338191 0.0409493611 0.0392287274
[26] 0.0372409836 0.0350539129 0.0327316483 0.0303329196 0.0279098757
[31] 0.0255074215 0.0231629878 0.0209066521 0.0187615256 0.0167443293
[36] 0.0148660914 0.0131329065 0.0115467110 0.0101060386 0.0088067297
[41] 0.0076425765 0.0066058952 0.0056880193 0.0048797146 0.0041715201
[46] 0.0035540195 0.0030180499 0.0025548579 0.0021562071 0.0018144483
[51] 0.0015225575 0.0012741483 0.0010634658 0.0008853653 0.0007352804
[56] 0.0006091856 0.0005035529 0.0004153084
> chisq.test(cnt, p=prb, rescale.p=TRUE, simulat
483:e.p.value=TRUE) Chi-squared test for given probabilities with simulated p-value (based on 2000 replicates) data: cnt X-squared = 65.504, df = NA, p-value = 0.2189 5未満の数があっても関係がない
484:１３２人目の素数さん
18/08/03 10:29:47.68 jRZN4rrt.net
>>467
警告が出るのは独立性検定の場合でしたね。
485:１３２人目の素数さん
18/08/04 02:59:15.87 8dNre62F.net
463です。
カイ二乗検定で、期待度が5未満が多いのは不適切なのかと思い込んでおりましたが、
それは独立性の検定の場合のことであり、
適合度（一様性）の検定の場合は、
気にしなくて良かったのか。。
ならカイ二乗検定でやります。
いろいろ教えて下さり、ありがとうございます。
486:１３２人目の素数さん
18/08/04 19:49:50.94 3IVOmTrE.net
なんでお前らSPSS使わないの？
URLﾘﾝｸ(www.cdleidyba.lt)
487:１３２人目の素数さん
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
普及度、パッケージ数、専門度、文献数の点でそれはない