裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

Mountain Lion 用の graphviz

2012年07月31日 | ブログラミング

2012/07/31 現在で,まだ用意されていないようだ

コメント

Mountain Lion でアプリケーションのインストール

2012年07月31日 | ブログラミング

“セキュリティ”環境設定でインストールが許可されているのは、Mac App Store と確認済みの開発元からのアプリケーションのみです。

といわれる。システム環境設定を弄ってもよいが,そこに以下の示唆が表示されているので,それに従うのが無難。

開発元が不明なアプリケーションを開きたい場合は、Control キーを押したままアプリケーションアイコンをクリックして“開く”を選択することにより、個別に許可することもできます。

コメント

X11 on Mountain Lion

2012年07月31日 | ブログラミング

X11 is not included with Mountain Lion, but X11 server and client libraries for OS X Mountain Lion are available from the XQuartz project: http://xquartz.macosforge.org. You should use XQuartz version 2.7.2 or later.

コメント

あまりにも,ひどすぎるんです

2012年07月07日 | ブログラミング

「ダメ出し」シリーズをお送りしておりますが,ほとんどの人は「なんで,そんなこというの?」,「どこに問題があるっていうの?いいんじゃない?」ということだと思います。

そんな事じゃないんです。悪いんです。悪いと思わないことが悪いんです。

とだけ,言っておきましょう。

分からない人は,今後も,ダメなプログラムを書き続けるだけのことでしょう。

嘆かわしい(^_^;)

コメント

ダメ出し:間違えてるよ~

2012年07月07日 | ブログラミング

R によるデータの統計的取り扱い/標準誤差」において

したがって、回帰直線の「切片」の95%信頼区間の上限値・下限値は

resultcoefficients[1] + qt(0.975, length(x)-1) * resultcoefficients[3]resultcoefficients[1] - qt(0.975, length(x)-1) * resultcoefficients[3]

同様に、説明変数 x の回帰係数の95%信頼区間の上限値・下限値は

resultcoefficients[2] + qt(0.975, length(x)-1) * resultcoefficients[4]
resultcoefficients[2] - qt(0.975, length(x)-1) * resultcoefficients[4]

で求めることができます。

うそで~~す。qt の自由度は length(x)-2 でなければなりません。

resultcoefficients は result$coefficients ですし。

confint を使うと,一発で答えが出るヨン。

この標準誤差は何のためにあるのでしょうか。答えは、ズバリ、信頼区間の計算のためです。

これも半分ウソ。それだけのためにある訳じゃあない。

===============================

追記 2012/08/01

コメントで,「length(x)-1 でいいんじゃないの?」ともらったので,そうではないということをはっきりさせよう。

元の記事のプログラムでは,以下のようになる。

> x <- c(61, 72, 84, 95, 97, 98, 100, 113, 126, 130)
> y <- c(83, 82, 99, 96, 115, 108, 95, 111, 114, 135)
> d.fr <- data.frame(x, y)
> RegModel <- lm(y ~ x, data = d.fr)
> summary(RegModel)

Call:
lm(formula = y ~ x, data = d.fr)

Residuals:
    Min      1Q  Median      3Q     Max
-10.362  -5.865   0.100   4.024  11.591

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  40.2680    12.3491   3.261 0.011514 *  
x             0.6509     0.1238   5.259 0.000765 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.104 on 8 degrees of freedom
Multiple R-squared: 0.7756,    Adjusted R-squared: 0.7476
F-statistic: 27.66 on 1 and 8 DF,  p-value: 0.0007653

> result <- summary(RegModel)
> result$coefficients[1]
[1] 40.26801
> result$coefficients[2]
[1] 0.6509425
> result$coefficients[3]
[1] 12.34913
> result$coefficients[4]
[1] 0.1237739

切片の信頼区間が

> result$coefficients[1] + qt(0.975, length(x)-1) * result$coefficients[3]

[1] 68.20369
> result$coefficients[1] - qt(0.975, length(x)-1) * result$coefficients[3]
[1] 12.33233

