確率・統計

(7) 標本の抽出と要約

今までは、おもに「確率論(Probability Theory)」の話題が中心でした。事象の集合から確率変数に置き換えて、確率密度関数として分布の形状を見ることによって、物事の起こりやすさを定量化することが、確率論の主な目的と考えることができます。元々、ギャンブルに有利に勝つ方法を数学を使って見つけるための手段として確率というものが誕生したわけですが、今では様々な分野で応用されている重要な学問の一つとなっています。

ここからは「統計学(Statistics)」の話題が中心となります。統計学は、たくさんのデータから意味のある情報を抽出するための手法を研究する学問です。そのための道具として確率論が取り入れられていることから、確率と統計はペアで紹介されることが多いわけです。まずはデータを集める手法から始めて、そのデータの性質を調べるための方法について紹介したいと思います。

(注) 数式などの記法について (ドキュメントの中で使用している数式の表現方法に関する注意点です)

1) 標本の抽出

対象となる標本が非常に大きい場合、全てのデータを抽出するのは非常にコストがかかり、実質不可能な場合もあります。平均と分散が存在すればどんな分布でも、標本の抽出を繰り返して標本平均を計算するとそれは正規分布に従い、抽出する数が多いほど分散は小さくなっていきます。よって、母集団から部分的に標本を抽出してデータ解析に利用する場合、標本の数は多い方が信頼性の高い結果が得られるわけですが、たくさんのデータを集めようとするほどコストは高くなっていきます。そこで同じ数のデータを抽出した時、できるだけ信頼性の高いような抽出方法を検討することは非常に重要になります。主な抽出方法としては次のようなものがあります。

単純無作為抽出法 (Simple Random Sample ; SRS)

いわゆる、完全にランダムな選択法です。N 個の母集団から標本を算出する場合、それぞれのデータに 1 から N までの範囲の番号を付けて、乱数によって抽出する番号を決める方法が「単純無作為抽出法(Simple Random Sample ; SRS)」です。

系統抽出法 (Systematic Sampling)

1 から N までデータに番号が付けられて並んでいた時、最初の一つだけをランダムに決めて、あとのデータは一定間隔で抽出する方法です。例えば 100 個のデータから 20 個のデータを抽出する場合、100 / 20 = 5 なので、最初のデータをランダムに決めて 23 だったとしたら、その後は 28, 33, ... と一定間隔で抽出します。データの末尾まで達したら、もう一度ランダムな値を決め直したり、データ数を引いて先頭側から抽出を続けるなどの方法が考えられます。

多段抽出法 (Multistage Sampling)

標本を抽出する対象が大きすぎて、単純無作為抽出法や系統抽出法ではコストが大きくなるような場合、抽出する標本の範囲を段階的に分ける場合があります。例えば、日本国内に住む人を対象にデータを抽出する場合、全ての人に番号を付けて単純無作為抽出法や系統抽出法を利用するのは大変なので、例えば最初に県を抽出し、次に市町村を抽出して、最後に個人を抽出するようなやり方をすることでコストを抑えることができます。このやり方を「多段抽出法(Multistage Sampling)」といいます。

層化抽出法 (Stratified Sampling)

母集団がいくつかのグループから構成されていた時、それぞれのグループに属している確率に応じて抽出する数を決める方法を「層化抽出法(Stratified Sampling)」といいます。例えば、二つのグループがあって、それぞれの要素数の比率が 7 : 3 であれば、抽出する標本の比率も 7 : 3 とします。

これらの方法は組み合わせて利用するものもあります。例えば、多段抽出法で抽出するグループが決まったら、その中からは単純無作為抽出法や層化抽出法を利用するような組み合わせが考えられます。しかし、母集団がいくつかのグループに分かれていた場合、全くランダムに標本を抽出する単純無作為抽出法に対し、層化抽出法ではグループによって抽出する比率を調節しているので、得られる結果も異なることが予想されます。どちらの方が信頼性の高い結果が得られるのでしょうか。

ここで、「信頼性が高い」というのは「分散が小さい」ということを意味します。簡単な例として、二つのグループからなるデータがあって、それぞれのグループは r1 : r2 の比率 ( 但し r1 + r2 = 1 ) であるとします。第 i グループの確率密度関数を pi(x)、平均と分散をそれぞれ μi, σi2 ( i = 1, 2 ) として、単純無作為抽出法と層化抽出法で N 個の標本を抽出した時の平均と分散を求めてみます。

まず、二つの分布を合わせた全体の母集団に対する確率密度関数を P(x) としたとき、ランダムに抽出したデータが第 i グループに属する確率が riであり、それが x である時の条件付確率が pi(x) となるので、

P(x) = r1p1(x) + r2p2(x)

になります。この母集団の平均 μ は

μ=∫{-∞→∞} x P(x) dx
=∫{-∞→∞} x{ r1p1(x) + r2p2(x) } dx
=r1・∫{-∞→∞} x・p1(x) dx + r2・∫{-∞→∞} x・p2(x) dx
=r1μ1 + r2μ2

であり、分散 σ2 は、

σ2=∫{-∞→∞} ( x - μ )2 P(x) dx
=∫{-∞→∞} ( x - μ )2{ r1p1(x) + r2p2(x) } dx
=Σi{1→2}( ri・∫{-∞→∞} ( x - μ )2pi(x) dx )
=Σi{1→2}( ri・∫{-∞→∞} { ( x - μi ) + ( μi - μ ) }2pi(x) dx )
=Σi{1→2}( ri{ ∫{-∞→∞} ( x - μi )2pi(x) dx
 + 2( μi - μ )∫{-∞→∞} ( x - μi ) pi(x) dx + ( μi - μ )2∫{-∞→∞} pi(x) dx } )

展開した各積分値のうち第一項は各グループの分散 σi2 を表し、第二項は

∫{-∞→∞} ( x - μi ) pi(x) dx=∫{-∞→∞} x pi(x) dx - μi∫{-∞→∞} pi(x) dx
=μi - μi = 0

また、第三項の積分値は 1 なので、

σ2 = Σi{1→2}( ri{ σi2 + ( μi - μ )2 } )

になります。ここで、

μ1 - μ=μ1 - ( r1μ1 + r2μ2 )
=( 1 - r11 - r2μ2
=r2μ1 - r2μ2
=r2( μ1 - μ2 )
μ2 - μ=r1( μ2 - μ1 )

となるので、

σ2=( r1σ12 + r2σ22 ) + r1{ r2( μ1 - μ2 ) }2 + r2{ r1( μ2 - μ1 ) }2
=( r1σ12 + r2σ22 ) + r1r2( r1 + r2 )( μ1 - μ2 )2
=( r1σ12 + r2σ22 ) + r1r2( μ1 - μ2 )2

が得られます。よって、母集団の平均は r1μ1 + r2μ2、分散は ( r1σ12 + r2σ22 ) + r1r2( μ1 - μ2 )2 になり、この母集団から N 個の標本を抽出した時の標本平均と標本分散はそれぞれ

m = r1μ1 + r2μ2

s2 = { ( r1σ12 + r2σ22 ) + r1r2( μ1 - μ2 )2 } / N

になります (「(5) 正規分布」の「4) 標本平均と標本分散」参照)。これはちょうど、単純無作為抽出法で抽出した標本に対する平均と分散を意味します。

層化抽出法の場合、それぞれのグループから比率に従って抽出することになるので、N 個の標本を抽出する場合、第 i グループからは riN 個の標本を抽出することになります。確率変数 x1, x2, ... xN がそれぞれ平均 μi, 分散 σi2 の分布に従うとき、

y = a0 + Σi{1→N}( aixi ) = a0 + a1x1 + ... + aNxN

の平均 μy は、

μy = a0 + Σi{1→N}( aiμi ) = a0 + a1μ1 + ... + aNμN

に、各確率変数が互いに独立ならば、分散 σy2 は、

σy2 = Σi{1→N}( ai2σi2 ) = a12σ12 + ... + aN2σN2

になるのでした(「(5) 正規分布」の「4) 標本平均と標本分散」参照)。よって、第 i グループから標本 ( x1(i), x2(i), ... xriN(i) ) を抽出した時、

y=Σi{1→2}( Σk{1→riN}( xk(i) ) ) / N
=Σi{1→2}( ri・Σk{1→riN}( xk(i) / riN ) )

の平均 μy は、各確率変数の平均が μi を平均とする確率分布に従うことから

μy=Σi{1→2}( ri・Σk{1→riN}( μi / riN ) )
=r1μ1 + r2μ2

になります(標本を抽出した時の標本平均 y に対する平均が μy となることに注意して下さい)。また、分散 σy2 は、分散が σ2 である確率分布に対して N 個の標本を抽出した時の標本平均の分散が σ2 / N になることから

σy2=Σi{1→2}( ri2・Σk{1→riN}( σi2 / riN ) )
=( r1σ12 + r2σ22 ) / N

と求めることができます。

「単純無作為抽出法」と「層化抽出法」それぞれの平均と分散を比較すると、平均は等しく、分散は「単純無作為抽出法」の方が r1r2( μ1 - μ2 )2 / N だけ大きくなります。グループが n 個になると、「単純無作為抽出法」の場合、

m = Σi{1→n}( riμi )

s2 = Σi{1→n}( riσi2 ) / N + Σi{1→n}( ri( μi - μ )2 ) / N

となり、「層化抽出法」の場合、

m = Σi{1→n}( riμi )

s2 = Σi{1→n}( riσi2 ) / N

となるので、やはり「単純無作為抽出法」の方が Σi{1→n}( ri( μi - μ )2 ) / N = Σi{1→n}( ri{ μi - Σk{1→n}( rkμk ) }2 ) / N だけ分散が大きくなります。和の中の部分を Δi としてこれを変形すると、

Δi=ri{ μi - Σk{1→n}( rkμk ) }2
=ri{ μi - ( r1μ1 + ... + riμi + ... + rnμn ) }2
=ri{ -r1μ1 - ... + ( 1 - rii - ... - rnμn }2

ここで、1 - ri = r1 + ... + ri-1 + ri+1 + ... + rn なので、

Δi=ri{ r1( μi - μ1 ) + ... + ri( μi - μi ) + ... + rn( μi - μn ) }2
=ri Σk{1→n}( { rk( μi - μk ) }2 )

になります。このことから、各グループに対する分布の平均が離れるほど、「単純無作為抽出法」での抽出では分散が大きくなっていきます。直感的には完全にランダムな抽出をした方が公正で信頼性の高い結果が得られるように感じるのですが、各グループの分布が離れている場合を想像してみると、「単純無作為抽出法」では精度が悪くなっていくというのは納得できるのではないかと思います。


2) ヒストグラム(Histogram)

