母集団が正規分布に従っているかを判断する手段は？

正規分布とは

正規分布とはデータを処理する上で非常に都合がいい分布です。

なので、できれば、標本のデータが正規分布に従っていると仮定できれば、統計処理する上で非常に都合がいい訳です。

自然界は正規分布に溢れていて、例えば国民の身長や雨粒の大きさなどが有名でいい例です。

図書館に行くと、統計学に関して様々な議論をしている教科書が山ほどありますが、現場で実際に統計を行う場面になると、

「どのように統計処理すればいいんだ？」

ってなると思います。

その理由の一つとして、このデータは正規分布に従うのか自分で判断する必要があるが、判断方法について教わってないからではないでしょうか。

正規分布を前提としてベラベラ語っている参考書は散見されるけど、正規分布として扱っていいのかについて判断する方法について考えている本ってあんまないですよね。

考えましょう。

分位点を用いた方法

Q-Qplot

qqplotのqはquantileのqで、分位点を示してます。

今回は、1871年から1970年間のナイル川の年間流量のデータを例で使います。

時系列でplotすると

これをヒストグラムにすると

累積密度関数にすると

となった。

この累積密度関数を拡大すると

こうなる。

標準正規分布でこれに対応する値を出してプロットしたものがqqplot.

つまり、

ナイル川のデータをx軸にすると、

(456,α1),(649,α2),......,(1370,α100)をplotしたものがqqplot.

じゃあそのα1〜α100はどう求めるかというと、

標準正規分布の累積密度関数で、それに相当する分位点の値。

例えば分位点の設定を

i/n にすると、

(456,z0.99),(649,z0.98).....(1370,z0)

となるが、z0は無限に発散するので、この分位点は利用できない。

(i-1)/n　も同様に、z100が-無限に発散するので利用できない。

なので、

ナイル川のi番目に相当する分位点αiは, (i-1)/n <α<i/n の範囲で任意に設定すればいい。

よく使われるのがtype5のやつ。　この場合は上記範囲の真ん中を取っている。

どうやってqqplotを再現するかというと、

正規分布の累積密度関数で分位点が0.005, 0.015, ....0.995　に対応する統計量Zをx

ナイル川のデータを昇順に並べたものをy

に収納して

これをプロットしたものがqqplotです。

_

いちいちこんな面倒なことせずに、qqnorm関数を使えば勝手にやってくれます。

分位点でいうtype5でplotしてやると

type8の分位点を採用すると

ってなります。

これが直線になっていれば正規分布と近似できるのではないか、と考える。

なんかモヤモヤするんですよね。

まず、x軸とy軸でデータ値の次元が違ったり(1メモリの値がお互い異なる)、範囲が違う(原点が0,01じゃない)のがモヤモヤ

ただ判断基準は1次関数になっているかどうかであり、傾きとか切片はどうでも良いから大丈夫って感じかな

あと、標本の経験分布関数に対する正規分布はどうせなら標準正規分布でなくて、

標本の平均と分散を同じにした正規分布の方がいいんじゃないかと思うんですよね。

plotすると

素人的にはこっちの方が見やすいと思うんですけど、、

拡大すると

この方法のデメリットは、主観的な判断になってしまうこと。

このナイル川のやつだと、(456,   )が外れ値で邪魔しているが、どう考えるかは検定者に委ねられる。

Shapiro-Francia検定

そこで出てくるのがこの検定方法。じゃじゃじゃじゃーん。

正規分布下でそれに相当する分位点での統計量の期待値(上でいうx)と

Nile側のデータ(上でいうy)

の決定係数(相関係数の2乗)を検定統計量W'とし,

そのW'の分布を用いて検定する。

RだとShapiroFranciaTestとかsf.testで検定できるらしいのですが

なぜか僕のRstudioだとできず、、

そもそもこの検定法の情報が少なくて、W'がどんな分布を呈するのかよくわからず、、

とりあえず決定係数は0.97 でした。

一言で言えば、qqplotしたものが本当に線型化しているのかを確認するために

単回帰分析して、決定係数で評価してるってことですね。

 

ちなみに、

対する正規分布を標準正規分布ではなくて、

平均と分散を標本分布に調整した正規分布で見てみると

決定係数は同じになりますね。

 

原理としてはこんなとこ。

LaTexが苦手でして、許して。

Shapiro-Wilk検定

shapirofranciaを進化させて検定。

xとyの値のばらつきが違うので、分散共分散行列で重みをつけた統計量を用います。

(むしろ、shapirofranciaはshapiro wilkの簡易版ってイメージらしい)

正規分布じゃないらしいです。

 

適合度を用いた方法

kolmogorov検定、Kuiper検定

母集団から導かれる実際の累積密度関数と、順位統計量から導かれる理想的な経験分布関数の差の大きさを用いて統計する方法です。

青い線は実際のデータの累積密度関数のグラフ。

赤い線は帰無仮説で仮定している分布の累積密度関数を示したグラフ。

今回は正規分布かを検討したい訳なので、

今回の場合だと、

平均値と分散がナイル川のデータと同じ正規分布

になります。

確率変数X1,X2....X100の点で、赤の関数F(x)と青の関数Fn(x)の値の差をみる。

青の関数Fn(x)の方が上にある時の差の大きさ　その内の最大値がDn +

赤の関数F(x)の方が上にある時の差の大きさ　その内の最大値がDn -

今回の例だと青の関数が上にいっているのが殆どのようだが、、

Dn +とDn-の最大値　max(Dn+,Dn-) がKolmogorov統計量

Dn+とDn- の和　(Dn+) ＋ (Dn-)　がKuiper統計量

になる。

このKolmogorov統計量を用いた正規性の検定として

Kolmogorov–Smirnov test　があり、

ナイル川のデータを実際にやってみると、

標準正規分布は否定されました。

ただ、平均値と分散が同じ正規分布で検定すると、否定できないという結果になりました。

分布の特性を用いた方法

正規分布の歪度と尖度はそれぞれ0,3なので、