回帰係数の信頼区間が

> result$coefficients[2] + qt(0.975, length(x)-1) * result$coefficients[4]

[1] 0.9309385
> result$coefficients[2] - qt(0.975, length(x)-1) * result$coefficients[4]
[1] 0.3709466

となるとしているのだが,間違いである。それは,confint を使ってみればわかる。

> confint(RegModel)
                 2.5 %     97.5 %
(Intercept) 11.7908555 68.7451654
x            0.3655195  0.9363656

値が違うのは,qt の自由度が length(x)-1 というのが間違えているから。これを length(x)-2 にしてやってみると,

> result$coefficients[1] + qt(0.975, length(x)-2) * result$coefficients[3]

[1] 68.74517
> result$coefficients[1] - qt(0.975, length(x)-2) * result$coefficients[3]
[1] 11.79086
> result$coefficients[2] + qt(0.975, length(x)-2) * result$coefficients[4]
[1] 0.9363656
> result$coefficients[2] - qt(0.975, length(x)-2) * result$coefficients[4]
[1] 0.3655195

となり一致することがわかる。

confint が間違ってんじゃね?」という意見は,即時却下。



コメント (3)

ダメ出し:R の特徴を生かしたプログラミング

2012年07月07日 | ブログラミング

割り付けがうまく行く確率」にて

k <- 5
p <- runif(k)
p <- p/sum(p)

library(gtools)

cmb <- combinations(k,3)

P <- rep(0,length(cmb[,1]))
for(i in 1:length(cmb[,1])){
    P[i] <- prod(p[cmb[i,]])
}

for を使わないようにというのではなく,コンパクトに書こうということ。
コンパクトになるだけでなく,なにを計算しているかがすぐにわかるというメリットがある。
ついでに,length(cmb[,1]) というのは,nrow(cmb)。これも,関数名から何をしようとしているのかすぐにわかる
以下のようにすれば,前もって P のメモリ確保をする必要もない(メモリ確保も P <- numeric(nrow(cmb)) とすべし)。1 行だけで終わる。

P2 <- apply(cmb, 1, function(x) prod(p[x]))

同じ結果になることも確かめておこう。

> identical(P, P2)
[1] TRUE



コメント

ダメ出し:他人の作った(ライブラリに入っている)関数の使い方はよく調べてから...

2012年07月05日 | ブログラミング

上昇トレンドの検定 ヨンクヒール・タプストラ検定」にて

SAGx というライブラリ(パッケージ)は普通にはないので,代替法を探していて,
http://aoki2.si.gunma-u.ac.jp/R/Jonckheere.html
にあるとか,cor.test で method="kendall" とするのと同じであるとかの収穫があった。

で,ブログ著者は,「この本の例をやってみる。」...「本ではものすごい小さいp値になるらしいけど1になる。保留。」としているが,

JT.test(y[1, ], rep(1, 7), alternative="increasing")

というのは,書式を間違えていると思う。だから,変な結果になっただけだろう。実際の所,以下のように確かに P 値はとっても小さな値になる。

> y <- rbind(c(1.06, 1.03, 1.10, 1.66, 1.72, 1.82, 1.91),
+            c(0.98, 1.02, 1.08, 1.58, 1.68, 1.78, 1.89),
+            c(1.09, 1.10, 1.11, NA, 1.71, 1.83, 1.92))
> y <- c(y)
> g <- rep(1:7, each=3)

> Jonckheere(y, g)

    ヨンキー検定

data:  y by g
J = 164.5000, E(J) = 85.5000, V(J) = 231.5250, Z-value = 5.1919,
p-value = 1.041e-07
alternative hypothesis: control <= treat-1 <= ... <= treat-n

> cor.test(y, g, method="kendall", alternative="greater")

    Kendall's rank correlation tau

data:  y and g
z = 5.1919, p-value = 1.041e-07
alternative hypothesis: true tau is greater than 0
sample estimates:
      tau
