文系のための「正規分布」（２）

正規分布の存在意義というのは、結局のところ、平均と標準偏差さえ分かれば、
「釣り鐘状」に値が分布するような関数を作ることができ、その関数を用いることで、
ある値を取る確率を簡単に計算できるということであった。

もちろん、釣り鐘状には分布しないこともあって、
その時には、残念ながら正規分布を前提にすることはできないのであるが、
別の様々な分布を前提とすることになる。

他の様々な分布を知ることは重要ではあるが、ここでは正規分布に限った話にする。
まずは、正規分布の意味が解っていれば、他の分布関数を理解する上でも、
多少の手助けにはなるだろう。

さて。前回までの話では、一般的な正規分布の話をしてきたが、
今回の話では、少し特殊な場合での正規分布について整理する。

どのように特殊なのかというと、つまり、
平均が「0」で標準偏差が「1」となるような正規分布である。
この種の正規分布のことを、一般的には、「標準正規分布」と呼ぶ。

この標準正規分布を用いると、色々と計算が便利になるので、様々なとこに出てくる。
なぜ、標準正規分布を用いるか、という話はとりあえず置いておいて、まずは基本的な話から。

さて、前回の話では、一般的な正規分布の関数の式を以下のように表した。



ここで、μは平均を表し、σは標準偏差を表したのであった。
標準正規分布では、μ=0、σ=1 となるので、以下のように表される。



確か、正規分布は、平均と標準偏差を与えれば良いだけなので、
これだけで、釣り鐘状の山を描くことができる。Rでは以下のようにする。

# Xの取り得る範囲を決める。とりあえず、-4～4くらいで。これはテキトー。
X.range <- seq(-4, 4, 0.01)

# 標準正規分布のX.range の間の f(x)の値を得る
X.dnorm <- dnorm(X.range)

# 山を描く
plot(X.range, X.dnorm, type="l", xlab="", ylab="", ylim=c(0,0.5))
title(main="Bell curve for Standard Normal Distribution")
title(xlab="Range of x value(Standard deviation)", ylab="Density")
abline(h=0)

# 平均を通る線を引く
abline(v=0, lty=2, col="red")
text(0, 0.45, expression(mu==0), cex=1.8)

この図において、x軸の値が標準偏差であり、y軸の値がその時の確率の密度を表す。
標準偏差が「0」の所で山の頂点があり、ここが平均を示している。
また、「±4」付近に山裾が存在しているのが確認できる。

標準正規分布において重要なことは、
平均からの距離を標準偏差の倍数で確率を推測するということであり、
区切りの良い距離として±1σ倍、±2σ倍、±3σ倍を用いることが多い。

とりあえず、この状況を可視化して整理してみる。まずは、±1σ倍の場合から。

なお、今回から、x軸の軸ラベルには、一般的な表記を用いるが、
平均「0」、標準偏差「1」の場合には、先ほどの図と同じ状況になる。

# 山を描く
plot(X.range, X.dnorm, type="l", xlab="", ylab="", xaxt="n", ylim=c(0,0.5))
axis(1, at=seq(-3, 3, by=1), labels=c("μ-3σ","μ-2σ", "μ-1σ", "μ", "μ+1σ", "μ+2σ", "μ+3σ"))
title(main="Bell curve for Standard Normal Distribution")
title(xlab="Range of x value(Standard deviation)", ylab="Density")
abline(h=0)

# 面積の部分を塗りつぶす
poly.x <- seq(-1,1,0.01)
poly.y <- c(0, dnorm(poly.x), 0)
polygon(c(-1, poly.x, 1), poly.y, col="lightblue")

# 平均を通る線を引く
abline(v=0, lty=2, col="red")

# キャプションを入れる。
prob <- paste(round(pnorm(1) - pnorm(-1), 3)*100, "%")
text(0, 0.1, prob, cex=1.8)

次は、±2σ倍の場合。

# 山を描く
plot(X.range, X.dnorm, type="l", xlab="", ylab="", xaxt="n", ylim=c(0,0.5))
axis(1, at=seq(-3, 3, by=1), labels=c("μ-3σ","μ-2σ", "μ-1σ", "μ", "μ+1σ", "μ+2σ", "μ+3σ"))
title(main="Bell curve for Standard Normal Distribution")
title(xlab="Range of x value(Standard deviation)", ylab="Density")
abline(h=0)

# 面積の部分を塗りつぶす
poly.x <- seq(-2,2,0.01)
poly.y <- c(0, dnorm(poly.x), 0)
polygon(c(-2, poly.x, 2), poly.y, col="lightblue")

# 平均を通る線を引く
abline(v=0, lty=2, col="red")

# キャプションを入れる。
prob <- paste(round(pnorm(2) - pnorm(-2), 3)*100, "%")
text(0, 0.1, prob, cex=1.8)

次は、±3σ倍の場合。

# 山を描く
plot(X.range, X.dnorm, type="l", xlab="", ylab="", xaxt="n", ylim=c(0,0.5))
axis(1, at=seq(-3, 3, by=1), labels=c("μ-3σ","μ-2σ", "μ-1σ", "μ", "μ+1σ", "μ+2σ", "μ+3σ"))
title(main="Bell curve for Standard Normal Distribution")
title(xlab="Range of x value(Standard deviation)", ylab="Density")
abline(h=0)

# 面積の部分を塗りつぶす
poly.x <- seq(-3, 3, 0.01)
poly.y <- c(0, dnorm(poly.x), 0)
polygon(c(-3, poly.x, 3), poly.y, col="lightblue")

# 平均を通る線を引く
abline(v=0, lty=2, col="red")

# キャプションを入れる。
prob <- paste(round(pnorm(3) - pnorm(-3), 3)*100, "%")
text(0, 0.1, prob, cex=1.8)

さて、ところで、平均「0」、標準偏差「1」、という状況から連想することはないだろうか？
確かにどこかで聞いたような。そうそう、標準化の話に似ている。
標準化の話では、標準偏差ではなく、分散が「1」と説明したが、結局は同じこと。

念のために。標準偏差を「」としたとき、分散は「」である。つまり、二乗ということ。
したがって、「」の場合、「」となる。

つまり、何が言いたいかというと、正規分布に従うデータがあって、
そのデータを標準化したデータは、標準正規分布に従うのである。
ここで、標準化の話を思い出してみる。確か、以下のような式であった。



重要なのは、右端の標準化の式である。
これを使えば、一般的な正規分布を標準正規分布に変換できる。

なぜ、これが便利なのか？

前回の話の最期に、数学嫌いの人にとって目を背けたくなる積分の式が登場した。
これを実際に「手計算」で解くのは、少々厄介なのである。
もっとも、現在は、コンピュータが計算してくれるので問題は無いのであるが...。

そういった状況があって、コンピュータが存在していない時には、
「標準正規分布表」というものを使っていた。

あらゆる正規分布の表を作ろうとすると、無限個の表を作らないといけないので、
要するに、標準正規分布の場合だけの表を作り、それを参照していたのである。
っということで、試しに標準正規分布表を作ってみる。

