学力テストの主成分分析のバイプロット

奥村先生の，学力テストの主成分分析のバイプロット

https://www.scoopnest.com/ja/user/h_okumura/903961406004805632-

に対して，

Taylor さんから，

> すごい。この図を一目見ただけで、視覚的に(←重要)とても多くのことが分かるってのはまさに主成分分析の威力か。あえて言うなら、第１軸、第２軸の寄与率が多少気になるところ。

というコメントがあってですね。

寄与率は
> d = read.csv("atest2017.csv")[,-1]
> rownames(d) = c("北海道", "青森", "岩手", "宮城", "秋田",
+ "山形", "福島", "茨城", "栃木", "群馬",
+ "埼玉", "千葉", "東京", "神奈川", "新潟",
+ "富山", "石川", "福井", "山梨", "長野",
+ "岐阜", "静岡", "愛知", "三重", "滋賀",
+ "京都", "大阪", "兵庫", "奈良", "和歌山",
+ "鳥取", "島根", "岡山", "広島", "山口",
+ "徳島", "香川", "愛媛", "高知", "福岡",
+ "佐賀", "長崎", "熊本", "大分", "宮崎",
+ "鹿児島", "沖縄")
> a = prcomp(d) # scale=FALSE
> print(a)
Standard deviations (1, .., p=8):
[1] 4.8600239 2.4245367 1.4255179 0.9540044 0.6784466 0.6582192 0.4617722 0.3711827

Rotation (n x k) = (8 x 8):
           PC1        PC2        PC3        PC4         PC5        PC6         PC7         PC8
小国A 0.2685353  0.2966870  0.3536034 -0.3922446 -0.47633410  0.5021088 -0.26976989  0.09587030
小国B 0.3559633  0.3486860  0.1517456  0.6129057 -0.36436550 -0.4430136 -0.09111218  0.12531363
小算A 0.2876291  0.5301904 -0.1242260 -0.4727719  0.47286916 -0.3866795 -0.13117396  0.08300222
小算B 0.3636727  0.3189432 -0.5000904  0.2562970  0.09810079  0.5004983  0.35105253 -0.25860909
中国A 0.3181836 -0.1678676  0.3342071 -0.2460678 -0.09212964 -0.1682090  0.80864133  0.08748155
中国B 0.3775471 -0.2166936  0.5331339  0.2060371  0.47944475  0.1294554 -0.22278065 -0.43308407
中数A 0.4343985 -0.4446556 -0.4062962 -0.2579537 -0.34817753 -0.2679320 -0.22925760 -0.36930797
中数B 0.3920535 -0.3670341 -0.1594133  0.1083912  0.21801803  0.1864553 -0.14962358  0.75480735

> x = a$sdev
> x^2/sum(x^2)
[1] 0.701197578 0.174510320 0.060326614 0.027018693 0.013664536 0.012861885 0.006330224 0.004090149
> cumsum(x^2)/sum(x^2)
[1] 0.7011976 0.8757079 0.9360345 0.9630532 0.9767177 0.9895796 0.9959099 1.0000000

ということで，第1，第2主成分の寄与率は，70.1% と 17.5% で，かなり違っているということです。

奥村先生の示したバイプロットは

> biplot(a, cex=0.6) # ,  pc.biplot=FALSE;  default

です。prcomp(d, scale=FALSE) であることに注意。また，biplot(a, cex=0.6, pc.biplot=FALSE) であることにも注意。

biplot(a2, cex=0.6,  pc.biplot=TRUE) とする場合には，目盛りが変わります。

> biplot(a2, cex=0.6,  pc.biplot=TRUE)

なぜ，目盛りが変わるかというと，

pc.biplot

If true, use what Gabriel (1971) refers to as a "principal component biplot", with lambda = 1 and observations scaled up by sqrt(n) and variables scaled down by sqrt(n). Then inner products between variables approximate covariances and distances between observations approximate Mahalanobis distance.

ということなのではありますが，この図からそんな情報を読み取るという意気込みを持つ人はいないと思われ，むしろ，単純に図を見て，「誤った情報をくみ取る」ことがあるのではないかと危惧されます。たとえば，「福井と中数Aが近い！！」とか。

さて，もう一方で，主成分分析は，変数を正規化するかどうか（分散共分散行列からスタートするか，相関係数行列からスタートするか）というのがあり，奥村先生は前者を取ったわけですが，後者を取った場合にはどうなるか，

> b = prcomp(d, scale=TRUE)
> print(b)
Standard deviations (1, .., p=8):
[1] 2.3667078 1.1682792 0.7013992 0.4755625 0.3574532 0.3142276 0.2452455 0.1704118

Rotation (n x k) = (8 x 8):
           PC1        PC2         PC3        PC4          PC5         PC6        PC7         PC8
小国A 0.3452283  0.3287580  0.45731718 -0.3634967 -0.621926560  0.06238154 -0.1916828  0.06356397
小国B 0.3643624  0.2958606  0.04867583  0.6916666  0.004563749 -0.50656879 -0.1764822  0.10562112
小算A 0.3093163  0.5172346 -0.12949306 -0.3864072  0.650100244  0.07318099 -0.1923574  0.07566833
小算B 0.3560964  0.2752739 -0.54310255  0.1375482 -0.279154784  0.33005870  0.4938781 -0.22969101
中国A 0.3772570 -0.2566438  0.34553984 -0.1992751  0.207000000 -0.33832347  0.6852228  0.07655498
中国B 0.3675927 -0.2743884  0.39488637  0.3058082  0.237866718  0.53923815 -0.1760226 -0.40223179
中数A 0.3433546 -0.4017514 -0.37533376 -0.2957626 -0.101814031 -0.39864743 -0.3525097 -0.44587388
中数B 0.3608830 -0.3981434 -0.24607800  0.0230828 -0.053227183  0.24706694 -0.1625342  0.74824164
> biplot(b, cex=0.6, pc.biplot=TRUE)

