裏 RjpWiki

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

nnet の entropy 引数

2010年07月14日 | ブログラミング
entropy 引数を指定するとエラーが出る。なんでかなぁ。

> nnet(Group~X1+X2, data=df, size=2, skip=FALSE, entropy=TRUE)
以下にエラー nnet.default(x, y, w, entropy = TRUE, ...) :
仮引数 ”entropy” が複数の実引数にマッチしました
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

hclust で使う dist

2010年07月14日 | ブログラミング
hclust(dist(x)^2, method="ward")
のように,ユークリッドの二乗距離でないとおもしろくない
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

tree と mvpart

2010年07月14日 | ブログラミング
デフォルトで走らせると,双方の結果が食い違うことが多い(というか,殆ど食い違う)。
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

tree の plot

2010年07月14日 | ブログラミング
ans <- tree(.....) のあと,
plot(ans) とやっただけでは,二進木しか書いてくれない。
変数名と分類基準を文字で書き加えるには,
text(ans) をやらないといけない。
ちゃんとそう書かれているのに,plot だけやって描かれないので,あれこれ悩んでいた。
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

R2.7.0 で pairs が落ちる

2010年07月14日 | ブログラミング
最初,日本語の変数名を含むデータファイルを指定したせいかと思ったが,plot(iris) でも落ちてしまう。なんてこった。

−−−−−コメント

pairs.default の cex.labels 引数を明示的に指定することでエラーを避けることはできる。 strwidth が Inf を返すようにエンバグされている?のが原因のようだ。
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

R 2.7.0 以降で legend が落ちる

2010年07月14日 | 裏 RjpWiki
いろいろ不具合が残っているようだ

−−−−−コメント

日本語に関する設定を外してください. 2.7.0のQuartzの日本語は旨くありません. 日本語が必要ならIPAフォントを導入してX11を使うとか... なんだかなあ
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

histogramがちょっと変?

2010年07月14日 | ブログラミング
# これは問題なく動く
for (i in 1:4) {
png("hist%i.png", width=400, height=300)
hist(iris[,i])
dev.off()
}
# histogram は lattice ライブラリに入っている
# これでは,できるファイルは長さが0だ。
# 当たり前だが,正当な png ファイルではない。
for (i in 1:4) {
png("histogram%i.png", width=400, height=300)
histogram(iris[,i])
dev.off()
}
# 個別にやるとちゃんとできるーーーなぜだ?
png("histogram-test.png", width=400, height=300)
histogram(iris[,1])
dev.off()

−−−−−コメント

for(i in 1:4) {
png(sprintf("test%i.png", i))
print(histogram(iris[,i]))
dev.off()
}
のように,print 関数を使わないといけないんだそうだ。なぜだ。
結局,
png("test%i.png")
lapply(iris, histogram)
dev.off()
でやった方がそんな変なことを気にしないでもよいし簡単。
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

斜交回転の因子間相関係数行列(3)

2010年07月14日 | ブログラミング
これでよいかな
因子間相関係数行列は factanal の返すオブジェクトには含まれない(単に表示されるだけ)
promax2 の digits 引数は,因子間相関係数を表示する小数点以下の桁数

promax2 <- function (x, m = 4, digits = 3)
{
sortLoadings <- function(z, cor)
{
order2 <- order(colSums(z^2), decreasing=TRUE)
z <- z[, order2, drop = FALSE]
cor <- cor[order2, order2]
neg <- colSums(z) < 0
z[, neg] <- -z[, neg]
cor[, neg] <- -cor[, neg]
cor[neg,] <- -cor[neg,]
return(list(z=z, cor=cor))
}

if (ncol(x) <2) return(x)
dn <- dimnames(x)
xx <- varimax(x)
x <- xx$loadings
Q <- x * abs(x)^(m - 1)
U <- lm.fit(x, Q)$coefficients
d <- diag(solve(t(U) %*% U))
U <- U %*% diag(sqrt(d))
z <- x %*% U
U <- xx$rotmat %*% U
cor <- solve(t(U)%*%U)
ans <- sortLoadings(z, cor)
z <- ans$z
cor <- ans$cor
dimnames(z) <- dn
class(z) <- "loadings"
colnames(cor) <- rownames(cor) <- colnames(x)
cat("因子間相関係数行列\n")
print(round(cor, digits))
list(loadings = z)
}
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

斜交回転の因子間相関係数行列(2)

