【R言語】統計解析フリーソフトR 第6章【GNU R】 [無断転載禁止]©2ch.net
■ このスレッドは過去ログ倉庫に格納されています
Run any R code you like. There are over three thousand R packages preloaded. https://rdrr.io/snippets/ Fisher test検定時に p<2.2e-16 と表示されるんですが、これより小さい値の指数桁数を正確に表記する方法教えて下さい。 例えば5.8e-35となるようにです。 fisher.test関数の返り値はlist型で、その中にp.valueという名前でp値が格納されているから$演算子を使って直接参照するか、broom::tidy関数に返り値を渡してdata.frame形式で出力してやれば見れる >>352 >>354 が言う直接的な参照 > fisher.test(matrix(c(1,120,130,2),2))$p.value [1] 1.691912e-69 352です。よく分かりました! ありがとうございます! >>356 技術的な助言をしたけど、学術的に言えば、 p < 0.01 は全て p < 0.01 として、具体的なp値を考える意味はないと思うよ。 一部の例外的な研究分野を除いて(e.g., 遺伝統計学)。 はい、まさにその例外的な分野で使おうとしてます。ありがとうございます。 Rのガンマ関数はいくつでオーバーフローするかやってみた。 > i=1 > while(gamma(i)!=Inf){ + i=i+1 + } Warning message: In gamma(i) : value out of range in 'gammafn' > i [1] 172 > gamma(172) [1] Inf Warning message: value out of range in 'gammafn' > gamma(171) [1] 7.257416e+306 matplot()で折れ線グラフ描いたときに、X軸をカテゴリで示したいのですが、 可能でしょうか? 例えばtemp <- c("0時間","8時間","24時間","48時間")として、 matplot()の引数にtempをとるやり方です。 他にもやり方あれば教えてください。 >>371 matplot(..., xaxt="n") axis(1, at=seq(along=temp), lab=temp) >>372 遅くなりましたがありがとうございました。 できました! 特定の長方形の中に複数の長方形を最小面積で敷き詰める平面充填に関するパッケージってありませんかね # jonckheereテストを書いてみた jonckheere <- function(L, alternative = c("two.sided", "increasing", "decreasing"), cat=TRUE){ # L : list of vectors A1,A2,...,Ak, with assumed tendency How.Many.Greater.Pairs <- function(A,B){ # How many pairs of A[i] > B[j], count as 0.5 when equal, n.a = length(A) n.b = length(B) how.many.greater.pairs = 0 for(i in 1:n.a){ for(j in 1:n.b){ how.many.greater.pairs = how.many.greater.pairs+ifelse(A[i]==B[j],0.5,A[i]>B[j]) } } return(how.many.greater.pairs) } Sum.of.Greater.Pairs <- function(L){ #L=list(A1,,,,Ak),A1 < A2 < A3,..,< Ak : vector k = length(L) comb = combn(1:k,2) # possible combinaition of pairs to compare n.comb = ncol(comb) # how many combinations J = 0 # sum of greater pairs for(i in 1:n.comb){ J = J + How.Many.Greater.Pairs(L[[comb[1,i]]],L[[comb[2,i]]]) } return(J) } J = Sum.of.Greater.Pairs(L) n = sapply(L,length) N = sum(n) EJ = (N^2-sum(n^2))/4 VJ = (N^2*(2*N+3)-sum(n^2*(2*n+3)))/72 Z = (J-EJ)/sqrt(VJ) alternative = match.arg(alternative) p.value = switch(alternative, 'two.sided' = 2 * min(pnorm(Z), pnorm(-Z), 0.5), 'increasing' = pnorm(Z), 'decreasing' = pnorm(-Z)) if(cat){ cat( 'p.value = ', p.value,'\n') cat('alternative hypothesis: ' ,alternative,'\n') } invisible(p.value) } >>388 実態調査か何か? <-と=は挙動が違う場合があるので、使い分けていますが、 代入はどっちかと問われたら、無論 <- または -> なお、 > 1 -> x これはエラーにならないけど、 > 1 = x 1 = x でエラー: 代入の左辺が不正 (do_set) です これはエラー >>389 俺は基本=派。 関数の定義は z.test <- function(x,n=16,sigma=1){ z=sqrt(n)*mean(x)/sigma 2*pnorm(abs(z),lower=FALSE) } と書いている。 >1 = x でエラー 当たり前 そんな使い方なんてするかよ 他言語と同じく=一文字の方がすっきりしてイイ こういうのが紛らわしいから、俺は = 推奨。 x <- 1 if(x <- 1) print('YES') if(x < -1) print('YES') やったことなかったので関数の初期値設定に<-を使うとどうなるかやってみた。 まず、= の場合 > z.test <- function(x,n=16,sigma=1){ + z=sqrt(n)*mean(x)/sigma + 2*pnorm(abs(z),lower=FALSE) + } > z.test(1:3) [1] 1.244192e-15 <- で初期値設定すると、エラー > z0.test <- function(x,n<-16,sigma<-1){ Error: unexpected assignment in "z0.test <- function(x,n<-" > z=sqrt(n)*mean(x)/sigma Error in mean(x) : object 'x' not found > 2*pnorm(abs(z),lower=FALSE) Error in pnorm(abs(z), lower = FALSE) : object 'z' not found > } 俺は見栄えがいいと思って関数定義には<-を使っているけど = でも通常に動作する。 > z.test = function(x,n=16,sigma=1){ + z=sqrt(n)*mean(x)/sigma + 2*pnorm(abs(z),lower=FALSE) + } > z.test(1:3) [1] 1.244192e-15 <- 推奨の人に聞きたいのだけど <- でないと動作しないってことあるのだろうか? RStudioつかってりゃ[Alt+-]で簡単入力 それと引数の指定は代入じゃないと思うんだが感覚が違うよかな? >>396 essだとアンダースコアを入れるを勝手に「<-」になる。 >>394 それは代入ではなく、関数の規定値の設定。 規定値の設定は「=」と決められているので、「<-」はアウト。 ただし、関数を実行するときには「<-」を使うことができる。 > mean(x<-1:10) [1] 5.5 > x [1] 1 2 3 4 5 6 7 8 9 10 =に拘り続けて、<-に変える瞬間の君たちがたまらないよ >>398 ありがとうございます。 そんな使い方ができたのですね。 こんなのしか知りませんでした。 > f <- function() { + x<<-1:10 + mean (x) + } > f() [1] 5.5 > x [1] 1 2 3 4 5 6 7 8 9 10 一画面しかグラフがないのに Hit <Return> to see next plot: と出るの鬱陶しいな、と思っていたら par(ask=FALSE) と設定しておけばいいんだな。 同名のdataがあるのでパッケージを::で指定してもうまくいかなかった。 data(netmeta::parkinson) でなくて data(parkinson, package = 'netmeta') とするんだな。 hist(rnorm(100),col=rgb(runif(1),runif(1),runif(1),runif(1))) だと、1色だけど hist(rnorm(100),col=sample(colours(),(sample(1:10,1)))) なら1〜10色で表示される。 > hist(rnorm(100),col=rgb(runif(1),runif(1),runif(1),runif(1))) 何色に増やしてもいいがな。 何をやりたいのかな? hist(rnorm(5000),col=apply(matrix(runif(80),4), 2, function(x){rgb(x[1],x[2],x[3],x[4])})) :::てどういう時に使うんだろ? ソース読みたくて library(BayesFactor) ttestBF_indepSample では表示されなかったが、 library(BayesFactor) BayesFactor:::ttestBF_indepSample だと出てきた。 >>410 グダグダ言っていないで、さっさとヘルプを参照すればいいじゃないか? ?':::' >>411 ありがとう。 stats:::t.test.default でt検定のソースが見られた 関数を値で上書きしてしまった場合に元に戻す方法を教えて頂きたく。 例) rm <- 1 としてしまった場合にrm関数を元に戻したいです。 よろしくお願いします。 >>413 直後ならRstudioでctrl+zで元に戻す。 再起動してたらこの方法は無理だろな >>413 自分で作ったのではない関数の rm だったら rm(rm) で数値ベクトルオブジェクトが消えて、関数の rm が見えるようになる。 >>413 > 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()で大丈夫だろう >>413 既に説明されてるけど上書きじゃなくて別オブジェクト扱いになるからrmで消すか明示的にパッケージを::演算子で関数を指定してやればおけ 自作関数だと上書きされちゃうけど 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) 同じ名前空間にあるなら上書きされちゃうよ .Gloalenv(だったかな?)に自作関数があり、そこに上書きして変数にしてるから消したら消えるだけで戻せない ちょっと古いけど、この辺りを読むとなぜそうなるかが分かると思う (Rの)環境問題について その1。 ttps://qiita.com/kohske/items/325bdf48f4f4885a86f1 >>418 戻せます。 > 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), [以下略] >>419 >.Gloalenv case sensitiveなRでそれはないのでは? > help(".GlobalEnv") >>414-422 みなさまご返信ありがとうございました。 例示したものがbase::rmだったので自作関数のことまでは想像だにしませんでした…(がとても良い勉強になりました。) 自環境でも上書きした関数を元に戻すことが出来ることを確認しました。 >>421 戻したいのが自作関数 rm = function(x) asin(sqrt(x)) だったら上書きされて無理じゃないの? ゴミ箱から消去したファイルを復活させるソフトもあるから PCに詳しければ可能かもしれないけれど Rの機能としては復活させるのは無理と思う。 自作のrm関数はbaseパッケージとは別に定義されるから base::rm(rm)で自作の方のrmは消えてbase::rmが「復活」します。 ああすまん 自作の方を上書きした場合は、そっちの復活は無理ですね。 >>425 ちょっと意味がわからない。 自作の関数なら定義を再実行するだけなのに、 消える消えないって何が問題なんだ? >>424 関数だけでなく既定値でも同じ。 > pi [1] 3.141593 > pi=2 > pi [1] 2 > rm(pi) > pi [1] 3.14159 >>306 知恵袋に漸化式があったので計算スクリプトを書いてみた。 表がでる確率pのコインを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(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] } >>306 この例で、答えとして何が返ってきて欲しいのか? 2で合ってる? >>433 1が5回以上出現しているからTRUE(数字なら1)が返り値。 >>434 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があるか? } >>435 whichとdiffを用いての解法のお返事ありがとうございます。 実行時間を>302のスクリプトと比べてみました。 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)回コインをなげてn(=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の使い方が勉強になりました。 時間を割いてスクリプトを作成していただいてありがとうございます。<(_ _)> Rstudioで日本語の入力ができないのですが、解決する方法知ってる方いませんか? 具体的には、IME自体をオンにすることができません。 日本語の表示やコピペは問題なくできます。 環境 Linux - Utubntu RStudio Version 1.0.143 RStudioのQtが古いのでパッチ当てないと日本語入力できないのが現状ですわ パッチはUbuntu 16.04.LTS用にしかないので他のバージョンならRStudioサーバー使った方が早いかと なお、パッチはRStudio 日本語とかでグクってみて windows版のRStudioでもたまに日本語入力ができなくなる。 別のソフトでIMEのON/OFFをやってから戻ると直るけど。 Windowsはストアアプリでも同じ症状でるから持病なんじゃないかと思う Ubuntuはfcitxのパッチ当てとけばWindowsのような症状は出ないので快適 これをやってみた。 http://img-cdn.jg.jugem.jp/d75/1405116/20140609_1105745.jpg スクリプトは http://imagizer.imageshack.com/img921/8020/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 問題 99人の囚人がいます。彼らの頭に1〜100までのナンバーカードが貼りつけられた帽子をランダムにかぶせます。 他人の帽子は見ることができても、自分の帽子は見ることができません。 帽子の数は全部で100なので、一つ使われずに余ります。 そのナンバーは囚人達にはわからないようにしておきます。 この状況で、囚人たちに一斉に自分のナンバーを宣言させて、全員が正解だったら釈放するという賭けをします。 囚人たちには帽子をかぶせられる前に相談タイムが設けられています。 どういう戦略を取れば、助かる確率を最も高くできるでしょうか? http://000013.blogspot.com/2010/12/99.html 解答を読んでも数理が理解できないが 確率が0.5になるのはシミュレーションできた。 # http://000013.blogspot.com/2010/12/99.html 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())) >>437 遅レスだが、Ubuntu 18.04LTSならデフォルトのIBus-mozcで日本語入力ができるそうだ。 詳しくはUbuntu Weekly Recipeの第527回を読んでくだされ 覆面算 RED + WHITE = COLOR どうもCでやったらうまくいかないのでRでやってみた。 # 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 0 2 0 5 # 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,] 大変初歩的な質問で恐縮ですが、ご教示いただければ幸いです。 PERT分布を学んでいるのですが、ベータ分布のスケーリングの仕方が分かりません。 最小値a、最頻値b、最大値c として、a, b, c から求めたパラメータをa1、a2 とします。 このとき、 dbeta(x, a1, a2)*(c-a)+a とすると、当然ながら横軸の値の範囲は0から1で、縦軸の値のみ大きくなってしまいます。 横軸の値の範囲がaからcで、縦軸の値は確率密度のままにするには、どうしたらよいのでしょうか? え?aだけ右に動かしたいならこれでいいのでは? dbeta(x-a, a1, a2) >>448 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) とすると、x の範囲が0から1で、yが0の一様分布のようなグラフになってしまいます。 ■ このスレッドは過去ログ倉庫に格納されています
read.cgi ver 07.5.5 2024/06/08 Walang Kapalit ★ | Donguri System Team 5ちゃんねる