知りたい標本分布の歪度と尖度を求め、それぞれが0,3から離れているかを検定する。

この分野の数学的解釈がまあ難しく、、、

数学的側面は断念して、もう道具として使います、、

中心極限定理による標本平均の分布が正規分布として推定できるサンプルサイズについての検討

中心極限定理とは

　ある程度統計を勉強していると中心極限定理という言葉を耳にします。

定義としては、

　「標本を抽出する母集団の分布が正規布分布に従っても従わなくても、サンプルサイズnが十分に大きくなれば、標本平均の分布は正規分布に近づいていく」

じゃあ正規分布に従うとみなせる程のサンプルサイズとは、ってなるんでしょうか。

正規分布に近似できるなら近似したいですよね。その方が結果を解釈しやすいし。

巷だとn≧20~30程度が目安と言われていますが、エビデンスってないと思うんですよ。

Rを使って実践的に調べてみましょう。

中心極限定理の原理

実際に実験して適切なサンプルサイズを考える

サンプルサイズについて実験する前に、中心極限定理という都合のいい定理がなぜ生じるのか考えます。

まずは正規分布N(0,1)から10個のサンプルを抽出し、その平均の分布について

ヒストグラムとして可視化するために、1000回繰り返しています。

1000回繰り返せば十分だと思っています、、

コード

するとヒストグラムは下図になりました。

当然正規分布っぽくなります。

手抜きで両軸のグラフ名を整えてないけど許して

データが正規分布しているかどうかの検定で、今回はshapiro-wilk test を使用しました。

余談ですがこの検定は

H0(帰無仮説)：データが正規分布に従う

H1(対立仮説)：データが正規分布に従わない

で検定しています。

巷だと、この検定でパラかノンパラどっちの手法が使えるか分かる

（p>0.05だと正規分布を示すからパラメトリック検定、p<0.05だと正規分布を示さないからノンパラメトリック検定）

と言っている人が多いのですが、

確かにp<0.05でH0が棄却されたら正規分布を示さないとみなせますが、

p>0.05の場合は正規分布に従わないとは言えない、って解釈だから正規分布に従っているとは言えないとおもうんですよね。

余談が長いと嫌いな上司みたいになっちゃうのでここまでにします。

二項分布で、probability=0.5,試行回数10回の確率変数を求める行為を10回(サンプルサイズ)繰り返し、その平均について調べます。

ヒストグラムにするために平均値を求める行為を1000回してます。

コード

ヒストグラム

正規分布とは言えないらしいです。

QQplotを見ると

正規分布っぽいけどね。

じゃあサンプルサイズを20にすると

正規分布にならないらしいです。

ヒストグラムやQQplotだと綺麗だけどね

サンプルサイズを30にすると

帰無仮説が棄却されませんでした。

次に確率を変えてみましょう。

probabilityを0.01に鬼下げしました。

サンプルサイズを10にすると

正規分布の面影すらないですね。

サンプルサイズを100にしても

みてられないです。

probabilityを0.001まで下げちゃうと、

まあそうなるよね。

サンプルサイズを1000まで増やすと

ましにはなるけど正規分布にはならないみたいです。

サンプルサイズを10000にしても

正規分布にはならない。

以上より、標本平均の分布はサンプルサイズが大きくなれば中心極限定理により必ず正規分布として扱えるが、必要なサンプルサイズの量は元の分布に依存することがわかります。

必要なサンプル量は？

上の実験より、正規分布として推定できるサンプルサイズは母集団の分布に依存することがわかります。

とはいえ母分散なんて知る由がない上に、正規分布として推定できるかを検定する客観的な方法はなく、ヒストグラムやQQplotから自分で判断するという主観的な物に頼らざるを得ません。

医療統計で用いる母集団ってそんなに汚い分布じゃないと思うんですよね。

わりかし左右対称の分布だけど個人差が大きいから標準偏差が大きい、みたいな。

便宜上だけど一様分布で試しました。

runif(10,0,100)

runif(20,0,100)

runif(30,0,100)

より正規分布との比較ができるようにbreaksを多めにして表示しました。

つまり、

サンプルサイズを多くして大まかな分布を把握し、母分布をR等でモデル化する。

そのモデルを元に得た標本平均のヒストグラムやQQplotの結果を主観的に読み取って、適切な

サンプルサイズを検討する。

となるんでしょうか。

ただ、こうせざるを得ない感が半端なくて　こんなの非現実的というか無理だと思います。

上で示し続けたグラフから分かるように、よほど母分布が変でない限りはサンプルサイズが30程度あればなんとなく正規分布っぽくなるんで、それも利用してt検定とか使っているんでしょうね。

やはり論文を読み解く際はp値だけを見るのではなく、こうした背景にも目を向ける必要があります。

だからこそ論文を読む際は統計学の理解を深めていないと適切な解釈が得られませんねって話です。

2018年NPBプロ野球の主要選手をクラスター分析してみた。

クラスター分析とは

　クラスター分析では複数の個体をグループ各に分類する方法です。

クラスター分析の種類・方法

　クラスター分析では複数の個体をグループに分けていきます。

当然ながら色々なグループの作り方がありそうですよね。

Rでheip(hclust)をrunしてやると

This should be (an unambiguous abbreviation of) one of "ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), "mcquitty" (= WPGMA), "median" (= WPGMC) or "centroid" (= UPGMC).

出典　R

だそうです。

クラスタ間の距離の選び方

• "ward.D" :   D(P,Q) = E(P ∪ Q) - E(P) - E(Q)

                              E(C)=Σd(C,gc):各点と重心の距離の2乗和の根(ユークリッド距離)

• "ward.D2":  "ward.D"の最新版らしい。
• "single" : 最近隣法。2つのクラスターの中で最も個体間の距離が短い値を代表値とする。
• "complete" : 最遠隣法。2つのクラスター中で最も個体間の距離が長い値を代表値とする。
• "average" : 群平均法。2つのクラスター間の各個体の全ての組み合わせの距離の平均を代表値とする。
• "mcquitty" : 以後精進します。
• "median" : 2つのクラスター間の各個体の全ての組み合わせで求まった距離の中央値を代表値とする。
• "centroid" : 重心法。2つのクラスターの重心間の距離を代表値とする。

