【R言語】統計解析フリーソフトR 第6章【GNU R】 [無断転載禁止]©2ch.net
■ このスレッドは過去ログ倉庫に格納されています
特定の長方形の中に複数の長方形を最小面積で敷き詰める平面充填に関するパッケージってありませんかね # 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の一様分布のようなグラフになってしまいます。 4パラべーたへの変換は下記かな? https://en.wikipedia.org/wiki/Beta_distribution#Four_parameters_2 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)) >>450 447です。 450様、まさにこれ、これです。どうもありがとうございました。 おかげさまで大変助かりました。 448様にも、ご助言下さいましたこt、厚く御礼申し上げます。 lmtestのパッケージをインストールしようとすると不正なマルチバイト文字がありますと出てくるのですがどうすればいいですか? >>452 Windows 10 (Rの日本語ヘルプはインストールしていない)+ R studioだけど、そのメッセージ出ないな。 職場のSPSS信者たちにRを認めてもらうにはどうしたら良いだろう 「安かろう悪かろう」「俺が知らないものは三流」的な考え方の人が多くて肩身が狭い… >>455 >>452 のようなエラーが出たら「無料のものはやっぱりダメね」 有償(SPSS) →信頼できる ボランティア→いい加減 と言う思考らしい >>454 ユーザー名はアルファベット1文字です studioインストールしてないんですけど、しないとダメですか? linuxなら環境変数LANGをCにすれば解決する問題だからwindowsでも似たような対処すれば行けそう。 これかな? > Sys.setlocale(locale="C") # locale の変更(USの設定にしないと表示が変になる) [1] "C" 元に戻すには > Sys.setlocale(locale="Japanese_Japan.932") すみません、教えて下さい。 フィッシャーの正確確率検定で、 適合度(一様性)の検定ってできますか? できるのは独立性の検定のみ? Rでいろいろ試したみたんですが、 2群ないとエラーになる。。 >>462 レス感謝です。 そうなんですが、 期待度数が5未満になってしまった場合は、 カイ二乗検定は向かないんですよね? >>463 こういうので 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 >>463 適合度を見るのに関係ないだろ。 > 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, simulate.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未満の数があっても関係がない >>466 警告が出るのは独立性検定の場合でしたね。 463です。 カイ二乗検定で、期待度が5未満が多いのは不適切なのかと思い込んでおりましたが、 それは独立性の検定の場合のことであり、 適合度(一様性)の検定の場合は、 気にしなくて良かったのか。。 ならカイ二乗検定でやります。 いろいろ教えて下さり、ありがとうございます。 >>469 いまSPSSを使う理由って何なの?煽りじゃなくマジで >>472 周囲のspssユーザを見ると「指導教員がspssしか使えないので、 自分もspssしか使えない、 新しい(学術系)ソフトは誰かに丁寧に根気よく指導してみらえないと無理」、 これの再生産。 ■ このスレッドは過去ログ倉庫に格納されています
read.cgi ver 07.5.5 2024/06/08 Walang Kapalit ★ | Donguri System Team 5ちゃんねる