2010年07月14日 | ブログラミング
sortLoadings がする,余計なこととは
(1) 因子の順序を入れ替える
(2) 因子単位で因子負荷量の符号を入れ替える
これがあるために,J.Foxが書いた処理の追加は注意が必要である
まだ,フィックスの方針のアナウンスもないが,取りあえずの解決法は以下のようにするとよいだろう
まず,proamx 関数をたとえば promax2 として,promax2 を編集する
promax2 <- promax
fix(promax2)
具体的な編集は,以下のように4行を追加する
promax2 <- function (x, m = 4)
{
if (ncol(x) <2)
return(x)
dn <- dimnames(x)
xx <- varimax(x)
x <- xx$loadings
Q <- x * abs(x)^(m - 1)
U <- lm.fit(x, Q)$coefficients
d <- diag(solve(t(U) %*% U))
U <- U %*% diag(sqrt(d))
dimnames(U) <- NULL
z <- x %*% U
U <- xx$rotmat %*% U
dimnames(z) <- dn
class(z) <- "loadings"
# 以下の4行を追加
print(z)
cor <- solve(t(U)%*%U)
colnames(cor) <- rownames(cor) <- paste("Factor", 1:ncol(x), sep="")
print(cor)
list(loadings = z, rotmat = U)
}
promax解を求めるときには,rotation 引数に promax2 を指定する
factanal(データ・フレーム, 因子数, rotation="promax2")
出力結果はたとえば以下のようになる。
promax2 の内部で書き出されるloadings と最終的に書き出される loadings が異なることがあることに注意が必要。J.Fox のコーディングには,loadings の変更が因子間相関係数行列に及ぼす影響を考慮していない点に注目。
つまりは,出力される因子間相関係数行列は promax2 の内部で書き出される状態の loadings の順序,符号に対するものである。最終的に書き出されるloadings の状態での因子間相関係数行列を出力するにはかなりの改造が必要(それを promax2 の内部で行うとしても)。

> factanal(GATB, 5, rotation="promax2")

Loadings:
Factor1 Factor2 Factor3 Factor4 Factor5
x1 0.678
x2 0.145 -0.102 -0.298 0.546
x3 0.102 -0.198 0.188 0.245 0.220
x4 -0.151 -0.342 0.965
x5 0.292 0.661
x6 -0.712 0.141 0.197
x7 0.119 -0.758 0.211 0.124
x8 1.084 0.151
x9 -0.166 -0.966
x10 0.343 -0.308 0.459
x11 0.146 -0.396 -0.503 -0.165

Factor1 Factor2 Factor3 Factor4 Factor5
SS loadings 1.399 1.693 1.186 1.040 1.645
Proportion Var 0.127 0.154 0.108 0.095 0.150
Cumulative Var 0.127 0.281 0.389 0.484 0.633
Factor1 Factor2 Factor3 Factor4 Factor5
Factor1 1.00000000 -0.52956558 -0.21173426 0.04101527 0.70168658
Factor2 -0.52956558 1.00000000 0.09140868 -0.16630606 -0.64960541
Factor3 -0.21173426 0.09140868 1.00000000 0.15334249 0.04571663
Factor4 0.04101527 -0.16630606 0.15334249 1.00000000 0.12556386
Factor5 0.70168658 -0.64960541 0.04571663 0.12556386 1.00000000

Call:
factanal(x = GATB, factors = 5, rotation = "promax2")

Uniquenesses:
x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11
0.558 0.573 0.663 0.232 0.431 0.360 0.458 0.005 0.231 0.355 0.413

Loadings:
Factor1 Factor2 Factor3 Factor4 Factor5
x1 0.678
x2 0.102 0.145 0.298 0.546
x3 0.198 0.220 0.102 -0.188 0.245
x4 0.965 -0.151 0.342
x5 0.661 -0.292
x6 0.712 -0.141 0.197
x7 -0.119 0.124 0.758 0.211
x8 -0.151 1.084
x9 0.966 -0.166
x10 0.459 0.343 -0.308
x11 0.396 0.146 0.503 -0.165

Factor1 Factor2 Factor3 Factor4 Factor5
SS loadings 1.693 1.645 1.399 1.186 1.040
Proportion Var 0.154 0.150 0.127 0.108 0.095
Cumulative Var 0.154 0.303 0.431 0.538 0.633

Test of the hypothesis that 5 factors are sufficient.
The chi square statistic is 5.49 on 10 degrees of freedom.
The p-value is 0.856
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

あるはずの関数がないと言われるとき