0.8788771

コメント

Jonckheere 検定の漸近近似検定は...

2012年07月05日 | ブログラミング

cor.test(x, g, method="kendall") だということ...

http://aoki2.si.gunma-u.ac.jp/R/Jonckheere.html

の例題

> x <- c(
+     153, 153, 152, 156, 158, 151, 151, 150, 148, 157,  # 第 1 群のデータ
+     158, 152, 152, 152, 151, 151, 157, 147, 155, 146,  # 第 2 群のデータ
+     153, 146, 138, 152, 140, 146, 156, 142, 147, 153,  # 第 3 群のデータ
+     137, 139, 141, 141, 143, 133, 147, 144, 151, 156   # 第 4 群のデータ
+     )
> g <- rep(1:4, each=10)
> Jonckheere(x, g, alternative="decreasing")

    ヨンキー検定

data:  x by g
J = 446.5000, E(J) = 300.0000, V(J) = 1705.0776, Z-value = 3.5479,
p-value = 0.0001942
alternative hypothesis: control >= treat-1 >= ... >= treat-n

> cor.test(x, g, method="kendall", alternative="less")

    Kendall's rank correlation tau

data:  x and g
z = -3.5479, p-value = 0.0001942
alternative hypothesis: true tau is less than 0
sample estimates:
       tau
-0.4391269

と同じだということ。知らなかった。

 

コメント

ダメ出し:整数定数,整数演算,演算順序

2012年07月04日 | ブログラミング

R で Excel っぽい色を出す」にて