距離の計算方法

　help(dist)を打ってやると

This must be one of "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski".

と言われます。

• "euclidean" : ユークリッド距離

　　                     n次元の点p,q間の距離　:　sqrt(sum(p-q)^2)　

　                            

• "maximum" : n次元の点p,qの各次元間の距離の絶対値の最大値  max(abs(p-q))
• "manhattan" : n次元の点p,ｑの各次元間の距離の絶対値の和　sum(abs(p-q))
• "canberra" :  sum(abs(p-q) / (abs(p)+abs(q)))  

                          

• "binary" : 2進法なんだろうけどいつ使うのだろうか
• "minkowski" : アインシュタインの特殊相対性理論の説明に使われているやつ。頭がよくなったら振り返ろう。

2018年NPBの主要選手をクラスター分析する

概要

選手が多すぎてもよくないと思うので規定打席(試合数×3.1)・規定投球回(試合数×1.0)を超えた選手に絞ってクラスターしていきます。

セ・リーグ打者のクラスター分析

順位	選　手	チーム名	打率	試合	打席	打数	得点	安打	2塁打	3塁打	本塁打	塁打	打点	盗塁	盗塁刺	犠打	犠飛	四球	敬遠	死球	三振	併殺打	長打率	出塁率
1	ビシエド	(中)	0.348	135	582	512	91	178	26	1	26	284	99	3	4	0	4	51	4	15	61	24	0.555	0.419
2	坂本　勇人	(巨)	0.345	109	502	441	87	152	27	2	18	237	67	9	5	0	0	61	10	0	83	4	0.537	0.424
3	平田　良介	(中)	0.329	138	568	493	83	162	26	5	9	225	55	8	7	2	3	67	4	3	69	7	0.456	0.41
4	青木　宣親	(ヤ)	0.327	127	567	495	85	162	37	3	10	235	67	3	4	0	2	51	1	19	48	13	0.475	0.409
5	アルモンテ	(中)	0.321	132	546	498	56	160	37	0	15	242	77	1	1	0	3	44	1	1	95	16	0.486	0.375
6	鈴木　誠也	(広)	0.32	124	520	422	86	135	32	2	30	261	94	4	4	0	5	88	2	5	116	4	0.618	0.438
7	雄平	(ヤ)	0.318	125	482	446	50	142	19	1	11	196	67	6	1	0	2	33	5	1	62	13	0.439	0.365
8	宮﨑　敏郎	(デ)	0.318	142	590	551	71	175	34	0	28	293	71	0	0	0	0	38	3	1	45	16	0.532	0.363
9	坂口　智隆	(ヤ)	0.317	139	595	508	64	161	22	4	3	200	37	9	7	7	2	75	5	3	60	15	0.394	0.406
10	山田　哲人	(ヤ)	0.315	140	637	524	130	165	30	4	34	305	89	33	4	0	3	106	5	4	119	8	0.582	0.432
11	ソト	(デ)	0.31	107	459	416	74	129	16	0	41	268	95	0	1	0	5	29	2	9	100	5	0.644	0.364
12	岡本　和真	(巨)	0.309	143	616	540	82	167	26	0	33	292	100	2	1	0	0	72	1	4	120	11	0.541	0.394
13	糸井　嘉男	(神)	0.308	119	509	419	60	129	24	0	16	201	68	22	3	0	5	77	4	8	63	9	0.48	0.42
14	丸　佳浩	(広)	0.306	125	566	432	109	132	22	0	39	271	97	10	10	0	1	130	8	3	130	5	0.627	0.468
15	松山　竜平	(広)	0.302	124	446	397	46	120	25	2	12	185	74	2	0	0	5	42	4	2	46	10	0.466	0.368
16	筒香　嘉智	(デ)	0.295	139	580	495	77	146	33	1	38	295	89	0	0	0	3	80	7	2	107	10	0.596	0.393
17	ロペス	(デ)	0.288	110	459	441	46	127	23	0	26	228	77	0	0	0	1	16	1	1	58	12	0.517	0.314
18	野間　峻祥	(広)	0.286	126	447	405	64	116	14	7	5	159	46	17	10	4	2	30	0	6	69	2	0.393	0.343
19	糸原　健斗	(神)	0.286	143	637	531	79	152	29	4	1	192	35	6	4	9	4	86	1	7	73	5	0.362	0.39
20	マギー	(巨)	0.285	132	547	499	65	142	28	0	21	233	84	2	0	0	6	41	1	1	92	12	0.467	0.336
21	福留　孝介	(神)	0.28	123	499	414	57	116	26	2	14	188	72	2	1	0	7	73	2	5	90	10	0.454	0.389
22	大島　洋平	(中)	0.274	141	645	588	92	161	20	7	7	216	57	21	9	1	5	47	2	4	80	2	0.367	0.329
23	バレンティン	(ヤ)	0.268	142	602	514	72	138	22	0	38	274	131	1	1	0	3	85	5	0	121	15	0.533	0.37
24	田中　広輔	(広)	0.262	143	675	572	92	150	19	10	10	219	60	32	13	6	5	75	3	17	118	6	0.383	0.362
25	福田　永将	(中)	0.261	133	487	440	50	115	22	1	13	178	63	0	0	0	5	37	1	5	119	13	0.405	0.322
26	梅野　隆太郎	(神)	0.259	132	455	386	45	100	27	1	8	153	47	5	1	28	1	39	3	1	67	6	0.396	0.328
27	高橋　周平	(中)	0.254	128	477	433	35	110	26	2	11	173	69	0	0	8	3	30	1	3	89	9	0.4	0.305
28	亀井　善行	(巨)	0.254	123	461	422	47	107	20	0	13	166	49	4	0	0	1	36	2	2	70	6	0.393	0.315
29	西浦　直亨	(ヤ)	0.242	138	550	479	57	116	28	0	10	174	55	1	2	20	3	40	1	8	88	10	0.363	0.309
30	京田　陽太	(中)	0.235	143	632	578	73	136	15	7	4	177	44	20	10	26	3	19	0	6	111	12	0.306	0.266
31	菊池　涼介	(広)	0.233	139	642	557	85	130	27	1	13	198	60	10	2	30	1	51	2	3	111	5	0.355	0.301