母集団から抽出した標本の分布を調べることは非常に重要な作業になります。以下に示したデータは、統計局ホームページから入手した、平成 21 年 10 月 1 日での年齢別人口を表にしたものです。

日本の年齢別人口(平成 21 年 10 月 時) (単位 千人)
年齢人口年齢人口年齢人口年齢人口年齢人口年齢人口
01,078201,302401,835602,26680978100以上48
11,092211,347411,800612,24781914
21,084221,388421,793622,13182847
31,072231,414431,407631,33583794
41,050241,463441,746641,43484703
51,088251,490451,636651,74785598
61,111261,494461,593661,69586522
71,145271,478471,541671,73587457
81,160281,490481,522681,68288393
91,180291,551491,534691,52689362
101,179301,589501,564701,32090256
111,193311,653511,521711,40291228
121,188321,698521,481721,42992193
131,183331,783531,560731,41793163
141,206341,869541,613741,34294129
151,208351,966551,614751,25195106
161,190362,002561,717761,2339680
171,212371,964571,812771,1819759
181,216381,918581,922781,1159840
191,253391,864592,068791,0329929

このデータをよく見ると、36 歳と 60 歳の二箇所でピークらしいものがあるのが分かりますが、全体の分布がどのような状態であるかまでを読み取ることは困難です。このようなデータを視覚的に分かりやすくするためによく利用されるのが「ヒストグラム(Histogram)」または「度数分布図」と呼ばれるグラフです。

ヒストグラム

ヒストグラムにすることで、分布の様子がはっきりと分かるようになります。70 歳くらいから、年齢が上がるごとに人口が減少していくのは当然としても、35 歳をピークに年齢が下がるにつれても人口が減少しています。これは、いわゆる「少子化傾向」を示しています。
ピークの他に、極端に人口が少ない年齢層があります。43 歳の人口が極端に少ないのは「丙午(ひのえうま)」の迷信から出産を控えた結果であり、63 〜 64 歳での減少は終戦前後であったことが原因です。二つのピークのうち、戦後のピークは第一次ベビーブームと呼ばれ、この世代は「団塊の世代」と呼ばれます。団塊の世代による出産から第二次ベビーブームが起こり、これが二つめのピークとなります。しかし、90 年代から 2000 年代で発生するはずの第三次ベビーブームは日本では発生しませんでした。

このように、分布を見ることでデータの様々な傾向を読み取ることができます。その傾向にはたいていの場合、理由があるので、その傾向の原因を類推し、分析することが統計学の主な目的であるともいえます。

ヒストグラムの横軸は「階級(Class)」と呼ばれ、一つの階級が、ある範囲のデータや一つの項目などを表します。縦軸は「度数(Frequency)」で、各階級に対するデータの数を示します。階級の決め方によってヒストグラムの形状は変わり、あまり細かく分割すれば各階級の度数が小さくなりすぎて、押しつぶして広げたような形状になってしまうし、逆に分割数が小さすぎると大雑把な分布しか見ることができなくなります。階級の数の目安としてはいろいろな考え方があるようで、例えば次のような方法があります。