excel.like.color <- function(n) {
  n <- as.integer(n)
    :
  residue <- (n - 1) %% 6 + 1
  quotient <- floor((n - 1) / 6)

という件があるが,整数定数,整数演算,演算順序を考えると以下の test2 で示すようにするのがよいだろう。

test <- function(n) {
    n <- as.integer(n)
    residue <- (n - 1) %% 6 + 1
    quotient <- floor((n - 1) / 6)
    return(c(i, residue, quotient))
}

test2 <- function(n) {
    n <- as.integer(n)
    quotient <- (n - 1L) %/% 6L
    residue <- n - quotient * 6L
    return(c(n, residue, quotient))
}

> for (i in 1:15) {
+     cat(test(i), "     ", test2(i), "\n")
+ }
1 1 0       1 1 0
2 2 0       2 2 0
3 3 0       3 3 0
4 4 0       4 4 0
5 5 0       5 5 0
6 6 0       6 6 0
7 1 1       7 1 1
8 2 1       8 2 1
9 3 1       9 3 1
10 4 1       10 4 1
11 5 1       11 5 1
12 6 1       12 6 1
13 1 2       13 1 2
14 2 2       14 2 2
15 3 2       15 3 2

コメント

ダメ出し:当たり前のことが当たり前にできるように

2012年07月03日 | ブログラミング

幾何平均」にて

max.g <- 0
min.g <- 0
for(i in 1:n.m){
    max.g <- max.g + max(D[[i]])
    min.g <- min.g + min(D[[i]])
}
max.g
min.g

当たり前に,

max.g <- sum(sapply(D, max))
min.g <- sum(sapply(D, min))

と書くようにしたいものだ。

コメント

ダメ出し:つまらない繰り返しをプログラムしない

2012年07月02日 | ブログラミング

幾何平均」で

AllelesId<-NULL
AllelesId[[1]]<-c("7",9,10,11,12,13,14,15,16,17,18)
AllelesId[[2]]<-c("27",28,28.2,29,30,30.2,30.3,31,31.2,32,32.2,33,33.1,33.2,34,34.2)
AllelesId[[3]]<-c("7",6,9,9.1,10,10.1,10.3,11,12,13,14,15)
AllelesId[[4]]<-c("7",8,9,10,11,12,13,14,15,16)
AllelesId[[5]]<-c("12",13,14,15,16,17,18,19,21)
AllelesId[[6]]<-c("5",6,7,8,9,9.3,10)
AllelesId[[7]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[8]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[9]]<-c("16",17,18,19,20,21,22,23,24,25,26,27,28)
AllelesId[[10]]<-c("9.2",10.2,11,11.2,12,12.2,13,13.2,14,14.2,15,15.2,16,16.2,17.2,18)
AllelesId[[11]]<-c("13",14,15,16,17,18,19,20,21,22)
AllelesId[[12]]<-c("7",8,9,10,11,12,13,14)
AllelesId[[13]]<-c("10",11,12,13,14,15,16,17,17.1,18,19,20,21,22,23,24,25,26,27)
AllelesId[[14]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[15]]<-c("17",18,19,20,21,22,22.2,23,23.2,24,24.2,25,25.2,26,27,28)
などとしているけど,そもそも,15 個の要素を持つ list を作りたいなら,vector("list", 15) とすべきだし,1000 個もあったらどうするつもり? AllelesId[[1]]AllelesId[[1000]] まで書くの?

プログラムというのは,「繰り返し」をそのまま書かないで済むようにすべき。だれも,1から1000までの和をとるときに 1+2+3+…+1000 なんて書かないでしょ?それと同じ。

以下のように書けばよいだけ。

AllelesId2 <- list(c("7",9,10,11,12,13,14,15,16,17,18),
c("27",28,28.2,29,30,30.2,30.3,31,31.2,32,32.2,33,33.1,33.2,34,34.2),
c("7",6,9,9.1,10,10.1,10.3,11,12,13,14,15),
c("7",8,9,10,11,12,13,14,15,16),
c("12",13,14,15,16,17,18,19,21),
c("5",6,7,8,9,9.3,10),
c("7",8,9,10,11,12,13,14,15),
c("7",8,9,10,11,12,13,14,15),
c("16",17,18,19,20,21,22,23,24,25,26,27,28),
c("9.2",10.2,11,11.2,12,12.2,13,13.2,14,14.2,15,15.2,16,16.2,17.2,18),
c("13",14,15,16,17,18,19,20,21,22),
c("7",8,9,10,11,12,13,14),
c("10",11,12,13,14,15,16,17,17.1,18,19,20,21,22,23,24,25,26,27),
c("7",8,9,10,11,12,13,14,15),
c("17",18,19,20,21,22,22.2,23,23.2,24,24.2,25,25.2,26,27,28))

なぜ,第1要素だけ文字列にしているのかな?全部文字列にするために全てに " " を付けるのをきらったからでしょう。

これで同じものができたかどうかは,identical 関数で調べれば分かる。

> identical(AllelesId, AllelesId2)
[1] TRUE

「完全に同じだ」ということがわかる。

ついでだからだけど,なぜ,以下のように書かないのかはわからない。単にテストデータだからか?だったらこそ,簡潔に書きたいと思うのでは?

AllelesId3 <- list(c("7",9,10,11,12,13,14,15,16,17,18),
c("27",28,28.2,29,30,30.2,30.3,31,31.2,32,32.2,33,33.1,33.2,34,34.2),
c("7",6,9,9.1,10,10.1,10.3,11,12,13,14,15),
c("7",8:16),
c("12",13:19,21),
c("5",6:9,9.3,10),
c("7",8:15),
c("7",8:15),
c("16",17:28),
c("9.2",10.2,11,11.2,12,12.2,13,13.2,14,14.2,15,15.2,16,16.2,17.2,18),
c("13",14:22),
c("7",8:14),
c("10",11:17,17.1,18:27),
c("7",8:15),
c("17",18:22,22.2,23,23.2,24,24.2,25,25.2,26:28))

当然だが,これも,最初の書き方でできるものと完全に同じだ。

> identical(AllelesId, AllelesId3)
[1] TRUE

ちなみに,【幾何平均をもって、「代表値」とすることの意義はあるのか、あるとしたらどのあたりにあるのか…はまだ不明】ということの真意はよく分からないけど,対数正規分布の代表値は幾何平均だよね。

コメント