かなり、長くなってしまうので、マイナス方向は省略。平均を中心に対称なので構わない。
単なる表を作るだけの作業なので、余裕の無い人はパスして構わない。

# 行数を数えるための変数を初期化する
nrow <- 0

# 結果を格納するための行列を準備する
N <- matrix(ncol=10, nrow=41)

# 標準正規分布表を作成する
for(i in seq(0, 4, by=0.1)){
    nrow <- nrow + 1    # 行数を一行進める
    ncol <- 0               # 列数を数えるための変数を初期化する
    for(j in seq(0, 0.09, by=0.01)){
        ncol <- ncol + 1  # 列数を一行進める

        # 平均から上の部分の面積を計算する。これが確立。
        N[nrow, ncol] <- round(pnorm(i+j) -pnorm(0), 4)
    }
}

# 列名と行名を入れる

rownames(N)<-seq(0, 4, by=0.1)

colnames(N)<-seq(0, 0.09, by=0.01)

そして、これを実行すると、標準正規分布表が出来上がる。

> N

         0   0.01   0.02   0.03   0.04   0.05   0.06   0.07   0.08   0.09

0   0.0000 0.0040 0.0080 0.0120 0.0160 0.0199 0.0239 0.0279 0.0319 0.0359

0.1 0.0398 0.0438 0.0478 0.0517 0.0557 0.0596 0.0636 0.0675 0.0714 0.0753

0.2 0.0793 0.0832 0.0871 0.0910 0.0948 0.0987 0.1026 0.1064 0.1103 0.1141

0.3 0.1179 0.1217 0.1255 0.1293 0.1331 0.1368 0.1406 0.1443 0.1480 0.1517

0.4 0.1554 0.1591 0.1628 0.1664 0.1700 0.1736 0.1772 0.1808 0.1844 0.1879

0.5 0.1915 0.1950 0.1985 0.2019 0.2054 0.2088 0.2123 0.2157 0.2190 0.2224

0.6 0.2257 0.2291 0.2324 0.2357 0.2389 0.2422 0.2454 0.2486 0.2517 0.2549

0.7 0.2580 0.2611 0.2642 0.2673 0.2704 0.2734 0.2764 0.2794 0.2823 0.2852

0.8 0.2881 0.2910 0.2939 0.2967 0.2995 0.3023 0.3051 0.3078 0.3106 0.3133

0.9 0.3159 0.3186 0.3212 0.3238 0.3264 0.3289 0.3315 0.3340 0.3365 0.3389

1   0.3413 0.3438 0.3461 0.3485 0.3508 0.3531 0.3554 0.3577 0.3599 0.3621

1.1 0.3643 0.3665 0.3686 0.3708 0.3729 0.3749 0.3770 0.3790 0.3810 0.3830

1.2 0.3849 0.3869 0.3888 0.3907 0.3925 0.3944 0.3962 0.3980 0.3997 0.4015

1.3 0.4032 0.4049 0.4066 0.4082 0.4099 0.4115 0.4131 0.4147 0.4162 0.4177

1.4 0.4192 0.4207 0.4222 0.4236 0.4251 0.4265 0.4279 0.4292 0.4306 0.4319

1.5 0.4332 0.4345 0.4357 0.4370 0.4382 0.4394 0.4406 0.4418 0.4429 0.4441

1.6 0.4452 0.4463 0.4474 0.4484 0.4495 0.4505 0.4515 0.4525 0.4535 0.4545

1.7 0.4554 0.4564 0.4573 0.4582 0.4591 0.4599 0.4608 0.4616 0.4625 0.4633

1.8 0.4641 0.4649 0.4656 0.4664 0.4671 0.4678 0.4686 0.4693 0.4699 0.4706

1.9 0.4713 0.4719 0.4726 0.4732 0.4738 0.4744 0.4750 0.4756 0.4761 0.4767

2   0.4772 0.4778 0.4783 0.4788 0.4793 0.4798 0.4803 0.4808 0.4812 0.4817

2.1 0.4821 0.4826 0.4830 0.4834 0.4838 0.4842 0.4846 0.4850 0.4854 0.4857

2.2 0.4861 0.4864 0.4868 0.4871 0.4875 0.4878 0.4881 0.4884 0.4887 0.4890

2.3 0.4893 0.4896 0.4898 0.4901 0.4904 0.4906 0.4909 0.4911 0.4913 0.4916

2.4 0.4918 0.4920 0.4922 0.4925 0.4927 0.4929 0.4931 0.4932 0.4934 0.4936

2.5 0.4938 0.4940 0.4941 0.4943 0.4945 0.4946 0.4948 0.4949 0.4951 0.4952

2.6 0.4953 0.4955 0.4956 0.4957 0.4959 0.4960 0.4961 0.4962 0.4963 0.4964

2.7 0.4965 0.4966 0.4967 0.4968 0.4969 0.4970 0.4971 0.4972 0.4973 0.4974

2.8 0.4974 0.4975 0.4976 0.4977 0.4977 0.4978 0.4979 0.4979 0.4980 0.4981

2.9 0.4981 0.4982 0.4982 0.4983 0.4984 0.4984 0.4985 0.4985 0.4986 0.4986

3   0.4987 0.4987 0.4987 0.4988 0.4988 0.4989 0.4989 0.4989 0.4990 0.4990

3.1 0.4990 0.4991 0.4991 0.4991 0.4992 0.4992 0.4992 0.4992 0.4993 0.4993

3.2 0.4993 0.4993 0.4994 0.4994 0.4994 0.4994 0.4994 0.4995 0.4995 0.4995

3.3 0.4995 0.4995 0.4995 0.4996 0.4996 0.4996 0.4996 0.4996 0.4996 0.4997

3.4 0.4997 0.4997 0.4997 0.4997 0.4997 0.4997 0.4997 0.4997 0.4997 0.4998

3.5 0.4998 0.4998 0.4998 0.4998 0.4998 0.4998 0.4998 0.4998 0.4998 0.4998

3.6 0.4998 0.4998 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999

3.7 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999

3.8 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999

3.9 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000

4   0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000

少々見難い出力となってしまったが、一行目は行名で、一列目は列名である。

ちょっと、特殊な見方になっていて、行方向は小数点第一位を示していて、

列方向は小数点第二位を示している。

つまり、0.36σ のときの面積であれば、

「0.3」+「0.06」なので、4行目の7列目の値を見る。

すると、「0.1406」となっている。これが、平均「0」からの面積となる。

さて、ここで、1σ 分を見てみると、11行目の1列目の値なので「0.3413」となっている。

この値を2倍すると「0.6826」となり、±1σの値と一致する。

つまり、標準正規分布表は以下のような状況を示している。

# 片袖の標準正規分布の図を描く

x <- seq(0, 4, by=0.01)

y <- dnorm(x)

plot(x, y, type="l", xlab="", ylab="", xaxt="n")

title(main="Bell curve for Standard Normal Distribution")

title(xlab="Range of x value(Standard deviation)", ylab="Density")

abline(h=0)

# 面積の部分を描く

poly.x <- seq(0, 1.8, 0.01)