スタージェスの公式(Sturges' Formula)
K = 1 + log2N ( K:階級の数 ; N:データ数 )
スコットの選択法(Scott's Choice)
h = 3.5s / N1/3 ( h:階級の幅 ; s:標本標準偏差 )
フリードマン-ディアコニスの方法(Freedman-Diaconis Rule)
h = 2( Q3 - Q1 ) / N1/3 ( Q1:第一四分位 ; Q3:第三四分位 )

「四分位(Quartile)」については後ほど説明します。他にも、データ数の平方根を利用するなど、いくつかの手法があるようです。


階級の決め方として、データの範囲で分ける方法(例えば身長別の分布など)や種類別で分ける方法(例えば国別の分布など)があることを先程簡単に書きましたが、収集したデータはその性質によっていくつかのカテゴリに分類され、それによって扱い方も異なってきます。この分類基準のことを「尺度水準(Level of Measurement)」といいます。尺度水準は、「スタンレー・スティーヴンス(Stanley Smith Stevens)」によって 1946 年の論文「測定尺度の理論について(On The Theory of Scales of Measurement)」の中で提唱され、以下の四つの尺度に分類されます。

比例(比率)尺度(Ratio Scale)
データが連続量であり、それぞれのデータの間隔に意味がある場合、データの大小関係によって順序を決めることができ、しかも各データの間隔からそれぞれがどの程度離れているかを知ることができます。さらに、各データの比にも意味があるような場合は「比例尺度」と呼ばれるカテゴリに分類されます。例えば、長さの場合、2m は 1m の 2 倍ということができるので、これは「比例尺度」にあたります。比率が意味を持つことから、データは全て負数とはならず、値がゼロとなる基準を勝手に決めることもできないので(例えば長さが「ゼロ」となる基準を決めることは通常は行いません)そのことを「絶対零点」を持つといいます。
間隔尺度(Interval Scale)
データどうしの比が意味を持たないことを除けば、比例尺度と同じ性質を持った尺度になります。比例尺度と間隔尺度を併せて「量的データ」と呼ばれ、いわゆる通常の測定データは量的データになります。間隔尺度の代表は温度で、ゼロとなる基準は「摂氏」と「華氏」で異なるようにある程度自由に決められ、マイナスの温度も存在することから比を求めることも意味を持ちません。しかし、絶対温度で考えれば比例尺度として扱うことができるようになります。
順序尺度(Ordinal Scale)
値の大小によって順序が決められるデータが「順序尺度」になります。例えば、100m 走のタイムは比例尺度になりますが、タイムによって付けられた順位は「順序尺度」となります。タイムの差に関係なく順位だけの尺度となるので、例えば 1 位と 2 位の差は 0.01 秒で、2 位と 3 位の差が 1 秒だったとしても、それぞれの間隔の情報は失われてしまうことになります。よって、各データは離散値となり、各データの間隔は意味を持たなくなります。
名義尺度(Nominal Scale ; Categorical Scale)
順序さえも意味を持たなくなった場合は「名義尺度」になります。例えば、国別のデータを調べた時に、各国に対して数値を割り振ったとき、その数値には順序関係はないので「名義尺度」になります。実際には、順序を付けることができても二つのみに分類される場合は名義尺度として扱うそうです。

それぞれの尺度に対する有効な統計量も異なるため、それについては後述します。ヒストグラムの場合、各データに対する度数はどの場合に対しても記述できますが、例えば名義尺度の場合は順序を持たないので、度数の大きい順に並べるなど他の尺度とは異なった表示をする場合もあります。

カテゴリ別にヒストグラム出力を行うためのサンプル・プログラムを以下に示します。

/*
  CategoryHistogram : カテゴリ別ヒストグラム
*/
class CategoryHistogram
{
  static const double MAX_CNT; // ヒストグラムの最大桁数

  map<string, unsigned int> _freq; // 階級と度数

  typedef map<string, unsigned int>::const_iterator freq_cit;

  // 度数によって降順ソートするための比較関数オブジェクト
  struct DescSort : public less<unsigned int> {
    bool operator()( const unsigned int& i, const unsigned int& j )
    { return( i > j ); }
  };

public:

  /*
    コンストラクタ

    const vector<string>& data : 分布を求める対象データ
  */
  CategoryHistogram( const vector<string>& data );

  // 度数を返す
  unsigned int operator[]( string category ) const
  {
    freq_cit cit = _freq.find( category );
    return( ( cit != _freq.end() ) ? cit->second : 0 );
  }

  // 度数の数を返す
  unsigned int size() const { return( _freq.size() ); }

  // ヒストグラムの画面出力
  void print() const;
  void printByFreq() const;
};

/*
  CategoryHistogram でのヒストグラムの最大桁数
*/
const double CategoryHistogram::MAX_CNT = 40;

/*
  CategoryHistogram コンストラクタ

  const vector<string>& data : 分布を求める対象データ
*/
CategoryHistogram::CategoryHistogram( const vector<string>& data )
{
  for ( unsigned int i = 0 ; i < data.size() ; ++i ) {
    if ( _freq.find( data[i] ) == _freq.end() )
      _freq[data[i]] = 1;
    else
      ++( _freq[data[i]] );
  }
}

/*
  CategoryHistogram::print : ヒストグラムの画面出力
*/
void CategoryHistogram::print() const
{
  if ( size() == 0 ) return;

  // 度数の最大値と階級の最大文字数を求める
  unsigned int maxFreq = 0; // 度数の最大値
  unsigned int maxLen = 0;  // 階級の最大文字数
  for ( freq_cit cit = _freq.begin() ; cit != _freq.end() ; ++cit ) {
    if ( maxFreq < cit->second )
      maxFreq = cit->second;
    if ( maxLen < ( cit->first ).length() )
      maxLen = ( cit->first ).length();
  }

  for ( freq_cit cit = _freq.begin() ; cit != _freq.end() ; ++cit ) {
    // 度数の表示桁数
    unsigned int num = MAX_CNT * (double)( cit->second ) / (double)maxFreq;

    // 階級の幅の表示
    printf( "%-*s |", maxLen - ( cit->first ).length() + 1, ( cit->first ).c_str() );

    // ヒストグラムの表示
    for ( unsigned int j = 0 ; j < num ; ++j )
      printf( "%c", '*' );
    printf( "\n" );
  }
}

/*
  CategoryHistogram::printByFreq : ヒストグラムを頻度順に画面出力(パレート図)
*/
void CategoryHistogram::printByFreq() const
{
  // 頻度順に並べた結果を保持する multimap
  multimap<unsigned int,string,DescSort> sortedFreq;
  typedef multimap<unsigned int,string,DescSort>::const_iterator mm_cit;

  if ( size() == 0 ) return;

  // 頻度順にソートしながらコピー
  unsigned int maxLen = 0;  // 階級の最大文字数
  for ( freq_cit cit = _freq.begin() ; cit != _freq.end() ; ++cit ) {
    sortedFreq.insert( pair<unsigned int,string>( cit->second, cit->first ) );
    if ( maxLen < ( cit->first ).length() )
      maxLen = ( cit->first ).length();
  }

  mm_cit cit = sortedFreq.begin();
  unsigned int maxFreq = cit->first; // 最大度数は先頭要素
  for ( ; cit != sortedFreq.end() ; ++cit ) {
    // 度数の表示桁数
    unsigned int num = MAX_CNT * (double)( cit->first ) / (double)maxFreq;

    // 階級の幅の表示
    printf( "%-*s |", maxLen - ( cit->second ).length() + 1, ( cit->second ).c_str() );

    // ヒストグラムの表示
    for ( unsigned int j = 0 ; j < num ; ++j )
      printf( "%c", '*' );
    printf( "\n" );
  }
}

クラス CategoryHistogram は、名義尺度や順序尺度に対して利用することを想定しています。階級は文字列(string)で表され、度数を返すメンバ関数 operator() も階級を表す文字列で指定することになります。
ヒストグラムは、階級の並び順で表示されます。名義尺度の場合、順序に意味を持ちませんが、順序尺度の場合は値の大小によって順序が決まるので、例えば階級を整数値で表すなどして意図した順番に並ぶようにする工夫が必要です。また、順序が意味を持たない場合は度数の大きい順に並べて表示したい場合も想定されるので、そのためのメンバ関数 printByFreq が用意されています。

サンプル・プログラムの中では、STL(Standard Template Library)で提供されているコンテナ・クラスの mapmultimap を利用しています。mapmultimap はともに「連想配列(Associative Array)」と呼ばれ、整数値以外の値を添字として利用できる配列です。添字はキーと呼ばれ、そのキーを指定することで値を抽出します。map はキーが一意になるのに対し、multimap はキーの重複を許した連想配列です。各階級の度数を保持するため、map<string,unsigned int> 型のコンテナ・クラス _freq がメンバ変数として定義されています。テンプレート引数にある string がキーの型、unsigned int が値の型を表し、それぞれ階級とその度数の型になります。従って、階級は文字列で表現し、その度数が unsigned int 型で定義されていることになります。コンストラクタの中で、map を利用した処理を行っていますが、これが典型的な利用方法になります。添字演算 (operator[]) によって、すでにキーが存在していればその値を返し、なければ新たなキーとして要素が追加されます。このとき、キーに対する値としてはその型が持つデフォルト値が代入されるので、値として利用する型は必ず初期値(デフォルト・コンストラクタ)を持たなければなりません。メンバ関数の find は、キーに対する反復子(Iterator)を返します。キーが存在しない場合は _freq.end() (末尾の要素の次を表す反復子) を返すので、これによってキーが存在するかチェックすることができます。キーが存在しなければ

_freq[キー] = 1;

として新たなキーを追加すると同時に 1 を代入し、存在すれば

++( _freq[キー] );

として値に 1 を加えます。このように、mapmultimap の添字演算は要素の追加を伴うので、CategoryHistogram のメンバ関数 operator[] のように、const 句を伴う定数メンバ関数では添字演算を利用することができなくなります。何故なら(たとえ find などによってキーの存在をチェックしていたとしても)、メンバ変数 _freq に対する変更を伴うことになるからです。このような場合は、サンプル・プログラムの中で行っているように、find を使って定数反復子(const_iterator)を取得した上で値を抽出するのが一般的な方法です。

メンバ関数 printByFreq の中では、multimap<unsigned int,string> 型のインスタンス sortedFreq_freq の要素をコピーしています。この時、_freq のキーを値に、また値をキーに入れ替えながらデータを移しているので、度数が重複する場合があることから multimap を利用しているわけです。これによって、度数による並べ替えを実現することができます。しかし、度数の大きい順に並べたいので、要素は降順にしておく必要があります。mapmultimap では、キーによって昇順にソートされた状態で要素が格納されているので、先頭の要素から順に読み込むと自動的にキーが並べられた状態になります。メンバ関数 print ではその性質を利用して、何もしなくてもキーによって並べ替えられた状態で出力されるのですが、printByFreq の場合は少し工夫が必要です。

mapmultimap を構築するとき、各キーの比較方法を別のものに定義することができます。サンプル・プログラムでは、関数オブジェクト less の派生クラスとして DescSort を用意して、逆順にソートされるようにしています。この方法は非常に強力で、他にも様々な利用方法があります。主成分分析のサンプル・プログラムでも利用しており、その中でも若干説明をしてありますので、併せてご覧ください。

度数の大きい順に並べたヒストグラムは、累積度数の比率による折れ線グラフと組み合わせることで「パレート図(Pareto Chart)」と呼ばれる品質管理(QC)用のツールとしてよく利用されます。


連続値に対して範囲別にヒストグラム出力を行うためのサンプル・プログラムを以下に示します。

/*
  extremeValues : 最大値と最小値を求める

  vector<T> x : データ
  T &min, &max : 求めた最大・最小値

  戻り値 : True ... 取得成功 ; False ... x の要素数がゼロ
*/
template<class T> bool extremeValues( vector<T> x, T& min, T& max )
{
  if ( x.size() == 0 ) return( false );

  min = max = x[0];
  for ( unsigned int i = 1 ; i < x.size() ; ++i ) {
    if ( min > x[i] ) min = x[i];
    if ( max < x[i] ) max = x[i];
  }

  return( true );
}

/*
  RangeHistogram : 範囲別ヒストグラム
*/
class RangeHistogram
{
  vector<unsigned int> _freq; // 度数
  double _min;    // 階級の最小値
  double _range;  // 階級一つあたりの幅
  int _precision; // データ精度

  static const string DEFAULT_FORMAT;

  // データを丸める
  double round( double data ) const;

  // データの登録
  void addData( const vector<double>& data, unsigned int freqCnt );

public:

  /*
    コンストラクタ

    const vector<double>& data : 分布を求める対象データ
    unsigned int freqCnt : 階級の数
    double range : 階級一つあたりの幅
    int precision : データ精度
  */
  explicit RangeHistogram( const vector<double>& data, unsigned int freqCnt, int precision = 3 );
  explicit RangeHistogram( const vector<double>& data, double range, int precision = 3 );

  // 度数を返す
  unsigned int operator[]( unsigned int i ) const
  { return( ( i < size() ) ? _freq[i] : 0 ); }

  // 階級の代表値(中央値)を返す
  double repVal( unsigned int i ) const
  { return( ( i < size() ) ? _min + _range * ( (double)i + 0.5 ) : NAN ); }

  // 階級の値(範囲)を返す
  bool range( unsigned int i, double& min, double& max ) const;

  // 階級一つあたりの幅を返す
  double range() const { return( _range ); }

  // 度数の数を返す
  unsigned int size() const { return( _freq.size() ); }

  // ヒストグラムの画面出力
  void print() const { print( DEFAULT_FORMAT ); }
  void print( const string& numFmt ) const;
  void print( const char* numFmt ) const;
};

/*
  RangeHistogram での数値表現デフォルト
*/
const string RangeHistogram::DEFAULT_FORMAT = "%.*e";

/*
  RangeHistogram コンストラクタ

  const vector<double>& data : 分布を求める対象データ
  unsigned int freqCnt : 階級の数
  int precision : データ精度
*/
RangeHistogram::RangeHistogram( const vector<double>& data, unsigned int freqCnt, int precision )
  : _min( 0 ), _range( 0 ), _precision( precision )
{
  double _max; // データの最大値

  if ( _precision < 0 ) _precision = 0;

  if ( freqCnt == 0 ) return;
  if ( ! extremeValues( data, _min, _max ) ) return;

  _range = ( _max - _min ) / (double)freqCnt;

  if ( freqCnt == 1 ) {
    _freq.assign( freqCnt, data.size() );
    return;
  }

  _range = round( _range );
  _min = round( _min );

  // 最大値が境界内に入るように _range を微調整
  double d = _range * pow( 10, -_precision );
  while ( _max > _min + _range * (double)freqCnt )
    _range += d;

  addData( data, freqCnt );
}

/*
  RangeHistogram コンストラクタ

  const vector<double>& data : 分布を求める対象データ
  double range : 階級一つあたりの幅
  int precision : データ精度
*/
RangeHistogram::RangeHistogram( const vector<double>& data, double range, int precision )
  : _min( 0 ), _range( 0 ), _precision( precision )
{
  double _max; // データの最大値

  if ( _precision < 0 ) _precision = 0;

  if ( range <= 0 ) return;
  _range = range;

  if ( ! extremeValues( data, _min, _max ) ) return;

  _min = round( _min );
  _max = round( _max ) + _range * pow( 10, -_precision );

  unsigned int freqCnt = ( _max - _min ) / _range + 1; // 階級の数

  addData( data, freqCnt );
}

/*
  RangeHistogram::round : データの丸め

  _precision に合わせてデータを丸める

  double data : 対象のデータ
*/
double RangeHistogram::round( double data ) const
{
  char s[_precision + 9]; // "-0.xE-000" + \0 : x は _precision の長さ

  snprintf( s, sizeof( s ), "%.*e", _precision, data );

  return( atof( s ) );
}

/*
  RangeHistogram::addData : データの登録

  const vector<double>& data : 対象のデータ
  unsigned int freqCnt : 階級の数
*/
void RangeHistogram::addData( const vector<double>& data, unsigned int freqCnt )
{
  _freq.assign( freqCnt, 0 );

  for ( unsigned int i = 0 ; i < data.size() ; ++i ) {
    unsigned int j = 0;
    for ( double d = _min ; d <= data[i] ; d += _range )
      ++j;
    if ( j > 0 && j <= freqCnt )
      ++( _freq[j - 1] );
  }
}

/*
  RangeHistogram::range : 階級の値(範囲)を返す

  unsigned int i : 階級
  double& min, max : 最小・最大値を得るための変数
*/
bool RangeHistogram::range( unsigned int i, double& min, double& max ) const
{
  if ( i >= _freq.size() ) {
    min = max = NAN;
    return( false );
  }

  min = _min + _range * (double)i;
  max = min + _range;

  return( true );
}

/*
  RangeHistogram::print : ヒストグラムの画面出力

  const char* numFmt : 数値用フォーマット
  string numFmt : 数値用フォーマット
*/
void RangeHistogram::print( const char* numFmt ) const
{
  if ( numFmt == 0 ) {
    print();
  } else {
    string s = numFmt;
    print( s );
  }
}
void RangeHistogram::print( const string& numFmt ) const
{
  const double MAX_CNT = 40; // ヒストグラムの最大桁数

  if ( size() == 0 ) return;

  if ( &numFmt == 0 ) {
    print();
    return;
  }

  // 度数の最大値を求める
  unsigned int maxFreq = (*this)[0];
  for ( unsigned int i = 1 ; i < size() ; ++i )
    if ( maxFreq < (*this)[i] )
      maxFreq = (*this)[i];

  // 数値出力用フォーマットの設定
  string fmt = numFmt + "-" + numFmt + " |";

  // 最初の階級の幅を取得する
  double min, max;
  range( 0, min, max );

  for ( unsigned int i = 0 ; i < size() ; ++i ) {
    // 度数の表示桁数
    unsigned int num = MAX_CNT * (double)( (*this)[i] ) / (double)maxFreq;

    // 階級の幅の表示
    if ( numFmt == DEFAULT_FORMAT )
      printf( fmt.c_str(), _precision, min, _precision, max );
    else
      printf( fmt.c_str(), min, max );

    // ヒストグラムの表示
    for ( unsigned int j = 0 ; j < num ; ++j )
      printf( "%c", '*' );
    printf( "\n" );
    min += range();
    max += range();
  }
}

連続値の場合は、各階級に対して範囲を決めて度数を求める必要があります。コンストラクタは、階級の数と、各階級の範囲のどちらかを指定することができるようにしてあります。対象のデータから最小値と最大値を求め、階級の数が引数として渡されている場合は各階級の範囲を、逆に各階級の範囲が引数として渡されている場合は階級の数を、それぞれ決定します。あとは、各データに対してどの階級に含まれるのかを求めながら度数に加算していく処理を繰り返していきます。
しかし、単純に上記処理を行うだけでは、各階級の範囲はキリのいい数値にならないので、下位を切り捨てる処理を行っています。その精度は precision を使って指定することができます。

下図は、CategoryHistogram を使って適当なサンプル・データをヒストグラムで出力した結果です。

A |****************************
B |****
C |********************
D |************************************
E |****************
F |********************************
G |********
H |****************************************

これを度数の大きい順で描画すると、下図のようになります。

H |****************************************
D |************************************
F |********************************
A |****************************
C |********************
E |****************
G |********
B |****

RangeHistogram を使った出力サンプルは下図のようになります。

38.0-40.7 |********
40.7-43.5 |**********
43.5-46.2 |************************
46.2-48.9 |*******************************
48.9-51.6 |**************************
51.6-54.4 |****************************************
54.4-57.1 |********************
57.1-59.8 |******
59.8-62.6 |*
62.6-65.3 |***

ヒストグラムの別種として「累積度数分布」というものがあります。これは、度数を累積したものを棒グラフで表示したもので、正規分布のように、中央の度数が大きく、そこから離れるに従って度数が小さくなっていく傾向がある場合、そのグラフはちょうど S 字によく似た分布を示すことから「シグモイド曲線 (Sigmoid Curve)」と呼ばれる形になります。


3) 要約統計量(Descriptive Statistics)

ヒストグラムは、分布の状態を視覚的に見ることができる重要な情報源となります。通常は、ヒストグラムを使って分布の特徴を調べ、それを数値で示すために最適な要約統計量を検討することが次の作業になります。要約統計量は大きく「分布の中心的位置を表すもの」と「分布のバラツキ度を表すもの」の二種類があります。「分布の中心的位置を表すもの」としては次のようなものがあります。

○ 平均値

平均値は、今まで何度も登場している代表的な要約値です。全データの和をデータの個数で割ったものは「標本平均(Sample Mean)」と呼ばれ、「中心極限定理」によって、任意の分布に対して標本平均は正規分布に従い、「大数の法則」によって、データの個数が大きくなるに従い標本平均の分散はゼロに近づくことは、「(5) 正規分布」の「3) 大数の法則(Law of Large Numbers)と中心極限定理(Central Limit Theorem)」で紹介しました。