2010年07月14日 | ブログラミング
stats にある factanal 関数を調べるために,関数本体はそのままにして別の名前を付けて実行しようとすると
> print(a(gatb,4), cutoff=0)
エラー: 関数 "factanal.fit.mle" を見つけることができませんでした
と言われる。factanal.fit.mle は stats の中にあるし,第一,factanal を実行しても,そんなエラーメッセージは出ない。
その回避策は,関数中で factanal.fit.mle を使っているところを,stats:::factanal.fit.mle とすればよい。つまり,「どこそこにある関数だよ」と教えるわけだ。
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

斜交回転の因子間相関係数行列

2010年07月14日 | ブログラミング
手こずった。
答えを出したのだけど,それが正しいかどうか,なかなか確認がとれなかった(参考にした本に誤りがあったりなんだかんだ)。
最終的には
http://tolstoy.newcastle.edu.au/R/devel/05/06/1414.html
にあったが,今の print.factanal はそこに呈示されている変更が加わっていないようだ。
ans <- factanal(df, 4) # デフォルトはバリマックス回転
ans2 <- promax(ans$loadings) # それをもとにしてプロマックス回転
a <- ans2$rotmat # 返される回転行列を用いて
solve(t(a)%*%a) # これが因子間相関係数行列

−−−−−コメント

へえ。 http://bugs.r-project.org/cgi-bin/R/wishlist?id=10931;expression=factanal;user=guest によれば,Mon, 10 Mar 2008 13:59:05 +0100 (CET) に,投稿があったようだ。偶然というか(今になって,これを見た) ----- COMMENT: AUTHOR: アール EMAIL: IP: 133.8.33.217 URL: DATE: 03/28/2008 06:44:44 PM 因子負荷量の符号を付け替えるsortLoadings が悪さをするので,上のやり方では相関係数の符号が逆になる組み合わせが出てくる。やはり,J.Fox のコメントを有効にするしかない。しかし,print.factanal を fix では,変更できないような。。。
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

久保さんの集計プログラム

2010年07月14日 | ブログラミング
久保さんの2008/03/11の集計プログラム
以下のようにしてみた(簡単明瞭?)
> set.seed(12345)
> d <- data.frame(n=sample(10, 500, replace=TRUE),
+ year1= sample(LETTERS[1:5], 500, replace=TRUE),
+ year2= sample(LETTERS[1:5], 500, replace=TRUE),
+ year3= sample(LETTERS[1:5], 500, replace=TRUE))
>
> calc <- function()
+ {
+ calc0 <- function(m)
+ {
+ x <- matrix(0, 5, 5) # year.i から year.j への推移行列
+ colnames(x) <- rownames(x) <- LETTERS[1:5]
+ x[rownames(m), colnames(m)]<- m[rownames(m), colnames(m)] # これがすごい
+ return(as.vector(x)) # ベクトルにして返す
+ }
+ ans <- xtabs(n~year1+year2+year3, d) # 3重クロス集計する
+ a1 <- apply(ans, 1:2, sum) # year1 -> year2 の推移行列
+ a2 <- apply(ans, 2:3, sum) # year2 -> year3 の推移行列
+ year1 <- calc0(a1) # 結果のyear1列
+ year2 <- calc0(a2) # 結果のyear2列
+ type.T1 <- rep(LETTERS[1:5], 5) # その他の列の定義
+ type.T2 <- rep(LETTERS[1:5], each=5)
+ type.T1T2 <- paste(type.T1, type.T2, sep="->")
+ result <- data.frame(type.T1=type.T1, type.T2=type.T2,
+ type.T1T2=type.T1T2,
+ year1=year1, year2=calc0(a2)) # データフレームを構成
+ return(result) # 結果を返す
+ }
> calc()
type.T1 type.T2 type.T1T2 year1 year2
1 A A A->A 101 89
2 B A B->A 146 81
3 C A C->A 52 110
4 D A D->A 169 71
5 E A E->A 104 92
6 A B A->B 93 142
7 B B B->B 57 139
8 C B C->B 98 130
9 D B D->B 81 71
10 E B E->B 138 108
11 A C A->C 147 76
12 B C B->C 122 43
13 C C C->C 123 176
14 D C D->C 130 134
15 E C E->C 99 225
16 A D A->D 84 126
17 B D B->D 115 134
18 C D C->D 102 89
19 D D D->D 113 92
20 E D E->D 126 111
21 A E A->E 106 139
22 B E B->E 127 70
23 C E C->E 189 116
24 D E D->E 187 172
25 E E E->E 88 161