これをRのstatsパッケージのhclust関数を用いてデフォルトで行いました。デフォルトではmethodは最遠隣法だそう。dist関数はデフォルトでユークリッド距離です。デンドログラムを作成すると

番号じゃなくて選手名にしたいのですがやり方が分からないので、取得次第修正します。1番近いのは3番と4番の選手のようです。これは平田選手と青木選手で、成績をみると確かに似ていました。ちなみにこの番号は打率の順番になっているのですが、クラスターになる際に近い番号同士になっている傾向がある気が僕はします。打率と選手の総合値にはそこそこの相関があるのでしょうか。それはそう。

パ・リーグ打者のクラスター分析

同様にパ・リーグの規定打席に到達した打者もクラスターにかけます。これもデフォルトで。

位	選　手	チーム名	打率	試合	打席	打数	得点	安打	二塁打	三塁打	本塁打	塁打	打点	盗塁	盗塁刺	犠打	犠飛	四球	敬遠	死球	三振	併殺打	長打率	出塁率
1	柳田　悠岐	(ソ)	0.352	130	550	475	95	167	29	5	36	314	102	21	7	0	5	62	4	8	105	8	0.661	0.431
2	秋山　翔吾	(西)	0.323	143	685	603	107	195	39	8	24	322	82	15	10	0	1	77	4	4	96	3	0.534	0.403
3	近藤　健介	(日)	0.323	129	555	462	59	149	29	3	9	211	69	5	0	0	5	87	8	1	90	16	0.457	0.427
4	吉田　正尚	(オ)	0.321	143	598	514	77	165	37	2	26	284	86	3	1	0	8	69	10	7	74	9	0.553	0.403
5	浅村　栄斗	(西)	0.31	143	640	565	104	175	27	0	32	298	127	4	2	0	5	68	4	2	105	18	0.527	0.383
6	中村　晃	(ソ)	0.292	136	580	506	57	148	28	1	14	220	57	1	1	3	6	60	4	5	68	5	0.435	0.369
7	井上　晴哉	(ロ)	0.292	133	548	476	59	139	26	2	24	241	99	1	0	0	6	63	0	3	106	7	0.506	0.374
8	島内　宏明	(楽)	0.292	103	453	394	53	115	16	3	11	170	53	11	5	3	3	47	0	6	45	10	0.431	0.373
9	外崎　修汰	(西)	0.287	119	510	453	70	130	24	3	18	214	67	25	9	6	1	47	0	3	102	9	0.472	0.357
10	中村　奨吾	(ロ)	0.284	143	639	552	82	157	30	3	8	217	57	39	15	0	5	60	2	22	94	11	0.393	0.374
11	山川　穂高	(西)	0.281	143	647	541	115	152	24	1	47	319	124	0	0	0	2	88	2	16	138	5	0.59	0.396
12	西川　遥輝	(日)	0.278	140	636	528	90	147	25	6	10	214	48	44	3	7	2	96	1	3	103	1	0.405	0.391
13	源田　壮亮	(西)	0.278	143	666	594	92	165	27	9	4	222	57	34	8	14	6	48	0	4	101	7	0.374	0.333
14	銀次	(楽)	0.276	139	552	492	45	136	16	5	5	177	48	1	1	6	2	48	3	4	47	7	0.36	0.344
15	今江　年晶	(楽)	0.276	127	465	421	44	116	17	2	10	167	49	0	3	0	4	32	5	7	76	19	0.397	0.334
16	森　友哉	(西)	0.275	136	552	473	67	130	34	2	16	216	80	7	2	0	7	70	2	2	105	3	0.457	0.366
17	上林　誠知	(ソ)	0.27	143	608	551	88	149	26	14	22	269	62	13	4	17	3	30	4	7	117	2	0.488	0.315
18	鈴木　大地	(ロ)	0.266	143	558	477	44	127	27	6	8	190	49	8	4	14	6	44	2	17	55	11	0.398	0.346
19	角中　勝也	(ロ)	0.265	112	470	411	44	109	23	2	7	157	57	3	3	0	6	48	2	5	63	9	0.382	0.345
20	中田　翔	(日)	0.265	140	599	540	61	143	32	0	25	250	106	0	1	0	13	43	1	3	81	24	0.463	0.316
21	田中　和基	(楽)	0.265	105	465	423	67	112	11	1	18	179	45	21	6	4	1	37	0	0	101	4	0.423	0.323
22	中島　卓也	(日)	0.261	132	449	391	57	102	11	3	1	122	23	29	5	22	1	35	0	0	88	4	0.312	0.321
23	松田　宣浩	(ソ)	0.248	143	580	517	72	128	21	3	32	251	82	3	3	0	3	56	4	4	113	12	0.485	0.324
24	田村　龍弘	(ロ)	0.239	143	476	415	32	99	14	7	3	136	35	3	3	16	2	40	1	3	68	13	0.328	0.309
25	デスパイネ	(ソ)	0.238	116	469	407	62	97	15	1	29	201	74	0	0	0	3	55	4	4	93	13	0.494	0.333
26	ロメロ	(オ)	0.237	119	501	443	63	105	18	1	25	200	63	7	2	0	6	42	2	10	114	10	0.451	0.313
27	レアード	(日)	0.233	120	505	450	47	105	14	2	26	201	65	0	1	0	4	44	5	7	124	10	0.447	0.309
28	藤岡　裕大	(ロ)	0.23	143	611	535	58	123	15	5	5	163	42	14	13	26	1	44	1	5	97	9	0.305	0.294
29	安達　了一	(オ)	0.219	140	514	465	44	102	12	3	3	129	41	20	5	16	5	25	1	3	64	8	0.277	0.261

