2番目の資料

●基本統計量（続編）
オーダー（桁）がほぼ同じ一次データ全体を代表する統計値に以下のものがある。
平均値（相加平均）
、中央値(中位数，median)、最頻値(mode)
1
1 N
( x1 + x 2 +  x N ) = å xi
N
N i =1
ただし、 xi ：i 番目のデータ、 x ：平均、 N ：全データ数
単独データの平均： x =
x=
度数分布で与えられたときの平均(加重平均)
但し、 f i ：度数、 n ：階級数、 xi ：階級値
1
(
N x1
f
1
+ x2
f
2
+  xn
f
n
)=
1
N
n
åx f
i
i =1
i
n
N = å f i ：全データ数、
i =1
度数でなく、確率密度関数 f ( x ) で表されている場合は、平均値：x =
ò x × f ( x )dx
データの広がり具合を表わす数値、散布度として、次のものがある。
範囲(range)：サイズ N のデータ（ｘ 1，ｘ2，
・・，ｘ N）において最大値と最小値の差
分散(variance)：
s
2
=
2
2
2
1
1
{( x1- x ) + ( x 2 - x ) +  ( x N - x ) } =
N
N
度数で与えられたとき： 2 = 1
s N
の分散
{(x1- x )
標準偏差(standard deviation)： s = s
平均偏差(mean deviation)：
2
2
f
+  + (xn - x )
2
1
f
N
å ( xi - x )
2
i =1
}= N1 å (xi - x )
n
n
2
i =1
f
但し、N =
i
N
å| x
i =1
i
-x|
変動係数(coefficient of variation)＝標準偏差／平均値
四分位数(4-quantile)：分布の左から累積相対度数 p が p=25%、50%、75%になる値で、
それぞれ第一、第二、第三四分位数と呼ばれる。第二四分位数は中央値(中位数)と同じ。
百分位数はパーセンタイル(percentile)と呼ばれる。
四分偏差(quantile deviation)：Q＝（(第三四分位数)－(第一四分位数)）／2
（例１）
３番目
8 番目
2, 2, 5, 6, 6, 6, 7, 8, 8, 20 という大きさの順に並べた 10 個のデータ群からでは、
５番目
６番目
全体の 1/4 の位置の 1/4 偏差 Q1 は、0.5＋10×1/4＝３番目のデータだから、Q1=5 となる。
全体の 3/4 の位置の 3/4 偏差 Q3 は、0.5＋10×3/4＝8 番目のデータの値だから、Q3=8 である。
したがって、四分偏差 Q=(8－5)/2＝1.5 となる。
10
åf
i =1
データの次元に戻したときのばらつき度合い
1
1
(| x1 - x | + | x2 - x | +  + | x N - x |) =
N
N
n
i
●正規分布関数
大数の法則
期待値がμであるような確率変数 xi (i = 1, 2, , n) の n 個の平均値 xn は、n を無限に
大きくしていくと限りなくμに近づく。すなわち、 lim xn = m である。
n® ¥
言い換えると、母集団からランダムに（無作為に）取り出した n 個のサンプル（標本）の
相加平均（標本平均）は、サンプルの数を大きくすると真の平均に近づく。
中心極限定理
母集団の真の平均と、標本平均と誤差は、多くの場合、分布の形に依存しなくて、近似的
に正規分布に従う。定量的には、
平均μ、標準偏差がσの母集団から n 個の標本による標本平均 xn のとる分布は、平均が
μ、標準偏差が
s
の正規分布になる。
n
0.4
正規分布関数の分布密度関数 f (x) の理論式は(2.1)式となる。
f ( x) =
0.3
( x - m )2
exp( )
2s 2
2p s
1
. . . . . . . .. . . (2.1)
(2.1)式において、変数ｘを(2.2)式にしたがって変数ｚに変換
0.2
すると、(2.3)式になる。
z=
x-m
s
j ( z) =
0.1
. . . . . . . . . . . . . . . . . . . . . . . (2.2)
1
exp(-
0
-4
z2
)
2
-3
-2
-1
0
1
2
4σ
3
図１標準正規分布関数の形
. . . . . . . . . . . . . (2.3)
どのような正規分布も(2.2)の無次元化、標準化の変数変換を施すことで、同一の(2.3)式の分布と
2p
なる。図１は、平均μ＝0、標準偏差σ＝1の正規分布関数(標準正規分布関数)の(2.3)式のグラフ
である。この正規分布の場合、±σの範囲に存在する確率は68.3%、±1.64σで90.0%、±1.96σ
で95%、±2σで95.4%、±2.58σで99%、±3.0σの範囲には、99.7%が含まれる。・・・(2.4)
●母平均の推定
正規分布をする母集団全てのデータを取り出すことは一般に困難なので、その中から無作為に
抽出した少数の n 個のデータ xi（標本(サンプル)という）を使って、母集団の母平均を推定する
ことができる。
標本平均 m の理論分布
n
1 n
1
標本平均： m = å xi ，標本標準偏差： s =
( xi - m) 2
å
n i =1
(n - 1) i =1
．
．．(2.5)
2
1.5
このとき、標本平均値の分布は、再び正規分布する性質がある。
（標本平均値の分布の平均値 m）＝（母平均μ）すなわち、 m = m
1
（標本平均値の分布の標準偏差 s）＝（母標準偏差σ）／√(標本サイズ)
すなわち、統計量： z = m - m
s/ n
．．．(2.6) が平均値 0、標準偏差 1 の
0.5
ｚ
0
-1
-0.5
0
0.5
1
-0.5
標準正規分布に従うことを意味する。
実際の m の値
11
（１）母標準偏差が既知であれば、上記２式より母平均の推定値の信頼区間を推定できる。
－z0＜z＜z0 の範囲、すなわちμ－z0･ s ＜ m ＜μ＋z0・ s 、または変形させて
n
n
m －z0・ s ＜μ＜ m ＋z0･ s の範囲に母平均が含まれる確率が正規分布から求まる。
n
n
例えば z0＝1 の場合は 68.3%となり、z0＝2 の場合は 95.4%となる。これが信頼区間となる。
母平均の信頼区間：
m - z0 ×
s
s
< m < m + z0 ×
n
n
．
．．(2.7)
正規分布
0.4
(例) z0=1.64 の場合は 90%信頼区間が求まる
ｔ分布(自由度３)
（２）母標準偏差が未知の場合
0.3
（標本平均値の分布の標準偏差）
0.2
＝（標本標準偏差）／√(標本サイズ)として、
（ａ）標本サイズが大きい(約３０以上)のとき、
0.1
標本標準偏差を母標準偏差とみなせる（正規分布と同じ）。
0
（ｂ）標本サイズが小さいとき
-3
-2
-1
0
1
2
平均値の分布がｔ分布するものとして信頼度を推定する。
これは、統計量： t =
3
t
(付録参照）
m-m
が自由度φ=ｎ－1 の t 分布に従うから。
s/ n
母平均の信頼区間は、 m - t ×
s
s
< m < m+t×
n
n
となる。・・・(2.8)
●分布と標本
標本と母集団の統計量
母平均(population mean) m 、母分散(population variance) s の正規分布をする大きな母集団か
2
ら、x1 , x2 ,, x N という１組の数値で与えられる標本抽出を行うことは、母集団と同じ確率的な構
造をもった、互いに独立な確率変数 X i がこれらの値を持ったと考える。
標本平均(sample mean)：m =
1
N
N
åX
i =1
i
、標本分散(sample variance)：V =
1 N
( X i - m)2 とする。
å
N i =1
標本平均の期待値： [ m] =
1
1
{[ X 1 ] + [ X 2 ] + [ X N ]} = × Nm = m ・・(2.9) ([ ]は期待値を表すとする)
N
N
標本分散の期待値： [V ] =
(付録参照）
1 N
N -1 2
[( X i -m) 2 ] =
s < s 2 となり母分散より小さく偏る。
å
N i =1
N
・・(2.10)
そこで、 s 2 =
N
1 N
2
2
V=
å ( X i - m)2 なる変数を考えると、その期待値は、 [ s ] = s となり
N -1
N - 1 i =1
母分散に等しく偏りが無くなる。この s2 を不偏分散(unbiased variance)とよぶ。
また、標本平均自身の、その期待値からのばらつきの程度を表わす分散の期待値は、
[(m - [m]) 2 ] ＝ [(m - m ) 2 ] ＝
s2
となる。
N
12
付録
t 分布の分布密度関数： f ( x, n ) =
æ
1
x2 ö
çç1 + ÷÷
æn 1ö
nø
n × Bç , ÷ è
2
2
è
ø
-
n +1
2
x2
，
lim f ( x, n ) =
n® ¥
1 -2
e
2p
自由度 n が無限大になると、標準正規分布関数となる。
20
ここで、B(a,b)はベータ関数と呼ばれるもので、
B (a , b ) =
1
G(a ) × G (b)
= ò t a -1 (1 - t )b -1 dt
0
G ( a + b)
で定義される。
Beta 関数の例を Mathematica でグラフ表示すると
右図のようになる。（下はコマンド例）
10
-4
-2
2
4
ベータ関数
-10
Plot[Beta[x,0.5],{x,-4,10},PlotStyle->Thickness[0.01]];
-20
またΓ(x)は、ガンマ関数と呼ばれるもので、
¥
G ( x) = ò t
0
G ( n) =
15
x -1 - t
e dt
で定義される。
10
5
G(n + 1)
となり、Γ(1)＝1 であり、
n
-4
n が正の整数のとき、Γ(n+1)＝n !（n 階乗）となる。
-2
2
4
-5
-10
Plot[Gamma[x],{x,-4,6},PlotStyle->Thickness[0.01]];
13
-15
ガンマ関数
6
付録
標本分散の期待値： [V ] =
N -1 2
1 N
s の証明
[( X i -m) 2 ] =
å
N i =1
N
å ( Xi - m)
= å {( Xi - m ) + (m - m )}
= å ( Xi - m ) + 2å ( Xi - m )(m - m ) + å (m - m )
= å ( Xi - m ) - 2(m - m )å ( Xi - m ) + N (m - m )
= å ( Xi - m ) - 2(m - m )( Nm - Nm ) + N (m - m )
= å ( Xi - m ) - 2 N (m - m ) + N (m - m )
= å ( Xi - m ) - N (m - m )
2
2
2
2
2
2
2
2
2
2
2
2
2
[
] は，
[å ( Xi - m ) ] = [(x - m ) ]+ [(x - m ) ]+  + [(x - m ) ]
第１項 å ( Xi - m ) の期待値 å ( Xi - m )
2
2
2
2
2
1
2
2
n
= s 2 + s 2 + + s 2
= Ns 2
第２項 N (m - m )2 の期待値は，
[
]
[
N (m - m ) = N ´ (m - m )
2
=N´
2
]
s2
（これは正規分布の性質）
N
=s2
従って
[
]
å ( Xi - m) = Ns 2 - s 2 = (N - 1)s 2 であり
結局
2
[V ] = éê 1 å ( Xi - m )2 ùú = N - 1s 2
ëN
û
N
14
【演習 9】ある母集団から任意抽出したサイズ 10 のデータが次のとおりであった。
27.0 22.3 23.5 26.2 24.5 26.0 20.9 25.5 24.4 24.7
白い面積が
全体の 95%
（母集団の標準偏差が判明しているので、正規分布しているものとして(2.7)式を使う）
(1) 母標準偏差がσ＝√3.40 と判明したとき、母平均μの 95%信頼区間を求めよう。
95%信頼区間は 95%の確率で母平均が存在する
範囲（右図の正規分布の両端の 2.5%ずつの範囲
0.4
を除いた領域）である。このときの z0 の値を
0.3
求めるには、Excel の関数を使う。
この面積が
全体の 2.5%
この面積が
全体の 2.5%
0.2
Excel では左の 2.5%の灰色面積部分の－z0 の値
0.1
は、－z0=NORMSINV(0.025)で求まる。
0
－z0
-4 -3 -2 -1
z
z0
0
1
2
3
4
(2) 母標準偏差が不明の場合の、母平均の 95%信頼区間を求めよう。
（母標準偏差が不明で、データサイズ n が小さいので、t 分布に従うとして(2.8)式を使う）
(2.8)式における t 値を求めるには、Excel の関数 t=TINV(1－0.95, n－1)を使う
【演習 10】さいころを振ったときに出る目の数の分布は一様分布となる。さいころを 10 回、
100 回、1000 回振ったときの模擬データを Excel 関数 RANDBETWEEN(1,6)を繰り返し
使って求め、それらの異なるデータサイズの場合の母平均の 95%信頼区間を求めてみよう。
【演習 11】中心極限定理では、ほとんどの確率分布関数において、平均μ、標準偏差がσの母集
団から n 個の標本による標本平均 xn のとる分布は、平均がμ、標準偏差が
るという。これを確かめるために、いま
x=1,2,3,･･･,50 のときの分布関数が
0.06
一様分布：f(x)=0.02、
0.05
指数関数分布：g(x)=exp(－0.05*x)/17.903
0.04
三角関数：h(x)=((1－cos(x*3.1416/25))/50
0.03
のときの値を求め、50 個のデータからサンプルとして、
0.02
ランダムに 3 個、10 個、30 個を取り出したときの
0.01
平均値の 95%信頼区間を求めてみよう。
s
の正規分布にな
n
一様
指数
三角
15
49
45
41
37
33
29
25
21
17
9
13
5
1
0
●データ解析環境 R の基本技（http://cse.naro.affrc.go.jp/takezawa/r-tips/r2.html 参照）
１． R の初期画面
最後にプロンプト＞が出てくる。
文字化けのときは「編集」から「GUI プリファレンス」を選択して Font を MSGothic 等へ。
２．関数電卓としての使い方
平方根演算（黄金比の値）
> (1+sqrt(5))/2
[1] 1.618034
> (sqrt(123))^2
[1] 123
> exp(log(2)+log(3))
[1] 6
> sin(3)^2+cos(3)^2
[1] 1
>
累乗演算
指数・対数
三角関数
演算子＋、－、＊、／、＾または＊＊(累乗)、％％(剰余)、％／％(整数商)
関数群
sin(x), cos(x), tan(x),
sinh(x), cosh(x), tanh(x),
asin(x), acos(x), atan(x),
asinh(x),
acosh(x), atanh(x), log(x):( 自然対数 ), log10(x):( 常用対数 ), log2(x):( 底が２の対数 ),
log1P(x):(log(1+x)のこと), exp(x), expm1(x):(exp(x)-1 のこと), sqrt(x), round(x):(四捨五入),
floor(x):(小数切捨), ceiling(x):(小数切上), trunc(x):(整数部分), sign(x):(x の符号)
３．変数と代入
(1)変数は英字が最初にくる
> x <- sqrt(2)
> x^2
[1] 2
＜－
(2)大文字、小文字は区別する
(3)数字も変数の一部に使える
の 2 文字で右辺の値を左辺の変数に代入する意味を持つ
４．ベクトルの生成
> x <- c(1,2,3,4,5)
> sqrt(x)
[1] 1.000000 1.414214 1.732051 2.000000 2.236068
ｃと（
）で数字を囲む
ベクトルに対して関数を適用
するとベクトルが返される
ベクトルを効果的に適用できる関数群
sum( ), mean( ), median( ), var( ), cor( ), max( ), min( ), range( ), cumsum( ):(累積和),
prod( ):(総積), diff( ):(前進差分), rev( ):(要素を逆順), pmax( ):(列最大値), pmin( ), sort( ),
rank( ):(各要素の順位), order( ):(各要素の元の位置)
> a <- c(1,3,5,7)
> b <- c(2,4,3,7)
> pmax(a,b)
[1] 2 4 5 7
> pmin(a,b)
[1] 1 3 3 7
> x <- c(1,2,3,4,5)
> x <- x^2
>x
[1] 1 4 9 16 25
> mean(x)
[1] 11
>x
[1] 1 4 9 16 25
> x[2] <- 0
>x
[1] 1 0 9 16 25
16
2 番目の要素を０にする
(インデックスは 1 から
始まる)
５．ベクトルの演算
> c(1,2,3)+c(4,5,6)
[1] 5 7 9
> c(1,2,3)*c(4,5,6)
[1] 4 10 18
> 1:5
[1] 1 2 3 4 5
> 3:-3
[1] 3 2 1 0 -1 -2 -3
> x <- 1:5
>x
[1] 1 2 3 4 5
> p <- c(1,2,3)
> q <- c(4,3,2)
内積演算
> p%*%q
[,1]
[1,]
16
ﾃﾝｿﾙ積演算
> p%o%q
[,1] [,2] [,3]
[1,]
4
3
2
[2,]
8
6
4
[3,]
12
9
6
>x
[1] 1 2 3 4 5
> sum(x)
[1] 15
> mean(x)
[1] 3
> var(x)
[1] 2.5
> range(x)
[1] 1 5
６．関数定義
継続を表す記号
(自動的に表示）
関数名
引数
> factorial <- function(n){
+ if( n < 1 ) return(1)
+ n*factorial(n-1)
返戻値を指定
+ }
> factorial(3)
再帰呼出
[1] 6
> factorial(5)
実行
[1] 120
条件分岐
if( 条件式 ) {
条件式が TRUE のときの実行式
｝else {
条件式が FALSE のときの実行式
条件指定の使用例
if( x == 1.0 )
# x が 1.0 と等しいならば
if( x != 2.0 )
# x が 2.0 でないならば
if( x >= 3.0 )
# x が 3.0 以上ならば
if( x < 4.0 )
# ｘが 4.0 未満ならば
｝
繰り返し文：for
> myfactorial <- function(n) {
+ i <- 1
+ for( j in 1:n ){
+
i <- i*j
+ }
+ return(i)
+}
> myfactorial(4)
[1] 24
> myfactorial(5)
[1] 120
論理演算子の使用例
if( (x <-1) || (x>1) ) # x<-1 または x>1 ならば
if( (x>=-1) && (x<=1) ) #
if( !(a>b) )
17
-1≦ｘ≦1 ならば
# a>b でなければ（a≦b ならば）
> matrix(1:6,2,3)
[,1] [,2] [,3]
[1,]
1
3
5
[2,]
2
4
6
> matrix(1:6,2,3,byrow=T)
[,1] [,2] [,3]
[1,]
1
2
3
[2,]
4
5
6
７．行列計算
R で行列を作るには、次の２段階で作成
(1)行列の要素をベクトルで用意する。
(2)関数 matrix( )でベクトルから行列にする。
・指定された要素は左の列から順に上から下へ
埋められる。
・引数 byrow=T を指定すると、要素を上の行から順に左から右に埋める（上図参照）
> (m <- matrix(1:6,2,3))
[,1] [,2] [,3]
[1,]
1
3
5
[2,]
2
4
6
行列ｍを生成、全体
に( )をつけると表示
１行２列成分を表示
> m[1,2]
[1] 3
1,2 行２列成分を表示
> m[c(1,2),2]
[1] 3 4
１行成分を全部表示
> m[1,]
２列成分を全部表示
[1] 1 3 5
> m[,2]
2 行目を非表示(負数で指定)
[1] 3 4
他の行の 1,3 列を表示
> m[-2,c(T,F,T)]
[1] 1 5
２列目を非表示
> m[,-2]
[,1] [,2]
[1,]
1
5
[2,]
2
6
行(row)ﾍﾞｸﾄﾙを結
合して行列生成
> (x <- rbind(c(1,2,3),c(2,4,6)))
[,1] [,2] [,3]
列 (colomn) ベクトル
[1,]
1
2
3
[2,]
2
4
6
を結合して行列生成
> (y <- cbind(c(1,3,5),c(7,9,3)))
[,1] [,2]
[1,]
1
7
[2,]
3
9
[3,]
5
3
18
> (a <- matrix(1:4,2,2)) # 行列生成表示
[,1] [,2]
[1,]
1
3
[2,]
2
4
> (b <- matrix(-2:1,2,2)) # 行列生成表示
[,1] [,2]
[1,]
-2
0
[2,]
-1
1
> a+b
# 行列の加算
[,1] [,2]
[1,]
-1
3
[2,]
1
5
> a%*%b
# 行列の積
[,1] [,2]
[1,]
-5
3
[2,]
-8
4
> a*b
# 各要素毎の積
[,1] [,2]
[1,]
-2
0
[2,]
-2
4
> diag(3)
# 単位行列の作り方
[,1] [,2] [,3]
[1,]
1
0
0
[2,]
0
1
0
[3,]
0
0
1
> diag(1,3)
[,1] [,2] [,3]
[1,]
1
0
0
[2,]
0
1
0
[3,]
0
0
1
> diag(rep(1,3))
[,1] [,2] [,3]
[1,]
1
0
0
[2,]
0
1
0
[3,]
0
0
1
> diag(2,3)
# 対角行列
[,1] [,2] [,3]
[1,]
2
0
0
[2,]
0
2
0
[3,]
0
0
2
三角行列の作り方
> (x <- matrix(1:9,3))
[,1] [,2] [,3]
[1,]
1
4
7
[2,]
2
5
8
[3,]
3
6
9
> x[upper.tri(x)] <- 0
>x
[,1] [,2] [,3]
[1,]
1
0
0
[2,]
2
5
0
[3,]
3
6
9
対称行列の作り方
三角行列を作る
>x
[,1] [,2] [,3]
[1,]
1
0
0
[2,]
2
5
0
[3,]
3
6
9
> y <- x+t(x)
> diag(y) <- diag(y)/2
>y
[,1] [,2] [,3]
[1,]
1
2
3
[2,]
2
5
6
[3,]
3
6
9
三角行列を作るに
はまず普通の行列
を作成
上三角成分を 0 と
することで三角
行列を作る
連立一次方程式の解は、関数 solve( )を使って
解くことができる
対角成分を半分に
する
array で行列作成
> (a <- array(c(1,1,-1,-1,0,1,1,1,0),c(3,3)))
[,1] [,2] [,3]
[1,]
1
-1
1
[2,]
1
0
1
[3,]
-1
1
0
> (b <- solve(a))
solve( )で逆行
[,1] [,2] [,3]
列を求める
[1,]
-1
1 -1
[2,]
-1
1
0
[3,]
1
0
1
> a <- matrix(c(0,1,2,3,4,5,6,7,4),3,3)
> b <- matrix(c(1,0,-2))
> solve(a,b)
3y + 6z = 1
[,1]
x
+
4y + 7z = 0
[1,] -1.0833333
2x + 5y + 4z = -2
[2,] -0.1666667
を ax=b で表現
[3,] 0.2500000
一般逆行列
逆行列の求め方
転置行列と元の
三角行列を加える
パッケージの呼び出し
> A <- matrix(c(1,2,3,4,5,6,7,8,9),3,3)
> B <- solve(A)
以下にエラーsolve.default(A) : Lapack routine dgesv: 線形方程式系は正確に特異です
> library(MASS)
> B <- ginv(A) # Moore and Penrose generalized inverse
>B
一般逆行列
> (a<-matrix(c(1,-1,4,2,0,-2,1,3,0),3,3))
[,1]
[,2]
[,3]
[,1] [,2] [,3]
[1,]
1
2
1
[1,] -0.6388889 -5.555556e-02 0.5277778
[2,]
-1
0
3
[2,] -0.1666667 -9.234353e-17 0.1666667
[3,]
4 -2
0
[3,] 0.3055556 5.555556e-02 -0.1944444
> (b<-c(2,3,4))
[1] 2 3 4
> a%*%b
行列×ベクトル
[,1]
[1,]
12
[2,]
10
ベクトル×行列
[3,]
2
> b%*%a
[,1] [,2] [,3]
[1,]
15 -4 11
> b%*%b
ベクトルの絶対値の平方
[,1]
[1,]
29
固有値と固有ベクトル
> a <- matrix(c(0,1,2,3,4,5,6,7,9),3,3)
> eigen(a)
３つの固有値
$values
[1] 13.985686 -1.169156 0.183470
３つの固有ベクトル
$vectors
[,1]
[,2]
[,3]
[1,] -0.4263524 -0.9366202 0.2264150
[2,] -0.5474684 -0.2042614 -0.8684128
[3,] -0.7200708 0.2846399 0.4411298
19

Download Report