和を個数で割った平均値は、「相加平均」や「算術平均(Arithmetic Mean)」の名で呼ばれることもあります。その他に、「相乗平均」または「幾何平均(Geometric Mean)」と呼ばれる次のような平均値が利用される場合もあります。

G = { Πi{1→N}( xi ) }1/N

データ数を N とした時、各データの積の N 乗根が相乗平均になります。この計算方法から分かるように、データは必ず正値である必要があります。

上式の右辺に対して対数をとると、

log( { Πi{1→N}( xi ) }1/N )=log( Πi{1→N}( xi ) ) / N
=Σi{1→N}( log( xi ) ) / N

と変形することができるので、相乗平均の対数は、データの対数に対する相加平均と等しくなります。相加平均と相乗平均の違いについて簡単な例を挙げると、2 と 8 の相加平均は ( 2 + 8 ) / 2 = 5 になるのに対し、相乗平均は ( 2・8 )1/2 = 4 になります。5 は、2 及び 8 に対する差が等しい位置にあるのに対し、4 は 2 及び 8 に対する比が等しい位置にあります。ある会社の採用者数が、去年は一昨年の 2 倍、さらに今年は去年の 8 倍 だったとき、採用者数は平均して何倍の増加であったかを考えると、例えば一昨年の採用者数が 1 人であれば、去年は 2 人、今年は 16 人となるわけですが、相加平均を使って毎年平均 5 倍ずつと考えると、今年の採用者数は 25 人となって実際より多くなってしまいます。相乗平均で考えれば今年の採用者数は 16 倍の 16 人となって正しい結果が得られます。このように、人口増加率のようなデータなどには相乗平均が用いられます。また、相乗平均の対数はデータの対数に対する相加平均と等しいので、指数関数的に増加・減少するようなデータに対して相乗平均が利用される場合があります。

その他の特殊な平均としては「調和平均(Harmonic Mean)」というものがあり、次のような式で表されます。

H = N / Σi{1→N}( 1 / xi )

つまり、データの逆数に対して相加平均を計算して、その逆数をとったものが調和平均になります。

調和平均は、相加平均では正しい平均値が得られないような場合において用いられることがあります。例えば、トラックで荷物を運搬するとき、行きは荷物のせいで 50 km/s の速度で、帰りは荷物がないので 70 km/s の速度であったとすると、速度の相加平均は ( 50 + 70 ) / 2 = 60 km/s になります。しかし、例えば片道の走行距離が 350km あったとすれば、行きは 7 時間、帰りは 5 時間掛かることになるので、合計 12 時間必要だったことになり、往復の走行距離は 700km なので、平均速度は 700 / 12 ≅ 58.3 km/s で、相加平均とはずれています。一方、調和平均を求めると、2 / ( 1 / 50 + 1 / 70 ) ≅ 58.3 km/s になって、正しい結果が得られます。速さは 距離 / 時間 で求めることができるので、時間が相加平均となるような早さを求めるには調和平均を使う必要があるわけです。実際、先程の例では走行時間の平均が 6 時間であり、350km の距離に 6 時間かかるならばその速さは 58.3 km/s になります。この計算式を変形していくと

350 / 6 = 350 / { ( 7 + 5 ) / 2 } = 2 / { ( 7 + 5 ) / 350 } = 2 / ( 1 / 50 + 1 / 70 )

となって、調和平均の求め方と等しくなります。

「相加平均」「相乗平均」「調和平均」は「ピタゴラス平均(Pythagorean Means)」と呼ばれ、それぞれの頭文字を取って A, G, H と略されるそうです。このとき、データが全て正値ならば

A ≥ G ≥ H

となり、等号はデータが全て等しい場合に成り立ちます(補足 1)。また、「一般化平均(Generalized Mean)」として

Mp = { Σi{1→N}( xip ) / N }1/p

があり、p = 1 ならば「相加平均」、p = -1 ならば「調和平均」、また、p → 0 の極限は「相乗平均」になります。さらに、p → ∞ での極限はデータの最大値、p → -∞ での極限は最小値になります(補足 2)。

○ 中央値

データをその大きさの順に並べた時、ちょうど中央に位置するデータのことを「中央値(Median)」といいます。例えば、

-3, 5, -1, 12, 4

という 5 つのデータがあった時、これを昇順で並べると

-3, -1, 4, 5, 12

となるので、中央値は 4 になります。データの数が偶数になると、ちょうど中央に位置するデータは存在しないので、中央に近い二つの値の相加平均を中央値にします。連続分布の場合、平均値は確率変数 x の期待値 E[x] で表されていましたが、中央値の場合、確率密度関数 p(x) に対して

∫{-∞→m} p(x) dx = 1 / 2 かつ ∫{m→∞} p(x) dx = 1 / 2

を満たすような m が中央値となります。つまり、中央値を境界として、それ以下になる確率とそれ以上になる確率がどちらも等しく 1 / 2 となることを意味します。

中央値は、外れ値と呼ばれる、データ全体が集まっている位置から極端に離れた値がある場合、代表値として平均値よりも適した値になります。例えば、

-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 100000

の平均値は 9090.5 で、全体の代表値としては不適切です。これは、100000 という値が他と比べて極端に大きいことが原因であり、通常は外れ値として除外する必要があります。しかし、中央値は 0 になり、これは代表値として適切な値になります。よって、測定時のミスなどが原因で外れ値が生じた場合、中央値であればそれらを除外する操作をしなくてもすむという利点があります。
また、年収を元に生活水準を調査するとき、年収が非常に高い層によって平均年収が引き上げられ、実態に合わない結果になるような場合もあります。以下の表は、平成 20 年の「民間給与実態統計調査結果」で、年収の分布を表したものです(統計元:国税庁)

年収人数(単位:万人)
男性女性合計
100万円以下82301.1383.1
100万円台196.2488.1684.3
200万円台341.5410.5752
300万円台500.6276.5777.1
400万円台477.1152.9630
500万円台355.379.4434.7
600万円台24437.1281.1
700万円台178.121.1199.2
800万円台121.713.1134.8
900万円台79.77.887.5
1,000越〜1,500万円152.613165.6
1,500超〜2,000万円32.4335.4
2,000超〜2,500万円9.50.810.3
2,500万円越11.1112.1
合計2781.81805.44587.2
年収の分布

度数分布表から平均値を求める場合、各階級の代表値とその度数について積を求め、その和を利用するのが一般的です。それを計算した結果は下表のようになります。