26番と27番が最も近そう。ロメロとレアード(納得)。7と16はアジャと西武の森(納得)。

NPBの打者を別の距離の計算法を用いてクラスター分析をする

まずはセ・リーグの選手に対してmethodを="average"にして群平均法にしてみると

並べてみると

completeかaverageかは群が大きくなるにつれて違いが生じるので、下の段階では房の付き方は似ているが上に行くほど違う構図がみられるのではないでしょうか。

重心法にすると

ちょっと意味が分からないです。

最後によく聞くウォード法で

 一応全てscale()関数で標準化しています。果たしてどれがいいんでしょうか、やっぱりウォード法ですかね。

dendrogram for Central league and Pacific league

2018年のMLB打者を主成分分析してみた。

主成分分析とは

　主成分分析とは多くのデータのまとめる分析方法です。難しい言葉でいうと、複数の変数を1次式のデータに収縮させる事です。

　本当は自学自習の為に簡単な説明を書きたいのですがそんな技力も暇力もないので工事中(施工日未定)

2018年MLBの打者データを主成分分析

 G~SOまでの13変数のデータを主成分分析にかけて1つの式に説明していきます。

www.baseball-reference.com

 試合数・打席数・打数・得点・安打・2塁打・3塁打・本塁打・打点・盗塁・盗塁死・四球・三振の13変数で主成分分析をかけると

 こんな感じに。

第一主成分=-0.33*G+(-0.32)*PA+(-0.32)*AB+・・・・・+(-0.27)*SO

第二主成分=(-0.01)*G+・・・・・+0.43*X3B+・・・・+  (-0.10)*SO

summary()関数で

こんな感じに。

上から順に標準偏差・寄与率・累積寄与率となります。(空いた時間でこの説明もしたい所)

第二主成分までで82％が表現できているので今回は第二主成分まで考えます。

全部の係数が負なのですがこれは負荷量の絶対値が重要なわけで、全部負なのであれば逆転して考えればいいはず。

第一主成分では値が高い程良い変数も悪い変数も全て同程度に考えられています。良い意味でも悪い意味でも1年間を通してどれ程目立ったかを示す打者の総合値でしょう。

第二主成分のうち係数の絶対値が大きくて重要そうなのは

X3B(3塁打)の0.43,HRの-0.29,RBI(打点)の-0.23,SB(盗塁成功)の0.55,CS(盗塁死)の0.56です。これより第2主成分は1年を通して足の速さを印象づける指標といえるでしょう。

そこで、各選手の得点をplotしていくと

 選手が多すぎて汚い。特定の選手を抽出する方法も知らないし。なんとなく目に留まったのが1291番と316番だったので調べてみるとそれぞれマレックス・スミス選手とクリス・デービス選手で、マレックス・スミス選手の10本の三塁打と12個の盗塁死はリーグ最高で、40盗塁を決めているスピード選手です。クリス・デービス選手は48本の本塁打王で123打点を打つ強打者で、盗塁盗塁死は0で典型的なパワーヒッターで有る事が分かります。

もう少し見やすく有ってほしいので

規定打席に到達した選手で再び主成分分析

をかけてみると

第一主成分と第二主成分の係数は変れど意味は同じでしょう。第一は総合値で第二は走力。

各選手のplotは

となりました。今回でも第２主成分に関しては1291(マレックス・スミス)と316(クリス・デービス)が目立ってます。PC1で最高得点を取っているのは791番の選手です。こちらはフランシスコ・リンドーア選手でリーグ最多の129得点と打率.277、38本塁打、92打点、25盗塁を記録しています。総合値が高くてめっちゃいい選手ですね。ただし第2主成分まででは60.4％までしか説明できない所をつかれると痛いですね。

Principal Component Analysis for MLB2018 batter

MLBとコンカフェに通うヲタクで理解する相関係数

 相関係数とは２つの値の関連の大きさを測る尺度です。一般的には大きい正の値になるほど正の相関を持ち、大きい負の値に成る程負の相関を持ちます。比例と反比例みたいな感じ。実は何種類かあるんです。

 Pearsonの相関係数

　なんといってもこいつです。誰もが幼い頃に習ったこいつは鉄板の相関係数でしょう。

      

　rが0に近いと相関関係は無い

　rが+1に近い程正の相関があり

　rが⁻1に近い程負の相関が有る

　正の相関~解釈添え~

　2018年のMLBにおいて20本以上のヒットを打った打者のうちx軸にヒット数、y軸に打数をplotしました。当然ヒットを打つには多くの打数が必要なので正比例します。これを正の相関と言います。Pearsonの相関係数を計算すると 0.9730145を示します。

　計算する式の分母は２乗の値であるので必ず正であり尺度の大きさを調整するものです。分子に注目すると(xi-x)(yi-y)となっています。つまり、平均からの差の大きさと符号が考慮されています。

　上の例の図からわかるように正の相関を示すという事はplot点が第1象限と第3象限に集積している事を裏返しで意味しています。この象限にいる点は分子の値が正となるので相関係数rも正を示します。

　負の相関

　負の相関とは第２象限と第４象限に点が集積している状態で、この場合各点で与えられる分子は負を示すので相関係数はrを示します。

　本当はMLBデータでやりたかったのですがいい負の相関が得られそうになかったので適当にpythonで作りました。

　相関無し

　2018年に規定打席に到達したMLBの選手の年齢と安打数のplotをするとこんな感じ。

異論はあるかもしれませんが年齢と安打数はあまり関係がなさそうです。打つ奴は何歳でも打つし打たない奴は何歳でも打たない。この相関係数は-0.1478702でした。

MLB correlation coefficient ~positive for H and AB ...

Spearmanの順位相関係数