あまり変わらないと言えば変わらないが。

しかし，いずれにしても，biplot は ? biplot で見るように

There are many variations on biplots (see the references) and perhaps the most widely used one is implemented by biplot.princomp.

ということで，どのような基準（やりかた）で書かれたものかを理解していないと，誤った判断をしてしまう可能性が大である。

ということで，主成分（重み）ではなく，主成分負荷量と主成分得点を別々に描く方がよいのではないかと思う次第ではある。

http://aoki2.si.gunma-u.ac.jp/R/prcomp2.html にあるような，主成分分析の結果を表すときによく使われる（？）普通の表現で，

> prcomp2(d)
                　PC1     PC2     Contribution
小国A            　0.817   0.384     0.815
小国B            　0.862   0.346     0.863
小算A            　0.732   0.604     0.901
小算B            　0.843   0.322     0.814
中国A            　0.893  -0.300     0.887
中国B            　0.870  -0.321     0.860
中数A            　0.813  -0.469     0.881
中数B            　0.854  -0.465     0.946
Eigen.values    　5.601   1.365
Proportion      　70.016  17.061
Cumulative.prop.　70.016  87.077

> x = t(b$sdev*t(b$rotation))
> #quartz("atest2017", width=5, height=5)
> plot(x[,1], x[,2], pch=19, xlim=c(-1,1), ylim=c(-1,1), asp=1, xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
> axis(1, at=-4:4*0.25, pos=0)
> axis(2, at=-4:4*0.25, pos=0)
> text(x[,1], x[,2], colnames(d), pos=c(2,4,4,1, 4,2,2,4))
> abline(h=0, v=0)
> theta = seq(0, 2*pi, length=1000)
> x = cos(theta)
> y = sin(theta)
> lines(x, y, lty=3, col="red")

これによって描かれるのは，主成分負荷量。赤の破線は主成分負荷量の二乗（寄与率）が 1 となる境界。これに近いほど，寄与率が高いということ。

主成分得点のプロットは
> plot(b$x[,1], b$x[,2], pch=19, cex=0.6, xlab="PC1", ylab="PC2")
> text(b$x[,1], b$x[,2], rownames(d), cex=0.6, pos=4)

最後の 2 つを 1 枚のグラフにまとめたのが biplot の一例。

しかし，2 枚の図を別々に見れば，変な誤解は生まれようがない。

まとめ

・prcomp で scale = FALSE / TRUE いずれにするか
・biplot で pc.biplot=FALSE / TRUE いずれにするか
・biplot をやめるかどうか

 

移動量が最小のハンカチ落とし

移動量が最小のハンカチ落とし

締め切りが 2017/09/05 10:00 AM なので，その 1 分後に投稿されるように予約

設問

n人に加えて鬼一人がハンカチ落としをしています。
ハンカチ落としでは、鬼以外の人が円になって座ったあと、鬼が円の外を走ります。
鬼が誰かの後ろでハンカチを落とすと、落とされた人は鬼が一周してくるまでの間に気付き、鬼を追いかけます。

なお、走る速さは全員が同じため、鬼が追い付かれることはないものとします。
また、落とされた人はすぐに気付いて追いかけ、落とした人はその場所まで一周回ってきて座るものとします。

このとき、鬼になった人は他の鬼と違う量だけ移動してハンカチを落とすことにします。
これを繰り返し、すべての位置に一度ずつハンカチを落とすことを考えます。
（すべての位置にハンカチが落ちた時点で終了します。）

最初に鬼はAの位置にハンカチを落とします。
例えば、4人の場合、1人分→2人分→3人分を順に移動して落とすと、すべての位置に一度ずつハンカチが落ちます。
同様に、3人の場合は1人分→4人分を順に移動、もしくは2人分→5人分を順に移動して落とす方法などが考えられます。

このとき、鬼の移動量が最小になるような移動方法を求め、その移動量の和を求めてください。
上記の通り、3人の場合は5(1+4)、4人のときは6(1+2+3)が最小となります。
標準入力から整数 n が与えられたとき、鬼の移動量が最小となる移動方法の移動量の和を標準出力に出力してください。
なお、n は9以下の正の整数が与えられるものとします。

【入出力サンプル】
標準入力
4

標準出力
6

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

f = function(n) {
    library(e1071)
    if (n == 3) {
        p = permutations(4)
        p = unique(p[p[, 1] == 1, 1:2])
    } else {
        p = permutations(n)
        p = cbind(1, p[, 1:(n - 2)] + 1)
    }
    p = p[order(rowSums(p)), ] # 移動量の合計が小さい順に試す
    for (i in 1:nrow(p)) {
        index = 0 # ハンカチが落とされる場所（n を法として）
        position = logical(n) # ハンカチが落とされた場所の記録
        position[index + 1] = TRUE
        ok = TRUE
        for (j in p[i, ]) {
            index = (index + j)%%n
            if (position[index + 1]) {
                ok = FALSE
                break
            } else {
                position[index + 1] = TRUE
            }
        }
        if (ok) {
            cat(sum(p[i,])) # 小さい順に探しているので，見つかったら探索終了
            break
        }
    }
}

# f(scan(file("stdin", "r")))

f(3) # 5
f(4) # 6
f(5) # 12
f(6) # 15
f(7) # 23
f(8) # 28
f(9) # 38

n が偶数のときは sum(1:(n-1)) である。
n が奇数のときは (4-n+n^2) / 2 であろう。