階級代表値代表値 x 度数
(単位:10万円)
100万円以下51915.5
100万円台1510264.5
200万円台2518800
300万円台3527198.5
400万円台4528350
500万円台5523908.5
600万円台6518271.5
700万円台7514940
800万円台8511458
900万円台958312.5
1,000越〜1,500万円12520700
1,500超〜2,000万円1756195
2,000超〜2,500万円2252317.5
2,500万円越2753327.5
合計195959
平均値42.72

一方、分布の累積比率を調べると、

階級人数累積人数累積比率
100万円以下383.1383.10.0835
100万円台684.31067.40.2327
200万円台7521819.40.3966
300万円台777.12596.50.5660
400万円台6303226.50.7034
500万円台434.73661.20.7981
600万円台281.13942.30.8594
700万円台199.24141.50.9028
800万円台134.84276.30.9322
900万円台87.54363.80.9513
1,000越〜1,500万円165.64529.40.9874
1,500超〜2,000万円35.44564.80.9951
2,000超〜2,500万円10.34575.10.9974
2,500万円越12.14587.21.0000

となるので、中央値は 200 万円 〜 300 万円の間にあることが分かります。度数分布表から中央値を求める場合、中央値が存在する階級が am-1 〜 am ならば、各階級の度数を fi ( i = 1, 2, ... )、度数の合計を N としたとき、中央値 M は

M = am-1 + ( am - am-1 ) { N / 2 - Σi{1→m-1}( fi ) } / fm

で求められます。線形補間法と同じような方法で、m - 1 までの累積度数と m までの累積度数との比率から、N / 2 にあたる値を求めて中央値としていることになります。

中央値での線形補間法

よって、年収の中央値は

20 + ( 30 - 20 )( 4587.2 / 2 - 1819.4 ) / 777.1 = 26.10 [ 単位 : 10万円 ]

になります。平均値と中央値の間で大きな差が生じたのは、年収の大きな階級に分布の裾が広がっているためです。その度数がわずかでも、度数の大きな階級の年収との差が大きければ、それだけ平均値は大きくなる傾向にあります。分布を見るかぎり、代表値として適切なのは中央値の方であると考えることができます。同様な現象は、例えば製品の歩留まり(Yield)、すなわち投入量に対して実際に製品として得られた数量の比率にもいえます。何らかの要因で歩留まりの低いロット(Lot ; 製品を管理する単位)が少量ある場合、それに引きずられて平均歩留まりは低くなる傾向にあります。

中央値は「分位数(Quantile)」の中の特別な値です。「q-分位数(q-quantile)」は、昇順に並んだ N 個のデータ列に対して、( N / q ) x i 番目 ( i = 1, 2, ..., q - 1 ) にある値になります。但し、二つのデータの間に位置する場合は、先ほどと同様に、二つのデータを使って線形補間によって適切な値にします。2-分位数の場合は中央値を意味し、4-分位数は「四分位(Quartile)」、100-分位数は「百分位(Percentile)」といって、よく利用される分位数です。四分位は、小さいものから順に「第一四分位(First Quartile)」「第二四分位(Second Quartile)」「第三四分位(Third Quartile)」と呼ばれ、第二四分位は中央値を意味します。これらをグラフ化したものは「箱ひげ図(Box Plot)」と呼ばれ、分布の大まかな形を知るのに役に立ちます。

箱ひげ図

上の図は、箱ひげ図の例を表しています。線(ひげ)で表されている部分は両端が最大・最小値であることを示し、箱の部分の両端が第一・第三四分位になります。また、箱の中の線は中央値を表します。その他に、平均値を併せて示すものや、外れ値を点で表しているもの(その場合、ひげの両端は外れ値を除外した時の最大・最小値を表します)などもあります。

分布が左右対称ならば、平均値と中央値は等しくなります。正規分布はそのような分布の代表です。

○ その他

その他の要約統計量として、「最頻値(Mode)」というものがあり、度数分布表の中で最も度数の高い階級の代表値を最頻値とするのが一般的です。いわゆる「多数決」によって決めた代表値と考えることができます。精度を上げるために、度数の最も高い階級に対して、その前後にある階級の度数を利用して、以下の式のような計算を行う場合もあります。

Mode = am + ( am+1 - am )・{ fm+1 / ( fm-1 + fm+1 ) }

m は度数の最も高い階級の番号を表し、am はその区間の最小値になります。また、fi は i 番目の階級の度数を表します。

「ミッドレンジ(Mid-range)」は最大値と最小値の相加平均によって得られる要約値です。ミッドレンジは、データ内に外れ値や異常値があると意味のない値になってしまうため、最大値と最小値を利用する代わりに第一四分位と第三四分位を利用する場合もあります。この時の値は「ミッドヒンジ(Mid-hinge)」と呼ばれます。


「分布のバラツキ度を表すもの」としては次のようなものがあります。

○ 分散と標準偏差

分散(Variance)も、平均値同様に今まで何回か紹介した重要な要約値です。まず、N 個の標本 xi ( i = 1, 2, ... N ) に対する標本分散(Sample Variance) s2 は、xi の標本平均を m としたとき、

s2 = Σi{1→N}( ( xi - m )2 ) / N

と定義されるのでした。s2 の期待値 E[s2] は

E[s2] = { ( N - 1 ) / N }σ2

となるので、s2 は不偏推定量ではありません。そこで、

u2 = { N / ( N - 1 ) }s2

として、u2 を不偏分散(Unbiased Variance)というのでした。詳しい内容については、「(5) 正規分布」の「4) 標本平均と標本分散」でも説明しています。

母集団の分散については、

σ2 = E[x2] - μ2

が成り立つことを、「(2) 確率空間」の「4) 期待値」で証明しています。つまり、(二乗の平均) と (平均の二乗) の差が分散と等しくなるわけです。標本分散の場合、

s2=Σi{1→N}( ( xi - m )2 ) / N
=Σi{1→N}( xi2 - 2mxi + m2 ) / N
=Σi{1→N}( xi2 ) / N - 2mΣi{1→N}( xi ) / N + m2Σi{1→N}( 1 / N )
=Σi{1→N}( xi2 ) / N - 2m2 + m2
=Σi{1→N}( xi2 ) / N - m2

となるので、ここでも (二乗の平均) と (平均の二乗) の差が分散と等しくなります。但し、ここでの平均は標本平均を指すことに注意して下さい。

分散の平方根は標準偏差(Standard Deviation)となるのでした。分散が二乗の形 σ2, s2, u2 の形で表されるのはそのためで、σ が母集団の標準偏差、s や u が標本標準偏差となります。s と u のどちらが標本標準偏差と呼ばれるかについては統一がされていないようですが、通常は u の方を指すようです。よって、データから標準偏差を求める場合、標本平均の差との二乗和を N ではなく N - 1 で割ったものを使うのが一般的です。u を不偏標準偏差と呼ぶことができればいいのですが、u は不偏推定量とはならないためそうは呼ぶことができません。このあたりについては別の章で詳しく紹介したいと思います。

○ その他

ミッドレンジと近い考え方で「範囲(Range)」というものがあり、最大値と最小値の差から得られます。また、ミッドヒンジと同様の考え方によって「四分偏差(Quantile Deviation)」というものがあり、これは第三四分位と第一四分位の差の半分から求められます。

範囲 R = MAX - MIN

四部偏差 QD = ( Q3 - Q1 ) / 2

但し、MAX : 最大値, MIN : 最小値, Qi : 第 i 四分位

その他に、「平均偏差(Mean Deviation)」というものがあります。これは、平均との差に対して二乗を計算する代わりに差の絶対値を計算して平均を求めた値です。

平均偏差 MD = Σi{1→N}( | xi - m | ) / N

この値は絶対値を含むため数学的な取り扱いが難しく、あまり利用されることはないようです。

また、バラツキの度合いを比較するような場合、両者の平均値や単位などが異なればそのまま比較することはできないので、「変動係数(Coefficient of Variation)」を利用する場合があります。これは、標準偏差を平均値で割った値です。

変動係数 CV = s / m

変動係数は通常、100 倍して % で表されます。


各尺度水準に対して利用できる要約値の一覧を以下に示します。○は利用可能なもの、△は条件付きで利用できるもの、- は利用できないものをそれぞれ表します。

尺度水準中心的位置を表すものバラツキ度を表すもの
相加平均相乗平均調和平均分位数最頻値ミッドレンジ分散/標準偏差平均偏差変動係数範囲
比例尺度
間隔尺度---
順序尺度---
名義尺度---------

負数が存在すれば相乗平均や調和平均は意味のある値にはならないので、これらが利用できる尺度は比例尺度のみとなります。また、比例尺度であってもデータにゼロがあれば、相乗平均は常にゼロであり、調和平均は計算することができないので、データが全て正数であるという条件が必要です。

変動係数も比例尺度のみで扱うことができます。例えば、摂氏 -10 度と摂氏 20 度の相加平均は摂氏 5 度で、標準偏差(ここでは不偏分散ではなく標本分散の平方根とします)は摂氏 15 度なので、変動係数は 300 % となるのに対し、絶対温度で(摂氏 0 度 = 273K として)考えると、相加平均は 278K、標準偏差は 15K で変動係数は 5.40 % になります。摂氏で考えた場合は間隔尺度になりますが、このときゼロ点は任意に決めることができるので、ゼロの位置によって平均値は変化してしまいます。よって、絶対零点を持つ比例尺度でなければ正しい変動係数を得ることはできないことがこの例から分かります。
変動係数は、複数の群に対してバラツキの度合いを比較することが目的なので、変動係数が有効であるかを見分けるには、各群の平均値と標準偏差を計算して、平均値に対して標準偏差が比例して大きくなるかを調べれば簡単にわかります。間隔尺度の場合、標準偏差は平均値に対して一定の値になるので、平均値が大きくなるほど変動係数は小さくなってしまいます。また、通常は比例尺度で扱えるものでも、場合によっては変動係数が意味をなさない場合があります。例えば身長計によって身長を測定した時、身長が異なっていても測定の誤差は変化しないと考えられます。すると、一人に対して複数回身長を測定してその測定誤差を見積もるような場合は、変動係数を利用することはできません。しかし、例えば子供と大人の群に分けてバラツキを評価する場合、身長が比例尺度であると考えれば変動係数を利用する方が正しい結果が得られることになります。