poly.y <- c(0, dnorm(poly.x), 0)

polygon(c(0, poly.x, 1.8), poly.y, col="lightgray")

text(0.8, 0.1, expression(integral(f(x)*dx,0,a)), cex=1.8)

現在では、すっかり、標準正規分布表を用いることは無くなったが、
かつては、この表を用いて確率を求めていたのである。

さて、重要なことは、要するに正規分布に従うあらゆるデータは、
そのデータを標準化することで標準正規分布に変換して考えることができる。

文系のための「正規分布」（１）

統計の最も基本的な知識でありながら、理解し難いものに「正規分布」がある。
「釣り鐘状の分布」などと言われているが、実際のところ、
何故、この分布が重要なのか？その理解に苦しむ人は多い。

例えば、庭の植木に生えている花の「花弁」を無作為に1,000枚採ってきたとする。
季節によっては、花など生えていないかもしれない。その時は同じ種類の「葉」でも良い。
とにかく、サイズを測りやすそうな植物の花弁か葉を1,000枚用意する。

そして、その長さを測るとしよう。まぁ、実際には幅でも良いのだが。

当たり前の話ではあるが、無作為に集めた1,000枚の花弁（葉）のサイズが
全て同じということはあり得ない。必ず、微妙な差が生じるのである。
では、この差がどのように現れるのか？

平均的なサイズが最も多く、それより大きくなるか、小さくなるほど数は少なくなる。
このような状況をヒストグラムで表現すると釣り鐘状に値が並ぶ。

このような分布をすることを「正規分布」と呼ぶ。

「自然界の至るところで見ることができる。」と、一般的には言われているが、
実際にそのように分布しているデータを探すのは意外に大変。
試験の結果も、身長のデータも、植物の花弁や葉も、実際には理論通りにはいかない。

では、不必要あるかというと、そういうわけでもない。

正規分布を前提とすることで、多くの統計的な処理が説明しやすくなり、
また、実用上の問題としては、正規分布に従っていなくても良いこともある。
もちろん、正規分布の前提が崩れると危険なものも存在するが…。

それでも、平均と標準偏差が分かれば、ある値の現れる確率を得れるという性質は重要。
それゆえに、現在においても、正規分布を前提とする手法は多い。

正規分布を盲目的に信用しそれを前提に置いてはならないという主張は確かにある。
しかし、だからと言って、知らないで良いわけでは断じてない！

さて、話は逸れたが、実際のデータを使いながら整理してみる。
Rには、様々なサンプルデータがあるので、今回はアヤメ（iris）のデータを使ってみる。
とりあえず、データを準備する。

# アヤメのデータを確認する。以下のコマンドを実行するだけ。
head(iris)

データの中身を見てみると、以下のようになっている。

