年齢の計算は暗算でもできる

中澤さんが，

> RのISOdate()関数は1899年以前の日付も扱えるので，例えば明治32年7月1日生まれ，
> 平成13年3月12日に亡くなった方の生存期間は，
> ISOdate(2001,3,12)-ISOdate(1899,7,1)
> とすれば，正しくTime difference of 37144 daysと返ってくる。この結果を365.24で
> 割ると，101歳で亡くなった方ということがわかる。

と，書いている。年齢を求めるのが目的ではなく，生存期間（日数）を求めるのが主目的。

一方，年齢計算は単純で，「生年 y1 と当年 y2 の差を取り，もし当月 m2 当日 D2 が誕生日 M1, D1 の前ならば 1 歳差し引く」

R で書く方が面倒なんだけど，やっていることは簡単で，閏年がどうのとか，365.24 などというマジックナンバーも無関係

age = function(y1, m1, d1, y2, m2, d2) {
    y2 - y1 - (m2 < m1 || (m2 == m1 && d2 < d1))
}

> (ISOdate(2001, 3, 12)-ISOdate(1899, 7, 1)) / 365.24
Time difference of 101.6975 days
> age(1899, 7, 1, 2001, 3, 12)
[1] 101

> (ISOdate(2001, 6, 30)-ISOdate(1899, 7, 1)) / 365.24
Time difference of 101.9987 days
> age(1899, 7, 1, 2001, 6, 30)
[1] 101

> (ISOdate(2001,7, 1)-ISOdate(1899, 7, 1)) / 365.24
Time difference of 102.0014 days
> age(1899, 7, 1, 2001, 7, 1)
[1] 102

ハートを描くプログラム

中澤さんがハートを描くプログラムを書いていた。
----- 引用開始
簡単なのでRでコーディングしてみた。

fg = function(x,y) sqrt(abs(x))+y*sqrt(5-x^2)
s = sqrt(5-1E-15)
curve(fg(x, 1), -s, s, ylim=c(-2,3), col=6)
curve(fg(x, -1), -s, s, col=6, add=TRUE)

で境界線のみの図ができるし，これに続けて，

vfg = Vectorize(fg)
x1 = seq(s, -s, length=100)
x2 = rev(x1)
plot(c(-3, 3), c(-3, 4), type="n")
polygon(c(x1, x2), c(vfg(x1, 1), vfg(x2, -1)), border=6, col=6)

とすれば，これら2つの関数で囲まれた領域を塗りつぶした図ができる。
----- 引用終了
0 を含むようにするなら，length=101 のように，奇数にする方がよい（curve の n のデフォルトは 101 になっている）。そうしないと，ハートのとんがりとへこみがシャープにならない。

Vectorize という関数は，知らなかった。しかし，今回のプログラムでは fg はベクトライズされている。

fg = function(x, y) sqrt(abs(x)) + y * sqrt(ifelse(5 < x^2, 0, 5 - x^2))
n = 101
x = seq(sqrt(5), -sqrt(5), length = n)
x = c(x, rev(x))
s = rep(c(-1, 1), each = n)
plot(0, 0, type = "n", xlim = c(-3, 3), ylim = c(-3, 4), asp = 1)
polygon(x, fg(x, s), border = 2, col = 2) # 概形だけ描くなら col = 2 を除く

分析プログラムにより結果が違う，そんな馬鹿な

Wikipedia の「主成分分析」に関するページ
http://ja.wikipedia.org/wiki/主成分分析　において。

主成分分析を実行するためのソフトウェアや関数によって、観測値の基準化の方法や数値計算のアルゴリズムに微細な差異が多く存在し、必ずしも全く同じ値が出るとは限らない（例えば、R における prcomp関数とFactoMineRのPCA関数の結果は異なる）。

とある。「全く同じ値」ということで，どの程度の違いを指しているのかはっきりしないが，これは，著者の誤解の可能性が大きい。

少なくとも，これを読んだ人が誤解してしまうことはあるようだ。別物を比較して違うといってもしようがない。以下に説明しよう。

--------------------------------------------------------------
1. FactoMineR の PCA は，デフォルトで変数を標準化する。
     scale.unit
       a boolean, if TRUE (value set by default) then data
       are scaled to unit variance

   prcomp では，デフォルトでは変数を標準化しない。普通は，標準化すべきと述べられている。
     scale.    （最後のドットは引数名に含まれる）
       a logical value indicating whether the variables should
       be scaled to have unit variance before the analysis
       takes place. The default is FALSE for consistency with S,
       but in general scaling is advisable. Alternatively, a
       vector of length equal the number of columns of x can be
       supplied. The value is passed to scale.