順序尺度や名義尺度では、平均値は意味を持ちません。また、平均がないことから分散や標準偏差なども利用できません。順序尺度の場合、平均値の代わりに中央値などの分位数や最頻値を利用することができます、名義尺度では順位も持たないために中央値も利用できず、最頻値のみが利用可能です。また、名義尺度では最大・最小値の概念もないので、ミッドレンジや範囲も利用できません。しかし、順序尺度の差に「意味がある」と見なして平均値を求めることはよく行われるようです。順位を付けるために使った実際の量の差がほぼ等しい(例えば 100m 走のタイム差はほとんど同程度)と考えれば、通常の測定値のように扱うことはできないことはないので、実際のデータを見て(また、そのデータがどのようなものかも判断材料として)臨機応変に対処するのが通例のようです。

尺度水準は、比例尺度・間隔尺度 > 順序尺度 > 名義尺度の順で情報量が少なくなっていきます。比例尺度の「身長」は順序尺度の「背の順」に置き換えることが可能で、さらにそれを二つのグループ「背の高いグループ」と「背の低いグループ」にすれば、名義尺度として扱われます。このように、情報量の多いものから少ないものへ変換することはできますが、その逆は不可能です。


今までに紹介した要約統計量を計算するためのサンプル・プログラムを以下に示します。

// 一変数の関数オブジェクト
template<class T> class UnaryFunc {
public:
  virtual T operator()( const T& x ) = 0;
};

// 底を e とする対数を求める関数
class Log : public UnaryFunc<double>
{
public:
  double operator()( const double& x )
  { return( log( x ) ); }
};
class LogL : public UnaryFunc<long double>
{
public:
  long double operator()( const long double& x )
  { return( logl( x ) ); }
};

// 逆数を求める関数
template<class T> class Reciprocal : public UnaryFunc<T>
{
public:
  T operator()( const T& x )
  { return( 1 / x ); }
};

// a との差の二乗和を求める関数
template<class T> class Deviation : public UnaryFunc<T>
{
  T _a;

public:

  Deviation( T a )
    : _a( a ) {}

  T operator()( const T& x )
  { return( ( x - _a ) * ( x - _a ) ); }
};

// a との差の絶対値を求める関数
template<class T> class AbsDeviation : public UnaryFunc<T>
{
  T _a;

public:

  AbsDeviation( T a )
    : _a( a ) {}

  T operator()( const T& x )
  { return( ( x > _a ) ? ( x - _a ) : ( _a - x ) ); }
};

/*
  sum : データの総和を求める

  vector<T>& x : データ

  誤差を軽減する方法は Kahanの加算アルゴリズムを利用

  戻り値 : 求めた総和(データ数がゼロならゼロを返す)
*/
template<class T> T sum( const vector<T>& x )
{
  T ans = 0; // 求める総和
  T r = 0;   // 積み残し
  for ( unsigned int i = 0 ; i < x.size() ; ++i ) {
    T y = x[i] - r;         // 計算値から積み残しを引いた値
    T buff = ans + y;       // 今までの合計値に y を加える
    r = ( buff - ans ) - y; // 積み残しを計算
    ans = buff;
  }

  return( ans );
}
template double sum( const vector<double>& x );
template long double sum( const vector<long double>& x );

/*
  sum : データを変換した上での総和を求める

  vector<T>& x : データ
  UnaryFunc<T>& func : データ変換用関数オブジェクト

  誤差を軽減する方法は Kahanの加算アルゴリズムを利用

  戻り値 : 求めた総和(データ数がゼロならゼロを返す)
*/
template<class T> T sum( const vector<T>& x, UnaryFunc<T>& func )
{
  T ans = 0; // 求める総和
  T r = 0;   // 積み残し
  for ( unsigned int i = 0 ; i < x.size() ; ++i ) {
    T y = func( x[i] ) - r; // 計算値から積み残しを引いた値
    T buff = ans + y;       // 今までの合計値に yを加える
    r = ( buff - ans ) - y; // 積み残しを計算
    ans = buff;
  }

  return( ans );
}
template double sum( const vector<double>& x, UnaryFunc<double>& func );
template long double sum( const vector<long double>& x, UnaryFunc<long double>& func );

/*
  sampleAverage : 標本平均を求める

  vector<T>& x : データ

  戻り値 : 求めた平均値(データ数がゼロならゼロを返す)
*/
template<class T> T sampleAverage( const vector<T>& x )
{
  T s = sum( x );
  if ( x.size() > 0 ) s /= (T)x.size();

  return( s );
}
template double sampleAverage( const vector<double>& x );
template long double sampleAverage( const vector<long double>& x );

/*
  sampleAverage : データを変換した上での標本平均を求める

  vector<T>& x : データ
  UnaryFunc<T>& func : データ変換用関数オブジェクト

  戻り値 : 求めた平均値(データ数がゼロならゼロを返す)
*/
template<class T> T sampleAverage( const vector<T>& x, UnaryFunc<T>& func )
{
  T s = sum( x, func );
  if ( x.size() > 0 ) s /= (T)x.size();

  return( s );
}
template double sampleAverage( const vector<double>& x, UnaryFunc<double>& func );
template long double sampleAverage( const vector<long double>& x, UnaryFunc<long double>& func );

/*
  harmonicAverage : 調和平均を求める

  vector<T>& x : データ

  戻り値 : 求めた調和平均(総和がゼロならゼロを返す)
*/
template<class T> T harmonicAverage( const vector<T>& x )
{
  Reciprocal<T> reciprocal; // 逆数を求める関数オブジェクト

  T s = sum( x, reciprocal );
  if ( s != 0 ) s = (T)x.size() / s;

  return( s );
}
template double harmonicAverage( const vector<double>& x );
template long double harmonicAverage( const vector<long double>& x );

/*
  median : 中央値(メディアン)を求める

  vector<T> x : データ

  戻り値 : 求めた中央値(データ数がゼロならゼロを返す)
*/
template<class T> T median( vector<T> x )
{
  std::sort( x.begin(), x.end() );
  unsigned int i = x.size();

  if ( i % 2 != 0 ) {
    return( x[i/2] );
  } else if ( i == 0 ) {
    return( 0 );
  } else if ( i == 1 ) {
    return( x[0] );
  } else {
    return( ( x[i/2-1] + x[i/2] ) / (T)2 );
  }
}
template double median( vector<double> x );
template long double median( vector<long double> x );

/*
  quantile : m-分位数を求める

  vector<T> x : データ
  vector<T>& q : m-分位数を格納する変数(サイズが m - 1 を表す)

  戻り値が False の場合はデータは変更されない

  戻り値 : True ... 取得成功 ; False ... x の要素数がゼロ
*/
template<class T> bool quantile( vector<T> x, vector<T>& q )
{
  // 等分する数を求める
  unsigned int qSize = q.size();
  if ( qSize++ == 0 ) return( true );

  // 階級の数(データの数 - 1) を求める
  unsigned int xSize = x.size();
  if ( xSize-- == 0 ) return( false );

  std::sort( x.begin(), x.end() );

  for ( unsigned int i = 1 ; i < qSize ; ++i ) {
    T t = (T)( xSize * i ) / (T)qSize;
    unsigned int j = t;
    t -= (T)j;
    q[i - 1] = x[j];
    if ( j < xSize ) q[i - 1] += ( x[j + 1] - x[j] ) * t;
  }

  return( true );
}
template bool quantile( vector<double> x, vector<double>& q );
template bool quantile( vector<long double> x, vector<long double>& q );

/*
  extremeValues : 最大値と最小値を求める

  vector<T> x : データ
  T &min, &max : 求めた最大・最小値

  戻り値 : True ... 取得成功 ; False ... x の要素数がゼロ
*/
template<class T> bool extremeValues( vector<T> x, T& min, T& max )
{
  if ( x.size() == 0 ) return( false );

  min = max = x[0];
  for ( unsigned int i = 1 ; i < x.size() ; ++i ) {
    if ( min > x[i] ) min = x[i];
    if ( max < x[i] ) max = x[i];
  }

  return( true );
}
template bool extremeValues( vector<double> x, double& min, double& max );
template bool extremeValues( vector<long double> x, long double& min, long double& max );

/*
  midRange : ミッドレンジを求める

  const vector<T>& x : データ

  戻り値 : ミッドレンジ( x の要素数がゼロの場合はゼロを返す)
*/
template<class T> T midRange( const vector<T>& x )
{
  T min, max;

  if ( ! extremeValues( x, min, max ) )
    return( 0 );

  return( ( min + max ) / (T)2 );
}
template double midRange( const vector<double>& x );
template long double midRange( const vector<long double>& x );

/*
  midHinge : ミッドヒンジを求める

  const vector<T>& x : データ

  戻り値 : ミッドヒンジ( x の要素数がゼロの場合はゼロを返す)
*/
template<class T> T midHinge( const vector<T>& x )
{
  vector<T> q( 3 );

  if ( ! quantile( x, q ) ) return( 0 );

  return( ( q[0] + q[2] ) / (T)2 );
}
template double midHinge( const vector<double>& x );
template long double midHinge( const vector<long double>& x );

/*
  sampleVariance : 標本分散を求める

  const vector<T>& x : データ
  bool isUnbiased : 不偏分散として計算するか

  戻り値 : 標本分散または不偏分散( x の要素数が足りない場合はゼロを返す)
*/
template<class T> T sampleVariance( const vector<T>& x, bool isUnbiased )
{
  if ( x.size() == 0 ) return( 0 );
  if ( isUnbiased && x.size() < 2 ) return( 0 );
  T m = sampleAverage( x );

  Deviation<T> dev( m );
  T d = sum( x, dev );
  T n = ( isUnbiased ) ? x.size() - 1 : x.size();

  return( d / n );
}
template double sampleVariance( const vector<double>& x, bool isUnbiased );
template long double sampleVariance( const vector<long double>& x, bool isUnbiased );

/*
  range : 範囲を求める

  const vector<T>& x : データ

  戻り値 : 範囲( x の要素数がゼロの場合はゼロを返す)
*/
template<class T> T range( const vector<T>& x )
{
  T min, max;

  if ( ! extremeValues( x, min, max ) )
    return( 0 );

  return( max - min );
}
template double range( const vector<double>& x );
template long double range( const vector<long double>& x );