　ノンパラメトリックです。n個のベクトル値(x,y)があるとします。i個目のベクトル

(xi,yi)のx,yの順位をそれぞれsi,tiとします。

 　分母は定数ですよね。分子は各点の平均とのズレを表しています。考え方はPearsonの相関係数と同じで、平均からの差の正負と値の大きさを包含しています。２変数が共に平均より大きいor小さいなら順位に正の相関が有りρは1に近づきそうです。逆に２変数の平均からのずれが各々正と負である点が多い場合(順位に負の相関)はρは-1に近づきそうです。無関係ならρは0になりそう。

　これを式変形すると(これを章末問題としよう)(章末問題にはページ数のみ)(解答は省略)

　となります。実際はプログラミングで計算してもらうと思うのでこの値の意味が分かる変形前の方が僕は好きなんですけれども、変形後の方がよく見かけます。 

　ぴんと来てないと思うので実際のデータで考えましょう。ヲタクがコンカフェの店舗に行って使ったコンカフェ代とチェキ代で考えます(コンカフェ代とは全額からチェキ代を引いた額とします)。

ヲタク	コンカフェ代	チェキ代
1	1800	2000
2	2600	3000
3	1000	8000
4	11800	1000

　ヲタク1は学生でしょうか。推しに使えるお金は少なくて他の強ヲタに引け目を感じている、といったところでしょうか。僕はこの類です。ヲタク2はちゃんとバイトして推し事をする学生か普通の社会人の方でしょう。偉い。ヲタク3はチェキが大好きなヲタクです。コンカフェの食事には目もくれず推しとのチェキを撮りまくるタイプ。帰りに仲間と鳥貴族でビールにチェキを浸しがち。4はコンカフェの特別メニューをめっちゃ頼むヲタクです。

　これを順位にすると

ヲタク	コンカフェ代	チェキ代
1	3	3
2	2	2
3	4	1
4	1	4

  となります。n=4なのでn^3-nは60. Σ(si-ti)^2は18. よってρ=1-108/60=-0.8.よって負の順位の相関があると考えられます。つまり、コンカフェ代を払う人ほどチェキ代は払わない傾向にあるという事です。

Kendallの順位相関係数

 　n個のベクトル(x,y)が与えられた時、取り出した2つの(x,y)のベクトル値の大小関係が一致している個数をK,不一致の個数をLとした時の相関係数は以下のように与えられます。よって正の相関を示すほどτは1に近づき、負の相関を示すほど-1に近づきます。

　これもヲタクのコンカフェ代とチェキ代で考えます。

　n=4のベクトル値があるので取り出す選び方は4C2で6通り。それぞれヲタクの番号で考えると

　(1,2),(1,3),(1,4),(2,3),(2,4),(3,4)です。

　例えば(1,2)のヲタクの各値は

(コンカフェ代,チェキ代)=(1800,2000),(2600,3000)で、コンカフェ代とチェキ代の大小関係は一致しています。従ってKに+1をします。

(1,3)ではどうでしょう。(1800,2000),(1000,8000)でコンカフェ代とチェキ代の大小関係は一致していません。従ってLに+1を加えます。この操作を６回分行うと

　K=1,L=5となります。代入すると

τ=-4/6=-0.66666　となります。よって負の順位の相関関係があることが読み取れます。

MLB2018の安打数は正規分布を描いているのかQ-Qplotで考察する

 　今回はRを用いてQ-Qplotに関する話をしたいと思います。Q-Qplotは準一級の範囲だと思います。統計好きの人が問題にしたそうな内容ですよね。知りませんが。

Q-Qplotとは？

　q-qplotとは簡単に言うと２つの分布（多くは得られた分布と理想的な分布）の類似の程度を視覚的に把握するplot図の事です。右肩上がりの直線に成る程良いとされています。使う機会を見かけるのは「この得られた分布は正規分布と言っていいのだろうか？」知りたい時とかですかね。

Q-Qplotの実例

　今回はQ-Qplotの中でも得られたデータが正規分布を示しているのかを検証します。

こぼ場合はqqnormを使います。

得られたデータが正規分布している時のイメージ

step1 得られた分布を累積密度関数に

　標準正規分布(mean=0,sd=1)に従う分布から100個の乱数を抽出し、得られた分布がきちんと正規分布に従っているかを確認しましょう。

　まずはN(0,1)に従う分布から100個の乱数を抽出します。

　それをx軸に並べ、y軸は累積密度にしたグラフが以下です。今回のサンプルサイズは100なので、1値毎に累積密度は0.01だけ上がっていくので階段状になっています。

　負と正が綺麗に50-50に分かれているのが怖いですが抽出された100個の乱数を昇順に並べると以下の通りです。x=-2.36が最も小さい値で1つ目なのでy軸の累積密度関数は0.01の値をとります。x=-2.09では0.02、x=-1.96では0.03…のようにplotしていき、最大値x=2.06で累積密度関数は1をとります。

step2 比較する関数を累積密度関数に

　今回は正規分布と比較するので、正規分布の累積密度関数を準備します。

　

　見たことある図だと思います。例えば、標準正規分布表によるとZ=1.96の面積は0.475ですよね。正規分布は対称なのでx=-1.96での累積密度関数は0.025を取っています。

step3 2つの累積密度関数において各分位数（同じ累積密度関数の値)における各xの値をプロット

　言葉ではわかりずらいので実際に視覚的に見てみましょう。

　左から順にstep1で得た100個の乱数の累積密度分布、step2で得た累積密度分布、そしてq-qplotです。

　q-qplotは同じ分位数当たりの値をカウントしていく訳です。例えば40個目の乱数、今回は乱数を100個抽出したのでつまり累積密度関数が0.4の時の各xの値は何でしょう？

　左のstep1においてｙ＝0.4と交わる点のxは-0.243です（図では小さくて分かりにくいですが、上図の乱数の昇順で40個目の値と同値なので⁾。