2. FactoMineR の PCA が返すのは，負荷量（主成分と元の変数の相関係数にあたるもの）である。以下の 2 つであるが，同じものである。(- dimensions が何かは知らない）
     $var$coord       coord. for the variables           
     $var$cor         correlations variables - dimensions

> PCA(iris[1:4])$var$cor
                  Dim.1      Dim.2       Dim.3       Dim.4
Sepal.Length  0.8901688 0.36082989 -0.27565767 -0.03760602
Sepal.Width  -0.4601427 0.88271627  0.09361987  0.01777631
Petal.Length  0.9915552 0.02341519  0.05444699  0.11534978
Petal.Width   0.9649790 0.06399985  0.24298265 -0.07535950

   prcomp が返すものを print メソッドで書いて示されるものは，Rotation と名前が付けられている（要素名は $rotation）が，軸の回転はすなわち変数に対する重みということである。scale.=TRUE とした場合の結果は，もとの変数の相関係数行列の固有ベクトルそのものである。

> p = prcomp(iris[1:4], scale.=TRUE)
> p$rotation
                    PC1         PC2        PC3        PC4
Sepal.Length  0.5210659 -0.37741762  0.7195664  0.2612863
Sepal.Width  -0.2693474 -0.92329566 -0.2443818 -0.1235096
Petal.Length  0.5804131 -0.02449161 -0.1421264 -0.8014492
Petal.Width   0.5648565 -0.06694199 -0.6342727  0.5235971

> e = eigen(cor(iris[1:4]))
> e$vectors
           [,1]        [,2]       [,3]       [,4]
[1,]  0.5210659 -0.37741762  0.7195664  0.2612863
[2,] -0.2693474 -0.92329566 -0.2443818 -0.1235096
[3,]  0.5804131 -0.02449161 -0.1421264 -0.8014492
[4,]  0.5648565 -0.06694199 -0.6342727  0.5235971

prcomp の $rotation から負荷量を計算するのは以下のようにする。$sdev は，固有値の平方根（固有値は主成分の分散であり，その平方根をとると標準偏差 standardized deviation ということである）

> t(p$sdev*t(p$rotation))
                    PC1         PC2         PC3         PC4
Sepal.Length  0.8901688 -0.36082989  0.27565767  0.03760602
Sepal.Width  -0.4601427 -0.88271627 -0.09361987 -0.01777631
Petal.Length  0.9915552 -0.02341519 -0.05444699 -0.11534978
Petal.Width   0.9649790 -0.06399985 -0.24298265  0.07535950

ということで PCA の $var$cor と同じものが得られる。

がしかし，符号が違うじゃないかと思う人も多いようなのだが，主成分単位で絶対値は同じだが符号が異なるのは，符号を付け替えればよいだけなので，「結果が違う」わけではない。固有値の符号はコンピュータプログラムの違い（アルゴリズム，プログラムのコンパイラ，コンピュータアーキテクチャなど）で変わりうる。

3. 主成分得点の分散は，固有値に一致する（平均値は 0）。
   しかし，PCA の主成分得点（$ind$coordは，分散（変動をサンプルサイズで割ったもの。すなわち不偏分散ではない）が固有値に一致する。
   prcomp の主成分得点は，不偏分散（変動を「サンプルサイズ-1」で割ったもの）が固有値に一致する。

> colSums(PCA(iris[1:4])$ind$coord^2)/150
     Dim.1      Dim.2      Dim.3      Dim.4
2.91849782 0.91403047 0.14675688 0.02071484

> colSums(prcomp(iris[1:4], scale.=TRUE)$x^2)/149
       PC1        PC2        PC3        PC4
2.91849782 0.91403047 0.14675688 0.02071484

ということなので，主成分得点の値は両者で若干異なる（が，変換可能である）。

はじめダメなら，最後までダメ

駆け出し統計家のプログラミング手帳
http://mizuki999aiu.blog.fc2.com/
の
χ2 (カイ二乗)検定
http://mizuki999aiu.blog.fc2.com/blog-entry-30.html

帰無仮説と対立仮説の設定が逆になっており，結論が逆。
しっかり！

>> カイ２乗検定はt検定や分散分析などとは仮説の設定が異なるので少し注意です。
>> この場合、
>> 帰無仮説＝2群は関連している。

>> data: travel
>> X-squared = 1.9863, df = 1, p-value = 0.1587

>> p-value>0.05なので帰無仮説を採択し、上記の２群は関連があるということがわかりました。

chisq.test と prop.test の関連性（等価な検定）から，以下と比較しても，筆者の誤解は明らか。

> x <- matrix(c(167, 185, 133, 115), 2)
> addmargins(x)
     [,1] [,2] [,3]
[1,]  167  133  300
[2,]  185  115  300
[3,]  352  248  600
> prop.test(c(167, 185), c(300, 300))

    2-sample test for equality of proportions with continuity correction

data:  c(167, 185) out of c(300, 300)
X-squared = 1.9863, df = 1, p-value = 0.1587
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.14199098  0.02199098
sample estimates:
   prop 1    prop 2
0.5566667 0.6166667

よくある質問

質問！ITmedia > 学問・教育 > 数学

４群の成績の比較。統計　検定、教えてください！
http://qa.itmedia.co.jp/qa8339547.html

この手の質問は結構ある。データもたくさん（240問への回答の正否）取ったので何とかしたいと思うのだろう。とりあえず，データを整理すると，

　　　　科目A 科目B 科目C 科目D
正解　　55    48   54    25  
不正解　5     12   6     35

のようになっていて，科目ごとに正解率が違うか検定したいと。

この表だけ見せられると，おっちょこちょいは「カイ二乗検定でイインジャネ？」とかいうかもしれない。

> chisq.test(matrix(c(55,5, 48,12, 54,6, 25,35), 2))

    Pearson's Chi-squared test

data:  matrix(c(55, 5, 48, 12, 54, 6, 25, 35), 2)
X-squared = 53.5657, df = 3, p-value = 1.389e-11

「オオ，ユウイダネ～。ヨカッタヨカッタ」

しかし，残念なことに「ある一人に問題を解いていただき」とある。その人についてだけ何かがいえても，何の知見が加わるというのだろうか。

統計センスがない。