【R言語】統計解析フリーソフトＲ 第6章【GNU R】 [無断転載禁止]©2ch.net
■ このスレッドは過去ログ倉庫に格納されています
>>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
エラーがでたｗ 正しいエラーｗはこっち
> 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
除算は第２引数をベクトルにすると誤動作
> 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検定って比較する両群に重なりがないとｐ値がＮＡになるんだが、
これってバグなのか？
> 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(θ) としても ｍ の分布を出せない。
JAGSだと
m 〜 dcat(θ）
とすると ｍの分布がだせる。
この理解でいいかな？ 以前、ボロノイ図について質問させていただいた初心者です。
質問後色々と試行錯誤したのですが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を使用しているので、宜しければ分かり易く教えて頂けると幸いです。 >>242
検証用のデータがないと誰も追試しないと思う。
>154は模擬データがあったのでやる気になった。 ■ このスレッドは過去ログ倉庫に格納されています