　真ん中のstep2でy=0.4と交わる点のx値といえば、標準正規分布表で面積が0.1になるZ値を負にしたものです。標準正規分布表を開くとZ=0.25で面積は0.0987, Z=0.26で面積は0.1026となっています。今回はZ=0.25を採用することにします。

　右のQ-Qplotでは縦軸にはstep1で得られた値で今回は-0.24が与えられます。横軸にはstep2で得た値で今回は-0.25が与えられます。この縦軸-0.24と横軸-0.25で与えられる点をplotします。

　もう1つやってみましょう。累積密度関数の値=0.05の時を考えます。step1の乱数側では昇順で5番目の値なので-1.62です。step2の正規分布の累積密度分布ではp=0.05でよく使うZ=1.64なので与えられる値は-1.64です。Q-Qplotで左から5番目の点を読み取ると座標は(-1.64,-1.62)です。

　これを累積密度関数値=,0.050.4だけでなく0.01から0.02,0.03…0.98,0.99,1.00の100個分でプロットすることでQ-Qplotが完成します。

Q-Qplot for rnorm(100) and normal distribution

Q-Qplotの見方

　しつこいようですがQ-Qplotでは同じ分位数(累積密度関数の値)において求められる各数値を取ってplotしたものです。

　今回は同じ正規分布から得たので数値のスケールは同じですが、普段は異なる機会が多いでしょう。ですがQ-Qplotでみたいのは2つの分布が類似しているかどうかだけ。重要なのは散らばり具合の類似度な訳です。同じ分位毎にplotしているので、どちらの値も尺度は違えど全体のスケールからみて同程度の値を取っていればy=xに近づきそうですよね。またy=xからの解離・ゆがみ具合でどのように分布が歪んでいるかを把握することが出来ますし、問題として作り易そうです。

（暇な時に問題作りたい)

MLB選手の安打数は正規分布を示しているかをQ-Qplotで検証

　2018年に出場したMLB選手の安打数は正規分布を示しているのかを検証したいと思います(してないと思うので失敗例として扱う事を前提にしています)。本当はNPBにしようとしたけどMLBの方が人数多いしデータ揃ってそうなので。ここから引用しました。選手の抽出や.CSVでdataを得られたり本当に便利。すごい。

www.baseball-reference.com

step1 2018年のMLB選手の安打数を累積密度分布に落とす

　まずはヒストグラムを描くと

　0~20本が多くて明らかに正規分布ではないです。これを累積密度分布に落とし込むと

　0が多すぎて初っ端から累積密度が0.4位から始まるという笑

どうしてこうなるかというと毎日レギュラーで活躍している人は一握りで、4割近くの人はマイナーを往ったり来たりする控え選手or投手だからでしょう。

 これでいいのかと考えた私はヒット数が20本以下は切り捨てて、21本以上の人に絞り込みました（本当は規定打席の方がいいんでしょうが)。すると

　正規分布ではないですが、まだ使えそうな分布が得られました。累積密度分布に変換

すると

y=√x ,log(x)のような関数になりました。せっかくなのでこの2つをそれぞれQ-Qplotで正規分布と比較していきます。

step2 正規分布を累積密度分布に

　これは前のと同じなので省略

step3 2つの累積密度関数において各分位数（同じ累積密度関数の値)における各xの値をプロット

　まずはヒット数で限局していない1つ目の累積密度分布と正規分布でかけてやると

　　

　直線とは程遠い曲線が得られました。3つ並べると

　左がヒット数の累積密度分布、真ん中が標準正規分布の累積密度分布、右がQ-Qplot.横軸が正規分布のx値で縦軸がMLB側のx値つまりヒット数です。初めは横軸が進んでも縦軸が全く縦に伸びていないのは0付近の値が多いからです。つまり半数近くは0付近のヒット数を誇っているのでしょう。

　次に、ヒット数21本以上に限局した場合のQ-Qplotでは

直線とは言い難いですが、1つ目の限局しない場合よりはマシだと思います。それでも極端でないとはいえ本数が少ない選手の方が相対的に多くなっているので最初は鈍く、途中から鋭く曲がるような曲線になってしまっています。

　よってMLBの選手の安打数は正規分布していないことが分かります（この程度だとhistで分かりますが）。

Q-Qplot as for the number of hit at(BAT) 2018MLB

Q-Qplotの正確な理解(未)

　上で説明したのが１番イメージしやすい考えだと思っているし、Q-Qplotってそんなに正確な情報を得るために使う物ではないと勝手に思っているので十分だと思うのですが、実際は違うみたいなんですよね～　頭いい人なら余裕で理解できるのでしょうが

　だってこの考えだと累積密度関数の値が1の時の正規分布側の値って無限になりません？　ずっと引っかかってました。wikipediaを拝見した所、

　正規確率プロットを使用する際には、標準正規分布の順序統計量の期待値の尺度であるランキットを使用する。 

引用元　Q-Qプロット - Wikipedia

　なんとなくわかった。でも確認のために標準正規分布の順位統計量を自力で求めてみればいいのだろうけどちょっと面倒くさい…それにRでhelp(qqnorm)を入れても全部英語だし。上で説明したのも的外れ？ 今後も勉学に努めます…

青木宣親から学ぶ 野球は統計学で説明できるのか？~2項分布編~

　みなさんは青木宣親選手をご存知でしょうか？

　Noriこと青木宣親選手は日本だけでなくMLBでも活躍できた日本の野球史に名を残す現役の名選手です。

　青木選手は1982年に宮崎県日向市で産声を上げました。

　それからというものの、早稲田時代ではリーグ4連覇を果たし、ヤクルトに入団後も順調に成績を重ね首位打者や最多安打・盗塁王といったタイトルを総なめにし、アナウンサーと結婚。そしてMLBに移籍しました。福留孝介選手や松井稼頭央選手など多くの名選手はMLBの成績が不安定で大幅に落としていた中で青木選手はMLBでも着実に成績を残していきました。その成績は