ちなみに,久保さんのプログラムと集計結果
> na.to.zero <- function(v) sapply(v, function(x) ifelse(is.na(x), 0, x))
>
> count.type <- function(data, col.count, col.years)
+ {
+ n <- length(col.years)
+ # prepare df.count
+ v.type <- sort(unique(unlist(data[, col.years]))) # 1, 2, 3, ...
+ df.count <- data.frame(type = v.type)
+ matrix.zero <- matrix(0, nrow(df.count), n)
+ df.count <- cbind(df.count, matrix.zero)
+ # count
+ colnames(df.count)[2:(1 + n)] <- col.years
+ for (y in col.years) {
+ v.count <- tapply(
+ data[, col.count],
+ factor(data[, y], levels = as.character(v.type)),
+ sum
+ )
+ df.count[, y] <- na.to.zero(v.count[df.count$type])
+ }
+ return(df.count)
+ }
>
> count.change <- function(data, col.count, col.years, sep = "->")
+ {
+ n <- length(col.years)
+ # prepare df.count
+ v.type <- sort(unique(unlist(data[, col.years]))) # 1, 2, 3, ...
+ df.count <- expand.grid(v.type, v.type) # all possilbe combinations
+ colnames(df.count) <- c("type.T1", "type.T2")
+ df.count$type.T1T2 <- apply(
+ df.count, 1, function(x) paste(x[1], x[2], sep = sep)
+ )
+ matrix.zero <- matrix(0, nrow(df.count), n - 1)
+ df.count <- cbind(df.count, matrix.zero)
+ colnames(df.count)[4:(2 + n)] <- col.years[1:(n - 1)]
+ # count
+ for (y in 1:(n - 1)) {
+ y1 <- col.years[y]
+ y2 <- col.years[y + 1]
+ v.change <- tapply(
+ data[, col.count],
+ factor(
+ apply(data, 1, function(x) paste(x[y1], x[y2], sep = sep)),
+ levels = as.character(df.count$type.T1T2)
+ ),
+ sum
+ )
+ df.count[, y1] <- na.to.zero(v.change[df.count$type.T1T2])
+ }
+ return(df.count)
+ }
> count.change(d, "n", c("year1", "year2", "year3"))
type.T1 type.T2 type.T1T2 year1 year2
1 A A A->A 101 89
2 B A B->A 146 81
3 C A C->A 52 110
4 D A D->A 169 71
5 E A E->A 104 92
6 A B A->B 93 142
7 B B B->B 57 139
8 C B C->B 98 130
9 D B D->B 81 71
10 E B E->B 138 108
11 A C A->C 147 76
12 B C B->C 122 43
13 C C C->C 123 176
14 D C D->C 130 134
15 E C E->C 99 225
16 A D A->D 84 126
17 B D B->D 115 134
18 C D C->D 102 89
19 D D D->D 113 92
20 E D E->D 126 111
21 A E A->E 106 139
22 B E B->E 127 70
23 C E C->E 189 116
24 D E D->E 187 172
25 E E E->E 88 161
同じになっているはず
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

エラーメッセージ

2010年07月14日 | 裏 RjpWiki
エラーメッセージが日本語になっても,意味がわからない,意味をわかろうとしないユーザがたんといるんだな。
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

R2.7の新機能(2)

2010年07月14日 | ブログラミング
The default graphics device in non-interactive use is now pdf() rather than postscript(). [PDF viewers are now more widely available than PostScript viewers.] ということ。未だに eps なんて古くさい。
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ロジスティック回帰

2010年07月14日 | ブログラミング
一変量ロジスティック回帰を,glm でやるときと nls でやるとき,答えが違うのはなぜなのか?
たとえば以下のようなデータ。
y x
1 0 0
2 0 0.1
3 0 0.3
4 0 0.4
5 0 0.4
6 0 0.4
7 0 0.4
8 0 0.4
9 0 0.4
10 0 0.4
11 0 0.7
12 0 1
13 0 1.2
14 0 1.2
15 0 1.4
16 0 1.4
17 0 1.4
18 0 1.4
19 0 1.4
20 0 1.5
21 0 1.6
22 0 1.6
23 0 1.7
24 1 1
25 1 1.4
26 1 1.7
27 1 1.8
28 1 1.9
29 1 2.5
30 1 2.7

−−−−−コメント
simplex 法や nls では残差平方和を最小にする解を求めるようになっているのに対して, glm は Deviance を最小にする解を求めるものであるため(普通は,後者を使うべし)
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

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

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