/*
  quantileDeviation : 四部偏差を求める

  const vector<T>& x : データ

  戻り値 : 四部偏差( x の要素数がゼロの場合はゼロを返す)
*/
template<class T> T quantileDeviation( const vector<T>& x )
{
  vector<T> q( 3 );

  if ( ! quantile( x, q ) ) return( 0 );

  return( ( q[2] - q[0] ) / (T)2 );
}
template double quantileDeviation( const vector<double>& x );
template long double quantileDeviation( const vector<long double>& x );

/*
  meanDeviation : 平均偏差を求める

  const vector<T>& x : データ

  戻り値 : 平均偏差( x の要素数がゼロの場合はゼロを返す)
*/
template<class T> T meanDeviation( const vector<T>& x )
{
  if ( x.size() == 0 ) return( 0 );
  T m = sampleAverage( x );

  AbsDeviation<T> absDev( m );
  T d = sampleAverage( x, absDev );

  return( d );
}
template double meanDeviation( const vector<double>& x );
template long double meanDeviation( const vector<long double>& x );

/*
  CategoryHistogram::mode : 最頻値(Mode)を求める

  戻り値 : 最頻値
*/
string CategoryHistogram::mode() const
{
  string candidateMode;     // 最頻値の候補
  unsigned int maxFreq = 0; // 度数の最大値

  for ( freq_cit cit = _freq.begin() ; cit != _freq.end() ; ++cit ) {
    if ( maxFreq < cit->second ) {
      maxFreq = cit->second;
      candidateMode = cit->first;
    }
  }

  return( candidateMode );
}

/*
  RangeHistogram::mode : 最頻値(Mode)を求める

  戻り値 : 最頻値(階級の代表値(中央値))
*/
double RangeHistogram::mode() const
{
  double candidateMode = NAN; // 最頻値の候補
  unsigned int maxFreq = 0;   // 度数の最大値

  for ( unsigned int i = 0 ; i < size() ; ++i ) {
    if ( maxFreq < (*this)[i] ) {
      maxFreq = (*this)[i];
      candidateMode = repVal( i );
    }
  }

  return( candidateMode );
}

平均を求める場合にはデータの総和が必要になるため、総和を求める関数 sum が用意されています。sum には一変数関数オブジェクト UnaryFunc 型のインスタンスを渡し、データに対して任意の変形を行った上で総和を求めることができるバージョンがあり、相乗平均や調和平均を求める時に利用されています。
相乗平均は、各データの総積から求める必要があるので、非常に大きな値になる可能性があります。よって、各データの対数に対する相加平均を求めることができるようにして、相乗平均そのものを計算するための関数は用意していません。そのため、対数を計算しながら総和が求められるように、Log, LogL といった専用の関数オブジェクトが用意されています。相乗平均の形にしたい場合は、求めた値を指数として e のべき乗を求めれば得ることができます。調和平均は、各データの逆数に対して相加平均を求めることによって計算するため、逆数を計算するための関数オブジェクト Reciprocal が用意されています。

分位数は、中央値を求めるための median と任意の分位数を求めるための quantile が用意されています。quantile では、等分する数を、分位数を得るための vector である q のサイズから得ています。等分数は、( q のサイズ + 1 )になることに注意して下さい。例えばサイズが 3 の場合、等分は 4 となります。また、各データの間の数は、データ数から 1 を引いたものになります。100 個のデータに対しては、データ間の数は 99 です。

いくつかの要約統計量を求めるためには最大・最小値が必要です。そのために extremeValues という関数が用意されており、ミッドレンジと範囲はこの関数が利用されています(同じ関数が範囲別ヒストグラム出力用サンプル・プログラムでも使用されています)。また、ミッドヒンジと四部偏差では quantile を使って四分位を求めています。

分散は sampleVariance で計算することができます。第二引数の isUnbiasedtrue ならば、不偏分散を求めるために偏差の二乗和を ( データ数 - 1 ) で割り、false なら標本分散としてデータ数で割っています。

最後に、頻度はヒストグラムから得ることになるため、前に紹介した CategoryHistogramRangeHistogram のメンバ関数として用意しています。


抽出したデータを元にその分布形状をヒストグラムで調べたり、要約統計量を求めてその状態を確認することによって、データが持つ性質を導き出すことができるようになります。今までに例として紹介したような、年齢別の人口や年収の分布などは、その国や都市などが抱えている問題(子供の減少や貧困層の増大など)を浮き彫りにします。「統計学」を意味する "Statistics" の語源はラテン語で「国家」を意味する "Statisticum" で、「ゴットフリート・アッヘンワル(Gottfried Achenwall」が、国の状態を調べるための学問である「国家の科学(Science of State)」として紹介したのがはじまりのようです。抽出したデータを単に眺めているだけでは何も分からないので、視覚的に分かりやすくしたり、代表となる値を求めるなどの工夫を行って、少ない情報でできるだけ多くの状態を知ることができるようにしたのが初期の統計学で、これを「記述統計学(Descriptive Statistics)」と呼びます。やがて確率論が誕生すると、それが統計学に応用され、「カール・ピアソン(Karl Pearson)」によって記述統計学は大成されました。


4) 積率(Moment)

確率分布の平均や分散を計算するときに便利な関数として「積率母関数(Moment Generating Function)」を以前紹介しました。

g(θ) = E[eθx]

積率母関数の r 階導関数 g(r)(θ) は、θ = 0 のとき E[xr] に等しくなります。

g(r)(0) = E[xr]

この性質を利用すると、平均値や分散を簡単に計算することができる場合があります。「積率母関数」の名の由来は、この関数が係数を E[xr] とする多項式で表されるため E[xr] の母関数となることによります。つまり、

g(θ) = E[1] + E[x]・x + E[x2]・x2 / 2! + ... + E[xr]・xr / r! + ...

とすれば g(r)(0) = E[xr] となり、この形は数列 ar = E[xr] を係数とする「指数型母関数」に一致します。この係数 E[xr] のことを r 次の「積率(モーメント ; Moment)」といいます。

「a のまわりの分散」は E[( x - a )2] で、特に a = μ = E[x] の場合を単に「分散」というのでした。これと同様に、a のまわりの r 次積率を

μr(a) = E[( x - a )r]

といい、特に a = μ1(0) = E[x] = μ の場合は「中心積率(Central Moment)」といいます。文献によっては中心積率のことを単に積率と表す場合もあるようです。

以上のことから、下記内容が導かれます。

  1. μ0(0) = 1
  2. μ1(0) = μ (平均)
  3. μ2(μ) = σ2 (分散)

標本から求めた積率は「標本積率(Sample Moment)」となります。これを mr とすると、

mr = Σi{1→N}( xir ) / N

になり、中心積率の場合は

mr(m) = Σi{1→N}( ( xi - m )r ) / N

と求めることができます。

平均や分散は、積率の中の特別なものと考えることができます。よって、積率や中心積率自体が要約統計量の一種と見なすこともできます。積率から得られる要約値として、他にも「歪度(Skewness)」と「尖度(Kurtosis)」があります。

確率変数 x を、次の式で変換する操作を「標準化(Normalization)」といい、変換した結果 z は「標準得点(Standard Score)」と呼ばれます。

z = ( x - μ ) / σ

これはちょうど、平均が 0、分散が 1 になるように確率変数を変換する操作と同じことになります。実際、確率密度関数 p(x) に対して E[x] = μ、E[(x - μ)2] = σ2 とすれば、

μz = E[z]=E[( x - μ ) / σ]
=E[x] / σ - E[1]・μ / σ
=μ / σ - μ / σ = 0
σz2 = E[( z - μz )2]=E[( x - μ )2 / σ2]
=E[( x - μ )2] / σ2
=σ2 / σ2 = 1

となります。μz と σz2 はそれぞれ、z に対する一次及び二次の中心積率を表しています。このとき、三次の中心積率 E[z3] と四次の中心積率 E[z4] をそれぞれ「歪度(Skewness)」、「尖度(Kurtosis)」といいます。また一般に、z に対する r 次の中心積率のことを「標準化積率(Standardized Moment)」といいます。

標準正規分布 N( 0, 1 ) = ( 1 / ( 2π )1/2 ) exp( -x2 / 2 ) の積率母関数 g(θ) は

g(θ) = E[eθx]=( 1 / ( 2π )1/2 ) ∫{-∞→∞} eθx exp( -x2 / 2 ) dx
=( 1 / ( 2π )1/2 ) ∫{-∞→∞} exp( -( x2 - 2θx ) / 2 ) dx
=( 1 / ( 2π )1/2 ) ∫{-∞→∞} exp( -{ ( x - θ )2 - θ2 } / 2 ) dx
=( 1 / ( 2π )1/2 ) exp( θ2 / 2 ) ∫{-∞→∞} exp( -( x - θ )2 / 2 ) dx

積分値において t = ( x - θ ) / √2 とすると、積分範囲は ( -∞, ∞ ) のままで dt = dx / √2 なので、

∫{-∞→∞} exp( -t2 ) √2・dt = ( 2π )1/2

となって(ガウス積分の計算方法参照)、

g(θ) = exp( θ2 / 2 )

となることから、

g(1)(θ)=θ・exp( θ2 / 2 )
g(2)(θ)=exp( θ2 / 2 ) + θ2・exp( θ2 / 2 )
g(3)(θ)=θ・exp( θ2 / 2 ) + 2θ・exp( θ2 / 2 ) + θ3・exp( θ2 / 2 )
=3θ・exp( θ2 / 2 ) + θ3・exp( θ2 / 2 )
g(4)(θ)=3exp( θ2 / 2 ) + 3θ2・exp( θ2 / 2 ) + 3θ2・exp( θ2 / 2 ) + θ4・exp( θ2 / 2 )
=3exp( θ2 / 2 ) + 6θ2・exp( θ2 / 2 ) + θ4・exp( θ2 / 2 )

よって、g(1)(0) = 0、g(2)(0) = 1、g(3)(0) = 0、g(4)(0) = 3 になります。g(3)(0) = E[x3] よりこれは歪度、g(4)(0) = E[x4] よりこれは尖度を表すので、標準正規分布の歪度は 0、尖度は 3 になります。