　2012  151試合 打率.288
　2013  155試合 打率.286
　2014  132試合 打率.285
　2015  93試合   打率.287
　2016  118試合 打率.283
　2017  110試合 打率.277 

　　　　通算打率.285

　ヒトは青木をこう呼びます

　ミスター2割8分

　2割8分に収束する男

　アヘアへ0.280マン

　青木宣親選手はシーズン中に波はありました。

　今年は無理だろうと巷で交わされていた時も数知れず。

　なのに終わってみると必ず2割8分台に乗せているのです。

　こんなことできるでしょうか？

　多くの選手は毎年少なからず成績が変動します。成績は体調面だけでなく私生活やモチベーション、相手チームの分析など多岐に依存するからでしょう。

　この異端な成績を統計学の側面から考えるとどのように評価できるのでしょうか？というのが今回の議題です。

 　野球と統計学

　ここに野球と統計の歴史を書こうと思ったのですが（マネーボールとか）、まあ要らないかなと思うので

　工事中(着工日未定)

　打率を統計学で考える

　今回は青木選手が毎年同じような打率で安打を重ねる事は統計学的にどういう意味を持つのかを考えます。ですが複雑に考えると大変になるので状況は簡潔にします。青木選手はどのような投手・状況・場面でも0.285の確率(打率)で安打を打つと仮定します。野球が好きな人ほど抵抗がある理論だと思いますが、打数が多くなれば案外イイ感じになるんです。

　この仮定において今回は、1年間の打率を2項分布で考えることが出来ます。

　この青木選手が1年間で1打席だけ立ったシーズンを1000回繰り返した(1000年間シーズンを過ごす)場合、1000シーズン分の打率はどうなるでしょうか？

このplotだけでは正確な数字は分かりませんが、1000×0.285=285年は打率が10割で715年は打率が0割と言う期待値に近いことは読み取れます。

※　table(a)で正確な数字を求めたところそれぞれ740と260でした。

では1シーズンで10打席に立った青木選手の年間打率を求める事象を1000回繰り返した場合、各シーズンの打率の値はどうなるでしょうか？

打率	シーズン数
0	34
0.1	142
0.2	280
0.3	242
0.4	175
0.5	96
0.6	26
0.7	2
0.8	3

 こんな感じです。

aoki_10bats added way to make table for html style ...

10打席しかないので打率も0から1まで0.1刻みですが流石に9割打者と10割打者は居ませんね。

次は、青木選手が1シーズンで100打席に立ったシーズンの打率を求める試行を1000回繰り返す(1000シーズン分)った時の分布は

打率	0.16	0.17	0.19	0.2	0.21	0.22	0.23	0.24	0.25	0.26	0.27	0.28	0.29	0.3	0.31	0.32	0.33	0.34	0.35	0.36	0.37	0.38	0.39	0.4	0.41	0.42	0.43	0.44	0.46
シーズン数	1	2	2	8	11	28	31	34	45	62	79	80	70	83	85	72	68	53	45	41	27	19	14	17	14	2	4	2	1

  打率0.31が85シーズンで最頻値となっています。一応、中心極限定理で正規分布を示していますが。

aoki_100bats added ways to change rows and comns

　では今回のメインです。1年間怪我無く無事に送った場合に立つ打席数は500打席程度でしょう。ですので、青木選手が1シーズンで500打席に立ったシーズンの打率を求める試行を1000回繰り返す(1000シーズン送)った時の分布は

	0.22	0.222	0.226	0.23	0.232	0.234	0.236	0.238	0.24	0.242	0.244	0.246	0.248	0.25	0.252	0.254	0.256	0.258	0.26	0.262	0.264	0.266	0.268	0.27	0.272	0.274	0.276	0.278	0.28	0.282	0.284	0.286	0.288	0.29	0.292	0.294	0.296	0.298	0.3	0.302	0.304	0.306	0.308	0.31	0.312	0.314	0.316	0.318	0.32	0.322	0.324	0.326	0.328	0.33	0.332	0.334	0.336	0.34	0.348
1	1	1	1	1	3	3	2	3	3	2	7	7	7	4	9	12	4	19	17	19	15	22	31	34	29	38	42	30	56	43	44	43	35	34	35	40	31	29	37	21	26	18	20	19	15	16	13	11	10	6	5	10	8	2	1	2	2	1	1

  期待値となる打率の0.282 , 0.284になるシーズン数はそれぞれ43,44シーズンで、2割8分台に収束しているシーズン数は221シーズンでした。

aoki_500bats

　つまり本来の打率が0.285の選手でもせいぜい500打席程度では選手の真の能力を測るには不十分であり、ぶれが生じるという話です。試しに両側95％信頼区間を計算してみると0.244≦打率≦0.326 であります。このように500打席程度では本来の能力とはかけ離れた数値を叩き出すことも容易に起こるのです。

　この結果から打率を2項分布で考えるのがナンセンスだ、確率論は野球に応用できないという旨の反論も当然あると思います。

　が

　一方で、特定の季節・相手投手・球場　etcで調子が悪いというのも実はこのような2項分布の偶然のぶれに起因するものであるという考えもできます。

　最終年だけ活躍するおかわり君や夏に強いゴジラも偶然に起因した結果なのかもしれません。また、2軍の数十打席の数字だけで1軍に上げるかどうかを判断する事にも疑問符が呈されるでしょう。

　というように、この打率の考え方では1シーズン当たりの打率にぶれは欠かさないので、より選手の能力を測れる変数を探すことになりました。その結果OPSなどが作られて現在ではそれらにスポットライトが当てられている訳なのですが、青木選手はこの自然の摂理の予想に反した結果を残している訳です。

　500打席あれば本来の力の打率付近に収束する確率は高いし今の言動は大袈裟かもしれませんが、少なくとも青木選手のアヘアへ0.280マンは2項分布の野球への応用に改善の余地がある可能性を示すものとなっているでしょう。

　※まだ不十分なのは承知なので、今後も暇な時に野球と統計の話を追加していきます。