> head(iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

Sepal はガクのことで、Petal は花弁のことである。
Speciesには、アヤメの種類が入っている。

まずは、全体の状況を可視化して確認したい。
beeswarm()関数とboxplot()関数を使って、アヤメの花弁の長さを可視化してみる。

library(beeswarm)
boxplot((iris$Petal.Length ~ iris$Species))
beeswarm((iris$Petal.Length ~ iris$Species), pch=16, col="red", add=TRUE)

三種類のアヤメ、setosa、versicolor、virginicaの花弁の長さの分布が示されている。

今回は、少々、変わった方法でデータを指定しているが、
これはモデル式と呼ばれる方法による指定。
Petal.Lengthという変数の値をSpeciesという変数で分けるという意味。

非常に有名なデータで、よく、チュートリアルデータとして用いられる。
このデータは、ロナルド・フィッシャーという有名な統計学者が、
判別分析という手法を紹介するためのサンプルデータとして紹介したものである。

元々は、エドガー・アンダーソンという人物が、
アヤメの形状を定量的に分析するために、
北アメリカのガスペ島という所で収集したものらしい。

とのかく、このデータは、3種類のアヤメについて、50株ずつ採集し、
ガクの長さと幅、花弁の長さと幅の4変数について測定されたものである。

さて、今回は、解りやすいようにsetosaの花弁のデータを用いることにする。
なんとなく、正規分布に従っていそうでもある。

X <- split(iris$Petal.Length, iris[,5])[1][[1]]
hist(X.sample)

ここでは、split()関数が新しく登場している。
この関数は、元のデータから、特定の列の値で切り分けるための関数である。

今回の場合は、花弁の長さのデータを、setosaで切り分けている。
最後に、[[1]]を付けているのは、list型からベクトルを取り出すためのオマジナイ。

さて、準備が出来たところで、早速、正規分布の特性について整理していく。

まず、正規分布は、ある値の発生する可能性を「確率」として表すことができ、
その確率を推定するための「関数」を「平均」と「標準偏差」によって決めることができる。

# setosaの花弁の長さの最小と最大をxの取る範囲とする
X.range = seq(min(X), max(X), 0.01)

# setosaの花弁の長さの平均と標準偏差で正規分布を描く
X.dnorm <- dnorm(seq(min(X), max(X), by=0.01),mean(X),sd(X))
plot(x=X.range, X.dnorm, type="l", col="red", xaxt="n", yaxt="n", ylim=c(0,2.6), xlab="", ylab="")

# setosaの元のデータの分布を重ね合わせる。
par(new=TRUE)
plot(density(X, cut=0), xlim=c(min(X), max(X)), ylim=c(0, 2.6))

Rでは、正規分布を描くためのdnorm()関数が用意されていて、
このパラメータに、x軸の値のベクトル、平均、標準偏差を与えてやると、
x軸の値に対応するy軸の値、すなわち、正規分布を描く関数によって導かれた値が出る。

今回は、x軸の値は、setosaの最小値〜最大値の範囲を等間隔に0.01刻みで与えている。
ある分位数を与えていると考えることもできる。Rのマニュアルではそのようになっている。
要するに、y軸の値を求めるときのx軸の値。刻み幅が小さいと鋭角な折れ線となる。

この図では、元のデータとの重ね合わせを行なっているが、今の段階では、
とりあえず、分布密度を取ったヒストグラムに対応するものと理解しておけば良い。
ヒストグラムとは異なり、ある関数を用いて推定を掛けている。

これを見てみると、赤い曲線で描いた正規分布の曲線と、
黒い曲線で描いた元のデータの推定は、なんとなく似た状況を示している。

さて、話を元に戻す。

そうそう、正規分布を描くと、ある値の取る可能性を知ることができる。
例えば、「花弁の長さが1.3センチよりも小さい」という状況は、
どの程度の確率で発生するのか？

そのようなことが、簡単に知ることができるのである。

イマイチ、ピンと来ないかもしれないが、
ある事象の起こり得る確率を求めるのは非常に難しい。
サイコロの目のように、離散的な値であれば良いが、
連続的な値の場合はもっと難しい。

正規分布というのは、「平均」と「標準偏差」だけ与えてやれば、
簡単に確率を教えてくれる。そういう、魔法がかかっているのである。

では、どうすれば、その「確率」とやらを教えてくれるのだろうか？
それを知るためには、まずは、正規分布の「関数」を知らなくてはならない。



ここからは、少し辛抱。厄介な式が出てくるが、要するに、「平均」と「標準偏差」、
そして、「xの値」が分かれば、釣り鐘状の山を描けるような関数であることが解れば良い。

唯一、意味不明なものが、「e」 の記号。これは、指数関数を表している。
説明をし始めると先に進まないので、とりあえずは、この記号は関数であって、
この関数のべき乗の部分が解れば、何かの値が出てくると理解しておく。

いずれ、指数関数の説明はしなければならない。

何度も繰り返すが、この関数は、平均と標準偏差だけが解れば良い。
この二つさえ知ることができれば、確率がでてくる。
xの値は分位数なので気にする必要は無い。

では、なぜ、確率が求められるかというと、
実は、この関数によって描かれる釣り鐘状の山の面積が「1」となるため。
つまり、面積の広さを確率として求めているのである。

したがって、確率を求めたいある範囲があって、それに対応するxの値が解れば、
この正規分布を描く関数を使って、面積を求めることができる。
そして、その面積こそが、真に知りたい確率となる。

ところで、正規分布において平均と標準偏差は何を意味するのか？

まず、この関数において平均は釣り鐘状の分布の最も高い地点に当たる。
そして、厳密に正規分布に従う場合この地点は中央値にも一致する。

そう言えば、対象数が多ければ、平均と中央値が近づくという話をした。
つまり、そのような状況は正規分布を前提としているのである。

一方、標準偏差というのは、要するに、バラツキの指標であって、
バラツキが大きいということは、それだけ、「山の裾」は広がる。

面積が「1」に固定された状態で、山裾が広がるということは、山の形は扁平になり、
逆に、標準偏差が低いということは、裾が狭まるので、山の形は切り立つ。

かなり、遠回りしたが、先程の例の場合の場合はどうなるか？つまり、
「花弁の長さが1.3センチよりも小さい」場合は、どのような確率であろうか？
実際には、かなり複雑な計算となるので、まずはRでの実行から。

pnorm(1.3, mean(X), sd(X))

Rでは、pnorm()関数を用いるとこの計算ができる。実行結果は以下の通り。

> pnorm(max(X), mean(X), sd(X))
[1] 0.1754524

この結果から、17.5%の確率で出現することが解った。
当然、この結果は、あくまで正規分布に従ったという前提の下ではあるが、
正規分布に従っているという前提を与えることで、容易に確率を求めることができる。

では、これはどのような状況なのか、可視化して確認してみる。

# 正規分布を描き直す
X.dnorm <- dnorm(seq(min(X), max(X), by=0.01),mean(X),sd(X))
plot(x=X.range, X.dnorm, type="l", col="red", xlab="", ylab="")

# タイトルと軸ラベルの設定
title(main="Probability of occurence for petal length within 1.3cm")
title(xlab="Length of Petal", ylab="density")

# 面積の部分を塗りつぶす
poly.x <- seq(1,1.3,0.01)
poly.y <- c(0,dnorm(poly.x, mean(X),sd(X)),0)
polygon(c(1,poly.x,1.3), poly.y, col="red")

# 理解しやすいようにオマケ
caption <- paste(round(pnorm(1.3, mean(X), sd(X))*100, 2), "%")
text(1.2, 0.25, caption, cex=1.3, col="white")

結局は、この赤色の面積のことを言っているのである。

ところで、この面積は、正規分布の積分で表され、
今回の例で言うと、1.3に当たる部分を一般化して「a」と置き、
「a」よりも小さい値が発生する確率を、 と置くと次のようになる。



積分は、丁度、f(x)とx軸で囲まれた部分の面積に相当するため、
この式によって、ある値「a」よりも下の値を取る場合の確率が求まる。
ちなみに、面積「1」の場合は、つまり全ての範囲は、以下のようになる。



実は、この積分を計算するのは、少々厄介で、しかも、どうにも解りにくい。
そこで、一般的には、標準正規分布と呼ばれる、
ある特殊な正規分布の形に変換してこの部分の面積を求めることができる。

この話については、次回に。

文系のための「二変数の関係」（１）

前回の話では、以下のことについて説明した。
「偏差」は、平均からの離れ具合であり、個々の個性を表した。
「分散」は、偏差の二乗の和であり、全体として個性の強さ、バラツキを表した。
そして、「標準偏差」は、分散の平方根をとったものであった。

ふむ。ようやく、データを観察する上での第一段階はクリア。
では、異なる変数（属性）同士の関係を観察するには、どうすれば良いか？
まずは、前回のデータの読み込みから。前回と同じ作業。

# Windows と Linux の人は次のコマンド
X <- read.table("clipboard", sep=",", header=TRUE)

# Mac の人は次のコマンド
X <- read.table(pipe("pbpaste") , sep=",", header=TRUE)

以下が、そのデータ。前回のものを再掲しただけ。

データの説明に関しては、平均の話を参照。

Oroshi_Kg, Taka, Naka, Yasu
Tai, 1336, 3150, 1217, 53
Chinu, 226, 630, 513, 105
Sawara, 212, 1575, 1259, 840
Yazu, 4834, 840, 315, 179
Suzuki, 618, 1575, 1146, 525
Akou, 21, 4725, 3129, 1575
Kochi, 28, 2100, 874, 210
Okoze, 28, 5250, 2635, 210
Ainame, 29, 3150, 621, 315
Hage, 1205, 1890, 383, 210
Konoshiro, 21, 2100, 1100, 105
Sayori, 12, 5250, 3916, 1260
Anago, 1637, 2625, 1721, 630
Mebaru, 330, 4200, 1680, 315
Tachiuo, 1211, 2310, 930, 105
Hamo, 1622, 3938, 592, 53
Karei, 41, 6300, 3178, 525
Hirame, 540, 2100, 1496, 210
Aji, 5149, 3150, 789, 210
Saba, 5413, 3675, 625, 53
Ika, 2414, 2888, 1083, 105
Tako, 1844, 1890, 1180, 525
Ebi, 560, 7350, 1420, 525
KurumaEbi, 222, 8925, 6508, 2625
Koiwashi, 1316, 1155, 617, 328
ObaIwashi, 2716, 499, 267, 175
Mentai, 678, 1155, 859, 394
Hamachi, 4772, 998, 632, 105
Isaki, 268, 2625, 1465, 840

現在、Xという行列の変数に、「魚の卸売価格」の情報が入っている。
このデータを使って、各変数間の関係を観察していくことにする。

さて、ある変数と別の変数との関係とは何か？
そもそも、二つの変数の間に、関係が「ある」のか？、「ない」のか？

もしも、関係があるとしたら、何らかの法則があるはず。
直感的に、考え得る法則は、二つある。
一つ目は、一方の値が高いときに、もう一方も高いという法則（連動）、
二つ目は、一方の値が高いときに、もう一方は低いという法則（反連動）。

どうすれば、この状況が理解できるか、少し考えてみる。

まずは、連動している場合を考えてみる。どのような状況であるか？
一方の値が「＋」に向いているとき、もう一方は「＋」に向く。
一方の値が「−」に向いているとき、もう一方は「-」に向く

ここで、小学校で習ったことを、思い出してみる。
「+」 × 「+」 = 「+」　であり　「-」×「-」 =「+」　また　「+」 × 「-」 = 「-」　
という、あの法則。確か、そのようなことを教わった。

・連動の場合、ある対象における二つの変数をかけ合わせると、
　その結果は、「+」方向に大きくなるし、
・反連動の場合、ある対象における二つの変数をかけ合わせると、
　その結果は、「-」方向に大きくなる。
・いずれでもない場合、「+」と「-」が打ち消し合って、「0」に近づく。

では、実際に、どのように計算するのかを整理する。
二つの変数、「」と「」があったとして、次のように考えることができる。

この計算結果の「正負の符号」を見れば、上の法則のいずれかが解る。

ところで、この式を観察してみる。
なんと、分散の式のΣの計算の部分と似ている。

もしも、＝ であれば、二乗となるので、この部分は分散と同じ。
したがって、分散の式に合わせて、次のように考えることもできる。



二つの変数の「バラツキ」を同時に見ているので、これを「共分散」と呼ぶ。
とにかく、この指標の「正」と「負」の大きさを観察すれば、
二つの変数間の全体として関係が解りそう。では、実際の計算を。

そうそう、行列Xと同じサイズの、平均が並んでいる行列が必要なのだった。
これは、前回にもやったので、説明は省く。処理の内容だけ。

# ステップ1：データの総数（n）を求める。

n <- nrow(X)

# ステップ2：卸売数量の平均を求める。転置が必要。

m <- t(colMeans(X))

# ステップ3：行列の引き算ができるように、平均値を総数個複製する。

M <- matrix(rep(m, each=n), nrow=n, ncol=4)
# ステップ4：元の行列から、平均値が複製された行列Mを引き算する。
A <- X-M

とりあえず、これで準備完了。気を取り直して、計算してみる。
全組み合わせは大変。今回は、卸売数量と高値の関係に注目してみる。
この二つの共分散を求める方法は以下の通り。ベクトルの掛け算。

(A$Oroshi_Kg %*% A$Taka)/(n-1)

そして、その結果。

> (A$Oroshi_Kg %*% A$Taka)/(n-1)
         [,1]
[1,] -1026368

ふむ。マイナスの結果が出てきた。反連動の関係のようだ。
っで、これは、大きいのか？それとも、小さいのか？
場合によっては、実は、反連動とは言えないかも。判断に困る。

可視化してみると、状況が解るかもしれない。
とりあえず、元の行列Xで確認してみる。

plot(X$Oroshi_Kg,X$Taka)

これを「散布図」と呼ぶ。横軸が卸売数量で、縦軸が高値。
共分散によって表される状況を可視化したものと言える。

さて、結果は...というと、何とも、微妙な...。
なんとなく、右肩下がりになっているような...。
とにかく、コメントしにくい結果である。正直言って、解釈に困る。

反連動と言えなくはない。が、やはり、その程度は不明。
そこで、この結果を、なんとか一般的にわかりやすいようにしたい。
どのようにすれば、良いか。考えてみる。

要するに、この問題は、共分散の上限と下限が決まっていないのが問題。
したがって、上限と下限が、固定されれば問題は解決できそうである。

わかりやすいように、偏差の値が、「-1〜1」の間に収まるようにすれば良いか。
つまり、何かの単位で揃えるということか。うん？単位にする？
単位にするということは、何かの逆数を掛ければ良かった。

では、二つの変数の共分散を考えた時、
正の方向に値が最も大きくなるのはどような状況か？

ここで、直感的に分かった人はセンスが良い。実は、二つの変数が同じ状況。
この状況は、分散と呼ばれる状況であった。では、分散で単位化すれば良いか？

それはでは「ダメ」！分散は、元の値よりも二乗分値が大きいのであった。
では？何で単位を揃えるべきか。そうそう。思い出した。
平方根をとって元と比較したもの。「標準偏差」を使おう。

まずは、標準偏差を求めたい。これも復習。sd()関数を使う。
今回は、結果を変数に入れておく。カッコで括っているのは、結果を表示させるため。

(X.sd.o <- sd(X$Oroshi_Kg))

(X.sd.t <- sd(X$Taka)) (X.sd.n <-  sd(X$Naka)) (X.sd.y <- sd(X$Yasu))

実行結果は、以下のようになる。

> (X.sd.o <- sd(X$Oroshi_Kg))
[1] 1677.734
> (X.sd.t <- sd(X$Taka))
[1] 2052.045
> (X.sd.n <-  sd(X$Naka))
[1] 1323.899
> (X.sd.y <- sd(X$Yasu))
[1] 553.1923

次に、行列Aを作成した時の「ステップ3」を応用して、
行列Xと同じサイズの標準偏差が入った行列を作成する。今回、結果は省略。

(S <- matrix(rep(c(X.sd.o, X.sd.t, X.sd.n, X.sd.y), each=n), nrow=n, ncol=4))

次に、偏差が入っている行列、すなわち行列Aを標準偏差で基準化する。

行列の「成分の割り算」であって、行列そのものの計算ではない。

間違っても、逆行列と混同してはいけない。

Z <- A/S

さて、ちゃんと基準化されているかを確認してみる。

各変数における分散が1となっていれば良いはず。

var((A/S)$Oroshi_Kg)
var((A/S)$Taka)
var((A/S)$Naka)
var((A/S)$Yasu)

実行結果は以下のようになる。確かに、上手くできている。

> var((A/S)$Oroshi_Kg)
[1] 1
> var((A/S)$Taka)
[1] 1
> var((A/S)$Naka)
[1] 1
> var((A/S)$Yasu)
[1] 1


ついでに、この時の平均についても確認しておく。平均は全て「0」。

colMeans((A/S))

そして、その実行結果。結果は、微細な計算誤差があるが限りなく「0」。

結局、上限が1となるように、上限と下限を圧縮しているので、

ちょうど、真ん中となる平均は、基準化する前と変わらず「0」のまま。

> colMeans((A/S))
    Oroshi_Kg          Taka          Naka          Yasu 
-3.840319e-17  1.007336e-16 -6.346695e-17  1.339924e-17 

以上の結果から、このように基準化を行うと、分散が1、平均が0となる。

このような基準化の方法を「標準化」と呼び、標準化された値を「Z値」と呼ぶ。

なお、標準化の英語は「scale」と呼ぶ。Rには、scale()関数がある。

かなり、面倒な手続きを踏んだが、以下のようにすれば、

元の行列Xから、簡単に標準化された行列を算出できる。

Z <- scale(X)

ここで、標準化の式も示しておく。以下の値は、個々の成分に作用する。

左から順に。左端は、「i」番目の変数「x」の「z」値という意味。
真ん中のは、「i」番目の変数「x」における実際の計算の式。
右端は、象徴的な式。ある変数xから平均を引いて、
標準偏差で基準化している状況を表している。

では、標準化されたデータを用いて共分散を計算してみる。
共分散は、ベクトルの掛け算で計算できた。解は、-1 〜 1 の間に収まるはず。

(Z$Oroshi_Kg %*% Z$Taka)/(n-1)

そして、その結果。

> (Z$Oroshi_Kg %*% Z$Taka)/(n-1)
           [,1]
[1,] -0.2981213

ふむ。この結果を見てみると、負の値となっているが、それほど大きくない。

卸売数量と高値の関係を、反連動というには難がありそう。

さて、このように標準化を行ったデータで共分散を見てみると、

相互の関係性が「-1〜1」の間に収まるので、見通しが明瞭になる。便利である。

こういった便利な方法には、大抵名前がついている。

これが、いわゆる「相関係数」と呼ばれる指標である。

これまでは、「連動」と「反連動」という言葉を使ってきたが、一般的では無い。
通常は、この相関係数から、「正の相関」と「負の相関」と呼ぶ。

余談ではあるが、「負」という言葉な否定的なイメージを持っているため、
「分析上、あまり良くない傾向が出ている」と勘違いする人がいる。
これは大きな間違い。「負の相関」が強く現れるデータにも意味がある。

ところで、今回は、元の値を標準化して、その共分散として相関係数導いた。
しかし、相関係数は別の方法で導くことも可能であり、
以下の式によって、計算することが可能である。
ここでは、二つの変数を「」と「」として、

少し複雑に見えるかもしれない。よく見れば、難しくない。

分子は、共分散のΣの中身になっている。そして、

分母は、「」と「」の分散のΣの中身の平方根になっている。

ここで、「=」の状況を考えてみる。つまり、2つの変数が同じ状況。

確か、平方根の二乗は、元の数に戻るのだったから、

となり、相関係数の値が出てくる。よくできている。

もちろん、Rには、相関係数を算出する関数がある。

通常は、このような回りくどいやり方をせず、次のようにする。

cor(X$Oroshi_Kg, X$Taka)

cor(X$Oroshi_Kg, X$Naka)

cor(X$Oroshi_Kg, X$Yasu)

実行結果は、以下の通り

> cor(X$Oroshi_Kg, X$Taka)

[1] -0.2981213

> cor(X$Oroshi_Kg, X$Naka)

[1] -0.4143816

> cor(X$Oroshi_Kg, X$Yasu)

[1] -0.3617937

ふむふむ。相関係数は、cor()関数の使うのである。

これは、便利。よく使う関数のひとつである。

上記の結果から、卸売数量と高値、中値、安値の関係を見てみると、
全て負の値をとっているので、反連動的な傾向を持っているようである。
一般的に、魚の希少価値が値段に反映すると言われているので、
そのような傾向が相関係数にも現れているのかもしれない。

ただし、これらの値を見てみると、いずれの値も-1よりも0に近く、
したがって、明らかに、「負の相関」があるとは言えない。

希少価値があっても、クセが強くて用途が限られる不人気な魚や、
逆に、希少価値は低いが、需要が高く、結果として高価となる魚もある。
今回のデータでは、そうした例外的な状況が影響しているようである。

最後に、直感的に強い関連がありそうな組み合わせもやってみる。
おそらく、中値と安値は、強い正の関係がありそう。どうなる？

> cor(X$Naka, X$Yasu)
[1] 0.8704575

結果を見てみると、0.87と1に近い値が出てきた。これは、正の相関と言えそう。
では、このような状況は、どのように可視化できるか？前にやったことを思い出す。

たしか、共分散の関係を標準化したのが相関係数だったので、
散布図で可視化すれば、この状況を上手く表現できるだろう。

plot(X$Naka, X$Yasu)

今度は、明らかに右肩上がりの図になっている。このような散布図を書いた場合、
「正の相関」が強いものは、右肩「上がり」の図となり、各々の点は対角線に集まり、
「負の相関」が強いものは、右肩「下がり」の図となり、各々の点は対角線に集まる。
2つの変数間の関係を観察するには、相関係数と散布図を見比べるのが良い。

なお、「0.7以上は相関がある」などと聞くが、その表現は、正しくない。
相関の基準は、誤差の許容範囲と対象数に依存するのであって、
一般的な基準値が存在しているわけではない。この話は、いずれ詳しく述べる。

文系のための「データのバラツキ」（１）

今回、少々、長い。画像やデータが多いということもある。
しかし、今回の話は、最重要。理解できないと、先に進めない。

早速、本題に入る。データの「真ん中」には、様々な定義があった。

一般的なものは、平均、中央値、最頻値、など。

「真ん中」とは何か、を突き詰めることは面白い。

しかし、本ブログの目的は、多次元データ解析。先に進む。

さて、そもそも、データの「真ん中」の話をした理由は、

各々の対象が「真ん中」から離れている度合いを知りたいからであった。

なぜ、そのようなことを知りたいのか？

「離れている度合い」によって、対象の「個性の強さを知る」ためである。

ふむ。解ったような、解らないような。では、簡単な質問をする。

対象全体としての「個性」の「バラツキ」は、大きい方が良いか？

それとも、対象全体としての「バラツキ」は、小さい方が良いか？

全体としての「バラツキ」が、小さいということは、

全体として「個性」が「無く」、観察できる視点は限定される。

したがって、答えは、「バラツキ」が大きい方。

何故、そのようなことが言えるのか、

冷静に考えてみれば、簡単に理解することができる。

例えば、教室Aに出席している100人の学生に「将来なりたい職業は？」と質問し、

100人全員が「地方公務員」と回答したとする。

そして、別の教室Bに出席している別の100人に同じ質問をして、

10人が「大学教授」、20人が「弁護士」、20人が「医者」、

30人が「公認会計士」、25人が「起業家」、残り5人が「宇宙飛行士」、

と答えたとする。個性的なクラスである。熱意と野望と夢に満ちている。

できれば、このような個性的なクラスを受け持ちたい。まぁ、良い。次に進む。

さて、教室Aは、全ての学生が「地方公務員」と答えた。夢の無いクラス。
一方、教室Bは、回答が大きく別れた。個性に溢れた夢のあるクラス。

分析対象として、興味が湧いてくるのは明らかに教室Bである。なぜなら、

• 性別の差は存在するのか？
• 年齢が微妙に異なっているのではないか？
• 出身地が影響のではないか？
• 両親の職業が関係している可能性は？

などなど。バラツキが大きい、すなわち、個性が豊かであるということは

新しい問題を発見する上で重要な手がかりとなる。

それゆえに、データ分析では「バラツキ」が重要となる。

かつて、「ゆとり教育」というものがあった。

あの政策が出てきた理由は「個性豊かたな人間を育む」とか、

たしか、そのような理由であった。評価はともかく、「理念」は間違っていない。

とにかく、分析において「バラツキ」というのは、「個性」の象徴であるので、

「バラツキ」の大きいことは良いことなのである。分析者としては。

では、いよいよ、バラツキについて本格的に考えていく。

とりあえず、平均の話で用いた「魚の卸売価格」のデータを用いて考える。

何度もやってきたので、データの読み込み方は大丈夫であろう。

# Windows と Linux の人は次のコマンド
X <- read.table("clipboard", sep=",", header=TRUE)

# Mac の人は次のコマンド
X <- read.table(pipe("pbpaste") , sep=",", header=TRUE)

以下が、そのデータ。前回のものを再掲しただけ。

データの説明に関しては、平均の話を参照。

Oroshi_Kg, Taka, Naka, Yasu
Tai, 1336, 3150, 1217, 53
Chinu, 226, 630, 513, 105
Sawara, 212, 1575, 1259, 840
Yazu, 4834, 840, 315, 179
Suzuki, 618, 1575, 1146, 525
Akou, 21, 4725, 3129, 1575
Kochi, 28, 2100, 874, 210
Okoze, 28, 5250, 2635, 210
Ainame, 29, 3150, 621, 315
Hage, 1205, 1890, 383, 210
Konoshiro, 21, 2100, 1100, 105
Sayori, 12, 5250, 3916, 1260
Anago, 1637, 2625, 1721, 630
Mebaru, 330, 4200, 1680, 315
Tachiuo, 1211, 2310, 930, 105
Hamo, 1622, 3938, 592, 53
Karei, 41, 6300, 3178, 525
Hirame, 540, 2100, 1496, 210
Aji, 5149, 3150, 789, 210
Saba, 5413, 3675, 625, 53
Ika, 2414, 2888, 1083, 105
Tako, 1844, 1890, 1180, 525
Ebi, 560, 7350, 1420, 525
KurumaEbi, 222, 8925, 6508, 2625
Koiwashi, 1316, 1155, 617, 328
ObaIwashi, 2716, 499, 267, 175
Mentai, 678, 1155, 859, 394
Hamachi, 4772, 998, 632, 105
Isaki, 268, 2625, 1465, 840

では、早速、このデータの「個性」を見てみる。

これまでの説明から、最初に、データの「真ん中」を算出し、

各対象の値からの「離れ具合」を求める。とりあえず、以下を順次実行していく。
今回は、とりあえず「平均」を「真ん中」にして考えることにする。
（実際には、中央値の方が適しているが、その議論はいずれ。）

# ステップ1：データの総数（n）を求める。

n <- nrow(X)

ステップ1は、大丈夫だろう。行数を変数nに代入しているだけ。

# ステップ2：卸売数量の平均を求める。転置が必要。

m <- t(colMeans(X))

colMeans()関数は、行列の列平均を算出してくれた。
便利であるが、数値の部分のみを取り出すには「転置」が必要。
これは、数学的な問題ではなく、Rというコンピュータ言語の仕様の問題。

# ステップ3：行列の引き算ができるように、平均値を総数個複製する。

(M <- matrix(rep(m, each=n), nrow=n, ncol=4))

ここでは、新しくrep()関数が登場している。Replicate（複製する）の略。
あるスカラーの値を特定のルールに従って複製することができる。
この場合は、「each = n」つまり「それぞれの値をn回複製する」という意味。
では、変数nには何が入っていたか？確か、行数（対象数）であった。
ここでは、rep()関数は、matrix()関数の中に入っていて、行数はnとなっている。
したがって、元の行列と同じサイズの、平均値が並んだ行列が作られた。

# ステップ4：元の行列から、平均値が複製された行列Mを引き算する。
(A <- X-M)

最後のステップは簡単。単なる行列の引き算。自信が無い人は復習コース。
少し複雑だったかもしれない。慣れれば、大したことは無い。

さて、最後の結果は、新しく作られた変数Aに格納されている。
この行列Aの各成分は、各対象について、平均が引かれたものになっている。
この指標が、各対象の「個性」の強さを表しているものであり、「偏差」と呼ぶ。

ここで少し数式も入れる。既にやったことを数式で表現しているにすぎない。



ここで、対象の数は 1 < i < n 個あり、変数は 1 < j < p 個ある。いつもの「np行列」。

まずは、記号の読み方。「X」の上に「〜」が書いてある。
読み方は「エックス・チルダ」。各成分について「列平均」を引いた行列を表す。
当然、列の数（p個の列）だけ「平均」あるので、列平均を「」のように表現する。

この行列のことを「偏差行列」と呼ぶ。偏差が並んでいる行列。頻繁に登場する。
元の多次元データから、各対象の「個性」だけの行列に作り替えたものである。

では、偏差行列を用いて何ができるのか？試してみる。

この状況では、少々、観察しにくいので、まずは、可視化してみる。可視化は重要。
今回は、結果を美しく見せるためにちょっとした工夫がしてある。説明は後で。

# 卸売数量の可視化

plot(A$Oroshi_Kg, type = "h", xaxt="n", xlab="", main="Oroshi")

abline(h=0)
axis(side=1, at = (1:n), labels = rownames(A), las=2)

# 高値の可視化
plot(A$Taka, type = "h", xaxt="n", xlab="", main="Takane")
abline(h=0)
axis(side=1, at = (1:n), labels = rownames(A), las=2)

# 中値の可視化

plot(A$Naka, type = "h", xaxt="n", xlab="", main="Nakane")
abline(h=0)
axis(side=1, at = (1:n), labels = rownames(A), las=2)

# 安値の可視化

plot(A$Yasu, type = "h", xaxt="n", xlab="", main="Yasune")
abline(h=0)
axis(side=1, at = (1:n), labels = rownames(A), las=2)

まずは、上記のコマンドの説明から。

基本的に、4つの図は同じ方法で作成されている。

plot()関数を使って、基本的な図を書き、

次に、abline()関数を使ってゼロを通る水平線を書いて、

最後に、axis()関数を使って軸ラベルを書いている。

各関数のパラメータを見てみる。

• plot()関数
• type: プロットの種類。棒で表示するので「h」。「histogram」の「h」
• xaxt: この段階では、x軸を描画しないので「n」。「No」 の「n」
• xlab: 今回は、x軸の軸ラベルを表示しないので空白。
• main: 図の上に表示させるグラフタイトル。
• abline()関数
• h: 水平線を通る線を指定。通る場所は「0」。
• axis()関数
• side: ラベルを配置する場所。下に表示するので「1」を指定。
• at: 対象のラベルを配置する位置。等間隔に「1〜総数個」。
• labels: 実際に配置する文字列。行列Aの列名を使うので「rownames(A)」。
• las: ラベルの向き。軸に対して垂直に置くので「2」。

Rでは、様々なパラメータによって図を整える。

上記以外にも、パラメータはあるが、今回は必要なもののみ。

余裕がある人は、Rのヘルプを参照すると良い。詳しい説明がある。

例えば、axis()関数のヘルプを参照するには、

コンソール上で、「?axis()」のように、前にハテナマークを付けるだけ。

次に、結果から何を読み取れるかを考えてみる。

基本的なこととして、平均は丁度0を通る水平線となっていて、

この線よりも上に出ているのは、平均値よりも高く、

この線よりも下に出ているのは、平均値よりも低いことを示している。

とりあえず、「車海老」を観察してみる。

車海老は、卸売数量が平均を下回り、高値、中値、安値が平均よりも高い。

したがって、希少価値が高くて、全体的に価格も高いと言える。

「鯵」は、卸売数量が平均を上回り、高値、仲値、安値ともに平均以下。

つまり、希少価値は低く、全体的に安いと言える。車海老の逆の傾向

面白いのは「チヌ」。卸売数量は、平均よりも低いが、

高値、中値、安値ともに、平均よりも低い。

マダイなど、他のタイ科よりもクセがあるから、

ニーズが限られるのかもしれない。

このように、平均からの離れ具合を見てやると、

それぞれの対象の個性が見えてくる。この感覚が重要である。

今回のデータから、卸売数量と価格の関係を観察したくなるが、
これは、もう少し先の話。もう少し、基礎的な話で我慢する。

さて、今回は、話が長くなる。もう少しの辛抱。
基本的な話は、時して退屈であり、また、そのような話を長々とするのも辛い。

話が逸れた。今度は、全体としての「個性」というのを考えてみる。
全体としての個性というのは、「偏差」を合わせたものと考えることができる。
すでに、「偏差」は変数Aの行列に格納されていた。
これを足しあわせてみるとどうなるか？実行する前に気づいただろうか？
Rでは、列和を次の方法で計算することができる。

colSums(A)

以下は、実行結果。

> colSums(A)
    Oroshi_Kg          Taka          Naka          Yasu
-1.818989e-12  5.002221e-12 -2.501110e-12  0.000000e+00

あれ〜？何かおかしい。若干の誤差があるが、全てが「0」となっている。
なぜか？この場合は、「平均」をデータの「真ん中」にしている。
したがって、偏差の負の側と偏差の正の側を足し合わせると、
結果として、両方を「平（たいら）に均（な）して」しまい、「0」となる。
この方法は、どうやら良くない。では、どうするべきか？

簡単な方法は、最初から、マイナスにならないようにする。
そのような方法は、これまでにやったことがあるだろうか？
思い出してみる。確か、二乗するという方法があった。これを試してみよう。

注意！これから行うのは、行列の計算ではない。行列の「成分の二乗」。

colSums(A^2)

この計算の結果は以下の通り。各列の「二乗和」が計算されている。

> colSums(A^2)
Oroshi_Kg      Taka      Naka      Yasu
 78814194 117904883  49075875   8568607

コンピュータの世界では、べき乗を「^（ハット）」で表す。Rでも同様。
さて、この状況には一つ問題がある。全体の個性を表すことはできているが、
これでは、新しい対象が加わる毎に増えていき、他の状況とは比較できない。
例えば、今回のデータの場合、別の日の「卸売価格」との比較ができない。

どうすれば良いのだったか？各対象の数に合致するように、
対象数の単位で基準化するのだった。つまり、対象数の逆数を掛ける。

colSums(A^2/n)

この計算結果は以下の通り。偏差の二乗の平均。

> colSums(A^2/n)
Oroshi_Kg      Taka      Naka      Yasu
2717730.8 4065685.6 1692271.6  295469.2

この状況も数式を使って一般化してみる。
この式の見方は、既に、平均の話でしている。省略。



「σ」の記号は「シグマ」と呼ぶ。「Σ」の小文字。
二乗和であるため、二乗が付いている。これ重要。

この式は、「母分散」と呼ばれる「バラツキ」の指標である。
統計学では、「母集団」と「標本」という考え方がある。
母集団というのは対象の全てであり、標本はその一部。
多くの場合、「母集団」を特定することは困難なので、
標本をとっていることになる。ふむ。そいうものであるのか。

確かに、今回のデータは、色々な意味で「全体の一部」であった。
ところで、全ての魚介類を一つの市場で扱うことは可能なことなのか？

まぁ、それは良いとして、なるほど、「標本」データである。
実は、この場合、母分散の推定値にはならないことが統計学的に証明されている。
実際の標本の分散は、母分散よりも少し大きくなる。

では、どうするべきか？そこで登場するのが「不偏分散」。
何やら、小難しい言葉が出てきた。ここは、我慢が重要。
それほど複雑ではない。以下のように式を変更する。


変更された部分は一箇所のみ。nの逆数を掛けていた部分を、n-1の逆数で掛ける。
そのような調整をしている。実際に計算してみる。

colSums(A^2/(n-1))

以下は、実行結果。

> colSums(A^2/(n-1))
Oroshi_Kg      Taka      Naka      Yasu
2814792.6 4210888.7 1752709.8  306021.7

実を言うと、R には不偏分散を計算するためのvar()関数がある。
では、上の結果とvar()関数の結果とを比較する。以下を実行。

var(X$Oroshi_Kg)
var(X$Taka)
var(X$Naka)
var(X$Yasu)

以下は、実行結果。計算誤差があるが確かに同じ。
このことから、var()関数は、不偏分散を求める関数であることが解る。

> var(X$Oroshi_Kg)
[1] 2814793
> var(X$Taka)
[1] 4210889
> var(X$Naka)
[1] 1752710
> var(X$Yasu)
[1] 306021.7

ところで、「分散」というのは、偏差の二乗の和となっている。
したがって、全体の「バラツキ」を考えると、実は、二乗分だけ多い。
これは、少々、不都合なことがある。どのような不都合か？
要するに、元の値や、平均値と比較することができない、という不都合。

元に戻すにはどうすれば良いのか、というのがここでの問題。

すでに実験したように、偏差を足しても「0」となる。
そこで、分散の平方根（ルート）を取ることにする。念のため。

　であり　

平方根というのは、二乗して元に戻るような値のことで、
「√（ルート）」という記号によって表したのだった。
Rでは、sqrt()関数というのが用意されているので、これを使う。

sqrt(var(X$Oroshi_Kg))
sqrt(var(X$Taka))
sqrt(var(X$Naka))
sqrt(var(X$Yasu))

そして、以下が実行結果。

> sqrt(var(X$Oroshi_Kg))
[1] 1677.734
> sqrt(var(X$Taka))
[1] 2052.045
> sqrt(var(X$Naka))
[1] 1323.899
> sqrt(var(X$Yasu))
[1] 553.1923

これが、対象数で単位を揃えた偏差の大きさとなる。
つまり、ここで計算された結果は「標準化された偏差」と言える。
一般的には、これを「標準偏差」と呼び、全体としての個性の指標としている。

実は、標準偏差に関しても、Rには関数が用意されている。sd()関数である。

sd(X$Oroshi_Kg)
sd(X$Taka)
sd(X$Naka)
sd(X$Yasu)

いつものように、実行結果を下に。確かに。不偏分散の平方根になっている。

> sd(X$Oroshi_Kg)
[1] 1677.734
> sd(X$Taka)
[1] 2052.045
> sd(X$Naka)
[1] 1323.899
> sd(X$Yasu)
[1] 553.1923

さて、標準偏差に関しても数式で表してみる。慣れが重要。
不偏分散の平方根を取るのだから、難しくはない。



今日の話は長かった。重要なことは以下の3つ。

• 偏差：平均からの離れている程度
• 分散：偏差の二乗の和。全体としてのバラツキの指標
• 標準偏差：分散の平方根をとって、元のデータに合わせた指標。

この3つのことを理解できていれば次の話が解るはず。