裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

ダメ出し:順列の逆関数

2011年10月24日 | ブログラミング

http://d.hatena.ne.jp/ryamada/20110622

これはちょっとすごい(ひどい)ですね。
定義された InvSeq は order 関数の結果と同じになります。

s が長いと,悲惨な結果になります。下の例では 1 万倍以上の速度比。

> s <- sample(10000) # sample(1:10000, 10000)  などとする必要はない

> system.time(invs <- InvSeq(s))
   user  system elapsed
 22.910   3.197  25.971

> system.time(invs2 <- order(s)) # 一瞬で終わる
   user  system elapsed
  0.002   0.000   0.002

> all.equal(invs, invs2)
[1] TRUE

> 22910/2
[1] 11455

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:ぐるぐる回るために

2011年10月24日 | ブログラミング

http://d.hatena.ne.jp/ryamada/20111002 について

Npt <- 10000
n <- 3
X <- matrix(rnorm(Npt*n), Npt, n)
X <- X/sqrt(rowSums(X^2))
Y <- rowSums(X[,1:(n-1)]*X[,2:n])

Npt<-10000
n<-3
X<-matrix(rnorm(Npt*n),Npt,n)

X<-X/sqrt(apply(X^2,1,sum))
#apply(X^2,1,sum)

Y<-rep(0,Npt)

for(i in 1:(n-1)){
    Y<-Y+X[,i]*X[,i+1]
}

より,20 倍速い

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:与えられた条件から、推定する~信頼区間?~

2011年10月24日 | ブログラミング

相変わらず for がお好き

http://d.hatena.ne.jp/ryamada/20110928

p<- seq(0,1,0.00001)
pp<- rep(0,length(p))
for(i in 1:length(p)){
    pp[i]<- dbinom(4,15,p[i])
}

より,35 倍速い。dbinom もベクトル化されている。使わないという手はない。

p<- seq(0,1,0.00001)
pp <- dbinom(4, 15, p)

コメント (1)
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:最小全域木とエッジの長さの和

2011年10月24日 | ブログラミング

http://d.hatena.ne.jp/ryamada/20110825

どうも,作者は for がお好きなようだが。

Npt<-100
k<-4
X<-matrix(0,Npt,k)
for(i in 2:Npt){
    X[i,]<-X[i-1,]+rnorm(k)
}

は,以下と同じことをしているはずだ。乱数発生もまとめて行う方が効率的。Npt=10000 にして計測すると,この部分だけだと 10 倍速いけど,それは全体の 5% 程でしかないので,トータルでは違いはないということになる。

Npt <- 100
k <- 4
x <- matrix(rnorm(Npt*k), Npt)
x <- matrix(apply(x, 2, cumsum), Npt)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:パーミュテーションテストのシャッフルの準備

2011年10月24日 | ブログラミング

http://d.hatena.ne.jp/ryamada/20110824

説明がよくわからないけど,プログラムを読む限りでは,行列の列を指定して,列ごとの要素をシャッフルするようなプログラムと言うこと。

sample 関数は,"与えられた要素をサンプリング(シャッフル)する" のだ。与える要素は 1:n でもよいし,x[1], ..., x[n] でもよい。しかし,1:n をシャッフルして,それに基づいて,x をシャッフルする必要などない。つまり,x <- c(4, 6, 2, 8) をシャッフルするなら,sample(x) だけでよい。よって,プログラムは簡単至極なものになる。

shuffle <- function(M, s=NULL)
{
    if (is.null(s)) {
        s < 2:ncol(M)
    }
    return(cbind(M[, -s], sapply(s, function(i) sample(M[, i]))))
}

実行例

> set.seed(888)
> n1 <- 3
> n2 <- 4
> (M <- matrix(sample(n1*n2), n1, n2))
     [,1] [,2] [,3] [,4]
[1,]    1    7    3    6
[2,]    4    9    2   11
[3,]   12   10    8    5
> shuffle(M, s=2:4)
     [,1] [,2] [,3] [,4]
[1,]    1    7    8    5
[2,]    4   10    2   11
[3,]   12    9    3    6

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:カーネル密度分布推定の中身を覗く

2011年10月24日 | ブログラミング

http://d.hatena.ne.jp/ryamada/20110820 は,ベクトルプログラミングで 90 倍速くなる

300 ポイント作成では system.time が測定しきれないので,3000 ポイントに増やして測定する。

Npt <- 30000
x <- y <- rep(0, Npt)
r <- 10
t <- runif(Npt)*2*pi/1
r2 <- r*(sin(t*200)+10)+rnorm(Npt)*0.1
x <- r2*cos(t)
y <- r2*sin(t)
x <- cbind(x, y)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:正八面体タイル

2011年10月24日 | ブログラミング

正八角形タイルか?

http://d.hatena.ne.jp/ryamada/20110717 も以下のようになる

速度はほぼ2,3倍。なぜにして,for (i ... が必要なのか線を濃くするだけ?無駄も多いし。for は大敵。

func <- function(ox, oy, or) # 円を描く関数
{
    lines(ox+xc*or, oy+yc*or)
}
depth <- 3
k <- 8
a <- 1
x0<-0
y0<-0
as <- seq(0, 2*pi, length=100)
xc <- cos(as)
yc <- sin(as)
ts <- seq(0, 2*pi, length=8)*(k-1)/k
plot(0, 0, type="n", xlim=c(-1.5, 1.5), ylim=c(-1.5, 1.5))
b <- a/2*tan(pi/k)
c <- sqrt((a/2)^2+b^2)
r <- c*tan(2*pi/k)
R <- sqrt(r^2+c^2)+c
func(x0, y0, a)
sapply(ts, function(theta) func(R*cos(theta), R*sin(theta), r))

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:双曲線幾何の直線その2

2011年10月24日 | ブログラミング

http://d.hatena.ne.jp/ryamada/20110713 のプログラムも読みにくく,効率が悪い。しかも途中で負数の平方根を求めるということになるのが不都合。

下のように書き直せば,100倍ほど速くなる。

func <- function(ox, oy, or) # 円を描く関数
{
     lines(ox+xc*or, oy+yc*or)
}
r1 <- 1
x2 <- 1.5
y2 <- 0
r2 <- 2
t <- seq(0, 2*pi, length=200)
xc <- cos(t)
yc <- sin(t)
plot(0, 0, type="n", xlim=c(-1.98, 3.5), ylim=c(-1.98, 3.5), asp=1)
func(x1, y1, r1)
func(x2, y2, r2)
k <- (x2^2+r1^2-r2^2)/(2*x2*xc)
r3 <- k^2-r1^2
mapply(function(x, y, z) if(z > 0) func(x*cos(y), x*sin(y), sqrt(z)), k, t, r3)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:双曲線幾何の直線

2011年10月24日 | ブログラミング

http://d.hatena.ne.jp/ryamada/20110709 のプログラム

もとのプログラムは,読みにくいこともさることながら,無駄なことをやっている。以下のように簡単になる。

## hyperbolic disk

t <- seq(0, 2*pi, length=100)
xc <- cos(t)
yc <- sin(t)
r <- tan(seq(from=0, to=pi/2, length=10))[-1]
func <- function(ox, oy, or) # 円を描く関数
{
    lines(ox+xc*or, oy+yc*or)
}
plot(0, 0, type="n", xlim=c(-1.5, 3), ylim=c(-1.5, 3), asp=1)
func(0, 0, 1)
sapply(c(r, -r), function(or) func(1, or, or))
abline(h=0)

短いプログラムはたいてい速い。案の定,元のプログラムの8倍ほど速くなっている。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村