平均 μ に対して、各確率変数との差 ( x - μ ) を求めてその総和を計算すると、結果は必ずゼロになります。確率分布が直線上に並んでいたとして、平均に位置するところに支点をおけば、ちょうど釣り合うようなイメージを持つとその意味が理解できると思います。平均に対してより小さなデータとより大きなデータがあって、それぞれが打ち消しあうことでバランスが保たれているわけです。しかし、分散は差の二乗の和を使っているため、全ての値が正値になります。平均から離れたデータがあるほど差の二乗は大きくなるため、分布がどの程度のバラツキを持っているかを知る指標として使われるわけです。
歪度は差の三乗との和から得られます。従って、各要素は正と負の値が混在しますが、差の総和とは異なり、差が大きいものほど三乗した値はより大きくなるため総和はゼロにはなりません。例えば、平均値より大きい側に分布が広がっていれば、総和は正の値をとるようになります。このことから、歪度は分布の対称性を見る場合に利用されます。正規分布の場合、平均を中心に左右対称となっていることから、歪度はゼロになるわけです。
尖度が高い分布は、通常は平均付近にあるピークが鋭く、周辺にある裾の部分は長く厚めになるのに対し、尖度の低い分布は逆に、ピークが緩やかで裾の部分が短く薄めになります。ピークが鋭く裾が厚ければ、差の小さい要素が少なく差の大きい要素は多い傾向となって、四乗和の総和は大きくなるので、逆にピークが緩くて裾が薄い分布と比べると、どちらも分散を 1 に正規化しているのなら尖度はより大きくなります。標準正規分布の尖度は 3 であり、正規分布の標準化が標準正規分布なので、正規分布の尖度は 3 です。正規分布の尖度を基準として、尖度から 3 を引いたものを「過剰尖度(Excess Kurtosis)」といい、過剰尖度が 0 のものを "mesokurtic"、正値のものを "leptokurtic"、負値のものを "platykurtic" と呼ぶそうです。

正規分布を標準としていることからも分かるように、歪度と尖度は正規分布として近似できるかどうかを判定するためによく利用されているようです。

標準化積率を求めるためのサンプル・プログラムを以下に示します。

// 標準得点のべき乗を求める関数
template<class T> class Moment : public UnaryFunc<T>
{
  T _m; // 平均
  T _s; // 標準偏差
  unsigned int _r; // 次数

public:

  Moment( T m, T s, unsigned int r )
    : _m( m ), _s( s ), _r( r ) {}

  T operator()( const T& x )
  {
    if ( _r == 0 ) return( 1 );
    T t = ( x - _m ) / _s;
    T ans = t;
    for ( unsigned int i = 1 ; i < _r ; ++i )
      ans *= t;

    return( ans );
  }
};

/*
  stdMoment : 標準化積率を求める

  const vector<T>& x : データ
  unsigned int r : 次数

  戻り値 : 標準化積率( x の要素数または標準偏差がゼロの場合はゼロを返す)
*/
template<class T> T stdMoment( const vector<T>& x, unsigned int r )
{
  if ( x.size() == 0 ) return( 0 );
  T m = sampleAverage( x );
  T s = sqrt( sampleVariance( x, false ) );
  if ( s == 0 ) return( 0 );

  Moment<T> moment( m, s, r );
  T ans = sampleAverage( x, moment );

  return( ans );
}
template double stdMoment( const vector<double>& x, unsigned int r );
template long double stdMoment( const vector<long double>& x, unsigned int r );

積率は標準得点のべき乗の平均値で得られるので、要約統計量のサンプル・プログラムと同じ方法で計算することができます。


今回は、標本の抽出とその分布の調査方法をメインに紹介しました。これらは「記述統計学」の分野であり今では古典の部類に入ります。次回から「推測統計学」の分野の内容を中心に紹介する予定です。


補足1) ピタゴラス平均

正値のみからなる任意のデータに対する「相加平均(A)」「相乗平均(G)」「調和平均(H)」の間には

A ≥ G ≥ H

の関係式が成り立ちます。データが一つだけの場合は明らかなので、二つの場合を証明してみます。対象のデータを x, y ( x > 0, y > 0 ) としたとき、

A = ( x + y ) / 2

G = ( xy )1/2

H = 2 / ( 1 / x + 1 / y )

になります。この時、A > 0, G > 0, H > 0 が成り立ちます。A2 - G2 を計算すると

A2 - G2=( x + y )2 / 4 - xy
=( x2 + 2xy + y2 ) / 4 - xy
=( x2 - 2xy + y2 ) / 4
=( x - y )2 / 4 ≥ 0

となるので、A ≥ G が成り立ちます(等号は x = y のとき)。また、x' = 1 / x, y' = 1 / y とすれば、1 / H と 1 / G はそれぞれ x' と y' の相加平均、相乗平均を表すので 1 / H ≥ 1 / G が成り立ち、従って、G ≥ H です。以上から、

A ≥ G ≥ H

が成り立ちます。次に、2p 個のデータ x[i] ( i = 1, 2, ... 2p ; xi > 0 ) までは関係式 A ≥ G が成り立つと仮定します。2p+1 個のデータに対する相加平均 A は

A=( x[1] + x[2] + ... + x[2p+1] ) / 2p+1
={ ( x[1] + x[2] + ... + x[2p] ) / 2p + ( x[2p + 1] + x[2p + 2] + ... + x[2p+1] ) / 2p } / 2

と変形することができます。データを二つに分けて、2p 個のデータによる相加平均を作ったわけです。それぞれの相加平均に対して A ≥ G が成り立つと仮定しているので

A{ ( x[1]・x[2]・ ... ・x[2p] )1/2p + ( x[2p + 1]・x[2p + 2]・ ... ・x[2p+1] )1/2p } / 2
{ ( x[1]・x[2]・ ... ・x[2p] )1/2p・( x[2p + 1]・x[2p + 2]・ ... ・x[2p+1] )1/2p }1/2
=( x[1]・x[2]・ ... ・x[2p+1] )1/2p+1 = G

となります。よって、帰納法によって 2p 個のデータに対しては A ≥ G が成り立つことが証明されます。任意の個数のデータについては、2p-1 < N ≤ 2p となる p に対して M = 2p - N ( M + N = 2p ) としたとき、S = Σi{1→N}( x[i] ) / N として

( x[1] + x[2] + ... + x[N] + Σi{1→M}( S ) ) / ( M + N ) ≥ ( x[1]・x[2]・ ... ・x[N]・SM )1/(M+N)

が成り立ちます。S を余分に M 回加えることで、データ数を 2 のべき乗にしているわけです。すると、x[1] + x[2] + ... + x[N] = NS、Σi{1→M}( S ) = MS なので左辺は S と等しくなり、

S ≥ ( x[1]・x[2]・ ... ・x[N]・SM )1/(M+N)

両辺を M + N 乗すると

SM+N ≥ x[1]・x[2]・ ... ・x[N]・SM

SM > 0 なので、SM で割ってから N 乗根を取れば

S ≥ ( x[1]・x[2]・ ... ・x[N] )1/N

S は相加平均の計算式そのものなので、A ≥ G が任意のデータ数に対して成り立つことが証明されました。

データ数が 2 の場合で示したように、調和平均と相乗平均の関係式は相加平均と相乗平均の関係式から証明することができるので、任意のデータ数に対して A ≥ G ≥ H であることが証明できたことになります。


補足2) 一般化平均

一般化平均は次の式で表されます。

Mp = { Σi{1→N}( xip ) / N }1/p

p = 1 のとき、

M1 = Σi{1→N}( xi ) / N

となって相加平均を表し、p = -1 のとき、

M-1={ Σi{1→N}( xi-1 ) / N }-1
=N / Σi{1→N}( 1 / xi )

なので調和平均になります。ここで、一変数の関数に対するマクローリン展開(Maclaurin Series)

f(x) = f(0) + f'(x) x + f''(x) x2 / 2 + ... + f(k)(x) xk / k! + ...

となり、x が 0 に非常に近ければ二次項以降を無視して

f(x) ≅ f(0) + f'(x) x

と近似することができます。f(p) = ln( Σi{1→N}( xip ) / N ) とすると、

f'(p) = ( Σi{1→N}( xip ) / N )-1・Σi{1→N}( xip・ln( xi ) ) / N

より f(p) のマクローリン展開から

f(p)ln( Σi{1→N}( xi0 ) / N ) + { ( Σi{1→N}( xi0 ) / N )-1・Σi{1→N}( xi0・ln( xi ) ) / N }・p
={ Σi{1→N}( ln( xi ) ) / N }・p

よって、p → 0 のとき ln( Σi{1→N}( xip ) / N ) / p → Σi{1→N}( ln( xi ) ) / N になります。これを利用すると、

Mp=exp( ln( Σi{1→N}( xip ) / N ) / p )
exp( Σi{1→N}( ln( xi ) ) / N ) ( p → 0 )
={ Πi{1→N}( xi ) }1/N

となって、p → 0 での極限は相乗平均になることが分かります。

xi の最大値を xMAX、最小値を xMIN とすると、

Mp={ Σi{1→N}( xip ) / N }1/p
{ Σi{1→N}( xMAXp ) / N }1/p = xMAX
{ Σi{1→N}( xMINp ) / N }1/p = xMIN

が成り立ちます。Σi{1→N}( xip ) = xMAXp{ 1 + Σi{1→N;i≠MAX}( ( xi / xMAX )p ) } より

Mp=xMAX{ 1 + Σi{1→N;i≠MAX}( ( xi / xMAX )p ) }1/p / N1/p
xMAX ( p → ∞ )

また、Σi{1→N}( xip ) = xMINp{ 1 + Σi{1→N;i≠MIN}( ( xi / xMIN )p ) } より同様に Mp → xMIN ( p → -∞ ) が成り立ち、一般化平均がデータの最大値と最小値の範囲を超えることはなく、極限において最大値・最小値となることがわかります。


<参考文献>
1. 「確率・統計入門」 小針あき宏著 (岩波書店)
2. 「統計数学入門」 本間鶴千代著 (森北出版)
3. 「プログラミング言語 C++ 第3版」 Bjarne Stroustrup著 (アジソン・ウェスレイ出版)
4. 我楽多頓陳館 --- 統計学入門
統計学について分かりやすく紹介されています。非常に参考になりました。
5. 「微分積分いい気分」 --- 相加平均と相乗平均の関係
相加平均と相乗平均の大小関係の証明はここから拝借しました。他の証明方法も詳しく紹介されています。
6. 「IIJIMASの日記」 --- 「平均」について
一般化平均の p → 0 での極限が「相乗平均」になることの証明が紹介されています。
7. Wikipedia

[Go Back]前に戻る [Back to HOME]タイトルに戻る
inserted by FC2 system