【R言語】統計解析フリーソフトR 第6章【GNU R】 [無断転載禁止]©2ch.net
■ このスレッドは過去ログ倉庫に格納されています
>>126 125です。レス遅くなったけど、やりたいことができました。ありがとう。 array(1:48, dim=c(6, 4, 2)) を,3行2列おきに平均して ,,1 [,1] [,2] [1,] 5 17 [2,] 8 20 ,,2 [,1] [,2] [1,] 29 41 [2,] 32 44 としたいんですが,効率的な方法はないでしょうか。 実際のデータは,10000x10000x400なんです。 >>154 >3行2列 さっぱりわからん。 例えば、結果の[1,1,1]の5は、どれとどれの平均? >>155 6行4列の左上の6つ1,2,3,7,8,9の平均5だと思う。 >>154 効率的かどうかはわからんが、やってみた (d=array(1:48, dim=c(6, 4, 2))) a=3 b=2 l=6 m=4 n=2 re=NULL for(k in 1:n){ for(j in 1:(m/b)){ for(i in 1:(l/a)){ re=append(re,mean(d[1:3+3*(i-1),1:2+2*(j-1),k])) } } } array(re,dim=c(m/b,l/a,n)) >>155 >>154 です。 説明が悪くてすみません。 >>156 の通りです。 >>157 ありがとうございます。 手元にデータがないので明日試してみます。 なるほど、じゃあその5,8,17,20はこういうことか。 > a <- array(1:48, dim=c(6, 4, 2)) > b <- expand.grid(1:(dim(a)[1] %/% 3), 1:(dim(a)[2] %/% 2)) > f <- function(i, j, k=1){mean(a[(3*i-2):(3*i),(2*j-1):(2*j),k])} > mapply(f, b[,1], b[,2]) [1] 5 8 17 20 あぁ、すでに >>157 が回答を書いていたorz > d=array(1:48, dim=c(6, 4, 2)) > a=3 > b=2 > l=6 > m=4 > n=2 > re=NULL > for(k in 1:n){ + for(j in 1:(m/b)){ + for(i in 1:(l/a)){ + re=append(re,mean(d[1:a+a*(i-1),1:b+b*(j-1),k])) + } + } + } > array(re,dim=c(m/b,l/a,n)) , , 1 [,1] [,2] [1,] 5 17 [2,] 8 20 , , 2 [,1] [,2] [1,] 29 41 [2,] 32 44 データは10000x10000x400なら l=10000 m=10000 n=400 でよくね? 最初の1行で > d=array(1:(10000*10000*400)) Error: cannot allocate vector of size 298.0 Gb エラーがでたw 正しいエラーwはこっち > d=array(1:(10000*10000*400),dim=c(6, 4, 2)) Error: cannot allocate vector of size 298.0 Gb そもそも、10000x10000x400だと10000行は3行にならないな。最終行を無視して9999行にするのかな? >160を参考に>157のforだらけのスパゲティー・ソースを(ちょっと一般化して)改変しました。 (変数名の重複を避けるため>160の変数名は変更してあります。) > a=3 > b=2 > c=1 > > l=6 > m=4 > n=2 > > A <- array(1:48, dim=c(l,m,n)) > B <- expand.grid(1:(l %/% a), 1:(m %/% b),1:(n %/% c)) > f <- function(i, j, k){mean(A[(a*(i-1)+1):(a*i),(b*(j-1)+1):(b*j),(c*(k-1)+1):(c*k)])} > re=mapply(f, B[,1], B[,2],B[,3]) > array(re,dim=c(m%/%b,l%/%a,n%/%c)) , , 1 [,1] [,2] [1,] 5 17 [2,] 8 20 , , 2 [,1] [,2] [1,] 29 41 [2,] 32 44 > 上級者のコードを読むのは勉強になります。(深謝) 先輩方,>>154 です。 巨細にわたってご教示・ご議論いただきありがとうございました。 >>166 のコードを参考にして,無事処理することができました。 実行速度も申し分なく,効率的に作業が進められそうです。 今回は私には目から鱗な解法で,勉強になりました。 質問して良かったと思っています。。本当にありがとうございました。 10,000行の扱いですが,3行ごとに処理して1行無視してもよいし, 4から7行程度で処理してもよく,フレキシブルです。 今回は4行にして余りなしにしました。 処理速度を比べてみた。 fs<-function(){ #スパゲッティ・ソース d=array(1:(l*m*n),dim=c(l,m,n)) re=NULL for(k in 1:(n/c)){ for(j in 1:(m/b)){ for(i in 1:(l/a)){ re=append(re,mean(d[1:a+a*(i-1),1:b+b*(j-1),1:c+c*(k-1)])) }}} array(re,dim=c(m/b,l/a,n)) } fe <- function(){ # Expertソース A <- array(1:(l*m*n), dim=c(l,m,n)) B <- expand.grid(1:(l %/% a), 1:(m %/% b),1:(n %/% c)) f <- function(i, j, k){mean(A[(a*(i-1)+1):(a*i),(b*(j-1)+1):(b*j),(c*(k-1)+1):(c*k)])} re=mapply(f, B[,1], B[,2],B[,3]) array(re,dim=c(m%/%b,l%/%a,n%/%c)) } a=4 b=2 c=1 l=400 m=200 n=2 > system.time(fs()) user system elapsed 2.17 0.00 2.18 > system.time(fe()) user system elapsed 1.14 0.00 1.14 >160氏のスキルに改めて感服。 # 既約剰余類群 (Z/nZ)* kjrg <- function(n){ nn=0:(n-1) f=function(x,y) (x*y)%%n names(nn)=paste0('n',nn) z=outer(nn,nn,f) idx=which(gmp::gcd(1:(n-1),n)==1)+1 return(z[idx,idx]) } kjrg(10) kjrg(15) kjrg(24) > kjrg(24) n1 n5 n7 n11 n13 n17 n19 n23 n1 1 5 7 11 13 17 19 23 n5 5 1 11 7 17 13 23 19 n7 7 11 1 5 19 23 13 17 n11 11 7 5 1 23 19 17 13 n13 13 17 19 23 1 5 7 11 n17 17 13 23 19 5 1 11 7 n19 19 23 13 17 7 11 1 5 n23 23 19 17 13 11 7 5 1 # 既約剰余類群での加減乗除 n=24 n_star=which(gmp::gcd(1:(n-1),n)==1) ; n_star tasu <- function(x,y) (x+y)%%n hiku <- function(x,y) (x-y)%%n # row - col kake <- function(x,y) (x*y)%%n g=function(x) n_star[which(x==1)] .M=outer(n_star,n_star,kake) G=apply(.M,2,g) gyaku <- function(x) n_star[which(G==(x%%n))] waru <- function(x,y) (x*gyaku(y))%%n # row ÷ col xx=yy=c(0,n_star) names(xx)=paste0('x',c(0,n_star)) names(yy)=paste0('y',c(0,n_star)) outer(xx,yy,tasu) # x + y outer(xx,yy,hiku) # x - y outer(xx,yy,kake) # x * y X=Y=n_star #outer(X,Y,waru) # WRONG!! a=expand.grid(X,Y) b=matrix(mapply(waru,a[,1],a[,2]),length(X)) rownames(b)=paste0('x',n_star) colnames(b)=paste0('y',n_star) b # x / y 名前も一致してidenticalなんだな。 > x=1:10 > y=1:10 > identical(x,y) [1] TRUE > names(x)=LETTERS[1:10] > identical(x,y) [1] FALSE > names(y)=LETTERS[1:10] > identical(x,y) [1] TRUE > names(y)=letters[1:10] > identical(x,y) [1] FALSE allの方は内容の全体が一致していれば順番は入れ替わってもTRUEなんだな。 > x=1:10 > y=1:10 > all(x,y) [1] TRUE > names(x)=LETTERS[1:10] > all(x,y) [1] TRUE > names(y)=letters[1:10] > all(x,y) [1] TRUE > all(x,sample(y)) [1] TRUE beki <- function(x,n,p){ # x^n%%p if(n==0) return(1) if(n==1) return(x%%p) re=numeric(n) re[1]=x for(i in 1:(n-1)) re[i+1] = (x*re[i])%%p return(re[n]) } p=1009 # フェルマーの小定理 xx=0:(p-1) n=p-1 y=sapply(xx,function(x)beki(x,n,p)) summary(y) 90歳女性のお看取り >120のコードから女性の平均余命をグラフにしてみた Age=85:105 (Life_Expectancy=sapply(Age,function(x) LE(f,x))) plot(Age,Life_Expectancy,pch=19,col=sample(colors())) abline(h=5,col='gray',lty=2) http://i.imgur.com/38x50fw.png 平均余命が5年を切るのは92歳なので、老衰にせずに 主治医指定の病名で診断書を作成。 カルテ記載のしっかりした患者の診療は負担なくこなせて(・∀・)イイ!! outerって使えない場面があるよね? 剰余系での加減乗除での実際 > n=5 # prime number > nn=1:(n-1) > tasu <- function(x,y) (x+y)%%n > hiku <- function(x,y) (x-y)%%n # row - col > kake <- function(x,y) (x*y)%%n > g=function(x) nn[which(x==1)] > .M=outer(nn,nn,kake) > G=apply(.M,2,g) > gyaku <- function(x) nn[which(G==(x%%n))] > waru <- function(x,y) (x*gyaku(y))%%n # row / col > waru(3,2) [1] 4 > xx=yy=c(0,nn) > names(xx)=paste0('x',c(0,nn)) > names(yy)=paste0('y',c(0,nn)) > outer(xx,yy,tasu) # x + y y0 y1 y2 y3 y4 x0 0 1 2 3 4 x1 1 2 3 4 0 x2 2 3 4 0 1 x3 3 4 0 1 2 x4 4 0 1 2 3 > outer(xx,yy,hiku) # x - y y0 y1 y2 y3 y4 x0 0 4 3 2 1 x1 1 0 4 3 2 x2 2 1 0 4 3 x3 3 2 1 0 4 x4 4 3 2 1 0 > outer(xx,yy,kake) # x * y y0 y1 y2 y3 y4 x0 0 0 0 0 0 x1 0 1 2 3 4 x2 0 2 4 1 3 x3 0 3 1 4 2 x4 0 4 3 2 1 > X=Y=nn > outer(X,Y,waru) # WRONG!! [,1] [,2] [,3] [,4] [1,] 1 1 1 1 [2,] NA NA NA NA [3,] NA NA NA NA [4,] NA NA NA NA > a=expand.grid(X,Y) > b=matrix(mapply(waru,a[,1],a[,2]),length(X)) > rownames(b)=paste0('x',nn) > colnames(b)=paste0('y',nn) > b # x / y y1 y2 y3 y4 x1 1 3 2 4 x2 2 1 4 3 x3 3 4 1 2 x4 4 2 3 1 >>174 そりゃあ、outer に渡す関数がまずいんだろ。 outer に渡す関数は、完全にベクトル対応になっていなければならないんじゃないの? でなければ、わざわざ outer 使うこともなかろう。 >>175 レスありがとうございます。 自作関数から自作関数を呼び出しているのが原因かと思っていましたが、 ベクトル対応していないとはこういうことでしょうか? 加減乗は引数のどちらをベクトルにしても動作 > kake(1:4,1) [1] 1 2 3 4 > kake(1,1:4) [1] 1 2 3 4 > tasu(1:4,1) [1] 2 3 4 0 > tasu(1,1:4) [1] 2 3 4 0 > hiku(1:4,1) [1] 0 1 2 3 > hiku(1,1:4) [1] 0 4 3 2 除算は第2引数をベクトルにすると誤動作 > waru(1:4,1) [1] 1 2 3 4 > waru(4,1:4) [1] 4 1 >>176 outer から呼び出されたときの引数を見てみればわかるよ。たぶん1度しか呼ばれず、できあがりの行列の要素すべてを計算するための引数がベクトルで渡されているはずだから。 >>177 ありがとうございます。 outerのソースをみたら、こんな記述があったので、X,Yを同時にベクトルで与えても動作する関数でないとダメみたいでした。 FUN <- match.fun(FUN) Y <- rep(Y, rep.int(length(X), length(Y))) if (length(X)) X <- rep(X, times = ceiling(length(Y)/length(X))) FUN(X, Y, ...) # greatest common divisor gcd <- function(x,y) ifelse(x|y,ifelse(x*y,max(intersect((1:abs(x))[abs(x)%%(1:abs(x))==0],(1:abs(y))[abs(y)%%(1:abs(y))==0])),max(abs(x),abs(y))),NA) 非正整数対応すると長くなるなぁ > gcd(24,15) [1] 3 > gcd(16,15) [1] 1 > gcd(5,0) [1] 5 > gcd(-24,15) [1] 3 > gcd(-24,-15) [1] 3 > gcd(0,0) [1] NA ユークリッドの互除法をつかうとこんな感じ > GCD <- function(x,y){ + if(round(x)!=x | round(y)!=y) stop('Not integer') + if(!x&!y) return(NA) + a=max(c(abs(x),abs(y))) + b=min(c(abs(x),abs(y))) + if(!a*b) return(a) + r=integer() + r[1]=a + r[2]=b + i=1 + r[i+2]=r[i]%%r[i+1] + while(r[i+2]!=0 & r[i+2]!=1){ + i=i+1 + r[i+2]=r[i]%%r[i+1] + } + return(ifelse(r[i+2]==0,r[i+1],1)) + } > GCD(0,0) [1] NA > GCD(24,15) [1] 3 > GCD(16,15) [1] 1 > GCD(5,0) [1] 5 > GCD(24,-15) [1] 3 > GCD(-24,-15) [1] 3 > GCD(2.3,4) Error in GCD(2.3, 4) : Not integer #.n発r中の狙撃手が.N発狙撃するときの命中数を返す Go13 <- function(.N, .n, r, k=10^3){ # k:シミュレーション回数 f <-function(S,N=.N,n=.n){ y=c(rep(1,S),rep(0,N-S)) sum(sample(y,n)) } xx=r:.N SS=NULL for(i in 1:k){ x=sapply(xx,f) SS=c(SS,which(x==r)-1+r) } print(summary(SS)) invisible(SS) } こういう問題を考えるのは楽しいね。 ゴルゴ13は100発100中 ゴルゴ14は10発10中 ゴルゴ15は1発1中 とする。 各々10000発撃ったとき各ゴルゴの命中数の期待値はいくらか? ドツボ13は100発0中 ドツボ14は10発0中 ドツボ15は1発0中 とする。 各々10000発撃ったときドツボの命中数の期待値はいくらか? ベクトルは再利用されるんだな > c(1,2,1,2)==c(1,2) [1] TRUE TRUE TRUE TRUE > all(c(1,2,1,2)==c(1,2)) [1] TRUE > equal <- function(x,y) sum((x-y)^2)==0 & length(x)==length(y) > equal(c(1,2,1,2),c(1,2)) [1] FALSE > x= list(length=5) > x $length [1] 5 > y=list() > length(y)=4 > y [[1]] NULL [[2]] NULL [[3]] NULL [[4]] NULL > z=list(5) > z [[1]] [1] 5 > w = vector(mode="list" , length= 5) > w [[1]] NULL [[2]] NULL [[3]] NULL [[4]] NULL [[5]] NULL Brunner-Munzel検定って比較する両群に重なりがないとp値がNAになるんだが、 これってバグなのか? > x=1:5 > y=6:10 > lawstat::brunner.munzel.test(x,y) Brunner-Munzel Test data: x and y Brunner-Munzel Test Statistic = Inf, df = NaN, p-value = NA 95 percent confidence interval: NaN NaN sample estimates: P(X<Y)+.5*P(X=Y) 1 > x=1:6 > y=6:11 > lawstat::brunner.munzel.test(x,y) Brunner-Munzel Test data: x and y Brunner-Munzel Test Statistic = 24.749, df = 10, p-value = 2.651e-10 95 percent confidence interval: 0.9423463 1.0298759 sample estimates: P(X<Y)+.5*P(X=Y) 0.9861111 突然ですが質問があります。最近Rソフトを利用し始めた者のです。単刀直入に、 ボロノイ図を制作した場合に分割された面積の平均を出力するには、どうすればよいのでしょうか? 現状としては住所をグーグルマップの緯度経度に変換し、分割まではできています。 宜しければどなたかお力添えをいただけないでしょうか? >>198 SpatialPolygonsクラスまでできるといると仮定して、 各ポリゴンの中にareaスロットがあるから、 それを抽出して、その平均を取れば良い。 ボロノイ図がSpatialPolygonsクラスになっていないなら、 まずは変換。 ボロノイ図の各領域の面積の平均値って、単に 描画領域の総面積÷点数 じゃないの? 何か暗黙の前提があるの? Stan って 離散量をparameter 扱いできないので m ? categorical(θ) としても m の分布を出せない。 JAGSだと m 〜 dcat(θ) とすると mの分布がだせる。 この理解でいいかな? 以前、ボロノイ図について質問させていただいた初心者です。 質問後色々と試行錯誤したのですがdir.areaを出せず悩んでおります。 library(deldir) > library(ggmap) > df <- read.csv("C:/Users/***/Desktop/***/***.csv",header=T) > loc <- c(min(df$lon),min(df$lat),max(df$lon),max(df$lat)) > GM <-ggmap(get_map(location=loc,zoom=11,source="google")) > GP <- geom_point(data=df,aes(x=lon,y=lat,colour=factor(type)),size=3) > dd <- deldir(df$lon,df$lat) >GS <-geom_segment(data=dd$dirsgs,aes(x=x1,y=y1,xend=x2,yend=y2),size=0.5,linetype=1) > GM+GP+GS と入力しボロノイ図は出せているのですが、 summaryと入力しても function (object, ...) UseMethod("summary") <bytecode: 0x06b5eef0> <environment: namespace:base> と表示されます。どの様に入力または再入力すればdir.areaが表示させられるでしょうか? 私自身初めてで、全く分からない状態でRを使用しているので、宜しければ分かり易く教えて頂けると幸いです。 ■ このスレッドは過去ログ倉庫に格納されています
read.cgi ver 07.5.5 2024/06/08 Walang Kapalit ★ | Donguri System Team 5ちゃんねる