●基本統計量(続編) オーダー(桁)がほぼ同じ一次データ全体を代表する統計値に以下のものがある。 平均値(相加平均) 、中央値(中位数,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 のデータ(x 1,x2, ・・,x 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 (例1) 3番目 8 番目 2, 2, 5, 6, 6, 6, 7, 8, 8, 20 という大きさの順に並べた 10 個のデータ群からでは、 5番目 6番目 全体の 1/4 の位置の 1/4 偏差 Q1 は、0.5+10×1/4=3番目のデータだから、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)式において、変数xを(2.2)式にしたがって変数zに変換 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 図1 標準正規分布関数の形 . . . . . . . . . . . . . (2.3) どのような正規分布も(2.2)の無次元化、標準化の変数変換を施すことで、同一の(2.3)式の分布と 2p なる。図1は、平均μ=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 z 0 -1 -0.5 0 0.5 1 -0.5 標準正規分布に従うことを意味する。 実際の m の値 11 (1)母標準偏差が既知であれば、上記2式より母平均の推定値の信頼区間を推定できる。 -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%信頼区間が求まる t分布(自由度3) (2)母標準偏差が未知の場合 0.3 (標本平均値の分布の標準偏差) 0.2 =(標本標準偏差)/√(標本サイズ)として、 (a)標本サイズが大きい(約30以上)のとき、 0.1 標本標準偏差を母標準偏差とみなせる(正規分布と同じ)。 0 (b)標本サイズが小さいとき -3 -2 -1 0 1 2 平均値の分布がt分布するものとして信頼度を推定する。 これは、統計量: t = 3 t (付録参照) m-m が自由度φ=n-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 という1組の数値で与えられる標本抽出を行うことは、母集団と同じ確率的な構 造をもった、互いに独立な確率変数 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 ) ] 第1項 å ( Xi - m ) の期待値 å ( Xi - m ) 2 2 2 2 2 1 2 2 n = s 2 + s 2 + + s 2 = Ns 2 第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 参照) 1. R の初期画面 最後にプロンプト>が出てくる。 文字化けのときは「編集」から「GUI プリファレンス」を選択して Font を MSGothic 等へ。 2.関数電卓としての使い方 平方根演算(黄金比の値) > (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):( 底 が 2 の 対 数 ), 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 の符号) 3.変数と代入 (1)変数は英字が最初にくる > x <- sqrt(2) > x^2 [1] 2 <- (2)大文字、小文字は区別する (3)数字も変数の一部に使える の 2 文字で右辺の値を左辺の変数に代入する意味を持つ 4.ベクトルの生成 > x <- c(1,2,3,4,5) > sqrt(x) [1] 1.000000 1.414214 1.732051 2.000000 2.236068 cと( )で数字を囲む ベクトルに対して関数を適用 するとベクトルが返される ベクトルを効果的に適用できる関数群 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 番目の要素を0にする (インデックスは 1 から 始まる) 5.ベクトルの演算 > 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 6.関数定義 継続を表す記号 (自動的に表示) 関数名 引数 > 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 ) # xが 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≦x≦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 7.行列計算 R で行列を作るには、次の2段階で作成 (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列成分を表示 > m[1,2] [1] 3 1,2 行2列成分を表示 > m[c(1,2),2] [1] 3 4 1行成分を全部表示 > m[1,] 2列成分を全部表示 [1] 1 3 5 > m[,2] 2 行目を非表示(負数で指定) [1] 3 4 他の行の 1,3 列を表示 > m[-2,c(T,F,T)] [1] 1 5 2列目を非表示 > 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) 3つの固有値 $values [1] 13.985686 -1.169156 0.183470 3つの固有ベクトル $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
© Copyright 2024