単回帰分析~コンビニエンスストア

概要

Pythonを使った統計計算と図示の練習のため、コンビニエンスストアで単回帰分析をやってみた。

コンビニの店舗数は「商業動態統計年報」の2016年データを使い、説明変数として2015年の国勢調査人口、2018年の国土地理院による都道府県面積、2017年道路統計年報の道路実延長データを使った。

計算コードは以下の通り。

人口との関係

コンビニ店舗数と人口の散布図と回帰式を以下に示す。

相関係数が極めて高いのは当然で、やはりコンビニの出店計画には人口ファクターが強く作用していることがわかる。

定数項bが店舗数のオーダーに比べてほぼゼロというのも興味深い。

係数aの値からは5万人弱で1店舗ということになるが、人口10万くらいの都市で2店舗しかないことになり、ちょっと少ないような気がする。もしかすると、一定規模以上の市町村単位や都市単位くらいで層化して出店計画を立てているのかもしれない。

面積との関係

店舗数と面積の関係

次に店舗数と都道府県面積の関係を見てみた。

この結果はかなりはずれで、データを見ても人口が少なく面積が群を抜いて大きい北海道の影響を大きく受けている。

ここで、面積が極端に大きい北海道(83423.83㎢~1位、2906店舗~5位)と面積が小さいが集積度が極端に高い東京都(2193.96㎢~45位、7003店舗~1位)の2つを除いて計算してみる。

これはさらにおかしな結果で、面積が小さいほど店舗数が多いことになる。

考えてみれば、集客数を期待するなら人口が集積している地域が有利だから、人口密度に比例する可能性を考えた方がいいのかもしれない。もし面積が小さい県の方が集積度が高いと想定すると、面積だけを取り出したときに逆の関係になるとも考えられるが、相関係数や決定係数が小さすぎるので考察は難しい。

人口密度

以下は人口密度との関係。

今度はかなりきれいに相関の高さが出ている。

直接的な計算式に入れているかどうかわからないが、GISなどで出店計画を立てるとしたら、人口密度の高いエリアを選んでいくだろうことが想定される。

ただ、店舗数は人口などといった売り上げに直結するデータから導かれるのが普通で、人口密度が高くても人口が少なければ出店インセンティブにはならない。

人口と人口密度の関係

試しに人口を説明変数、人口密度を被説明変数として両者の関係を見てみると、驚くことに「人口が多いほど人口密度が高くなる(あるいはその逆)」という関係になる。

ここから先は人口論や地域論になりそうなので置いておくが、少なくとも日本においては、「狭いところほど人が集まっている傾向がある」ということになりそうである。

もちろんこれは他の国でも一般に当てはまることかもしれないが、朝のラッシュ時に特定の車両に無理やり乗り込んでいる割に離れた車両がすいているとか、1本電車を遅らせたらガラガラだったとか、そのあたりの行動パターンを見ていると、何となく日本に特有のような気がする。

道路延長との関係

最後に道路延長との関係を見てみる。

東京のように稠密な都市は例外とすると、概ね関係はありそうである。ただし相関係数、決定係数は高くない。

コンビニ店舗が道路の利便性に依っていることは推測できるが、やはり人口という売り上げ直結のデータに比べると関係は弱い。

高速道路の延長についても見てみたが、こちらはほとんど関係は見られなかった。

ただ、高速道路の延伸に伴ってコンビニエンスストアの店舗数が伸びているようであり、マクロな延長というよりも物流上のインパクトが大きいことは予想される。

 

相関係数

相関係数の定義

相関係数は、多数のデータの組がどの程度線形に近い性質を持つかを表す値で、以下で定義される。

(1)    \begin{eqnarray*} r &=& \frac{{\rm Cov}(X, Y)}{\sqrt{V(X) \cdot V(Y)}} \\ &=& \frac{E(X - \overline{X})(Y - \overline{Y})} {\sqrt{ ( E\left[ (X - \overline{X})^2 \right] E\left[ (Y - \overline{Y})^2 \right] }} \\ &=& \frac{\displaystyle \sum_{i=1}^{n}(x_i - \overline{x})(y_i - \overline{y})} {\left[ \displaystyle \sum_{i=1}^{n} (x_i - \overline{x})^2 \displaystyle \sum_{i=1}^{n} (y_i - \overline{y})^2 \right]^{1/2}} \end{eqnarray*}

相関係数の線形変換

変数を線形変換した場合の相関係数は分散・共分散の性質から、以下のようになり、元のままか符号が反転する。

(2)    \begin{eqnarray*} r' &=& \frac{{\rm Cov}(aX+b, cY+d)}{\sqrt{V(aX+b) \cdot V(cY+d)}} \\ &=& \frac{ac{\rm Cov}(X, Y)}{\sqrt{a^2 V(X) \cdot c^2 V(Y)}} \\ &=& \frac{ac}{|ac|} \cdot r \end{eqnarray*}

完全線形関係の場合の相関係数

XとYが完全な線形関係にある場合、相関係数は1または-1になる。このとき、傾きの大きさや、平行移動は相関係数に影響しない。

(3)    \begin{eqnarray*} r &=& \frac{ {\rm Cov}(X, aX+b) } { \left[ V(X) \cdot V(aX+b) \right]^{1/2} } \\ &=& \frac{ a {\rm Cov}(X, X) } { \left[ V(X) \cdot a^2 V(X) \right]^{1/2} } \\ &=& \frac{a}{\sqrt{a^2}} \cdot \frac{ {\rm Cov}(X, X) } { \left[ V(X) \cdot V(X) \right]^{1/2} } \\ &=& \frac{a}{|a|} \cdot \frac{ V(X) } { V(X) } \\ &=& 1 \; {\rm or} \; -1 \end{eqnarray*}

いろいろな分布の相関係数

完全な線形関係

以下のコードで確認。

相関係数は傾きや平行移動に対して影響を受けず、増加関数なら1、減少関数なら-1になることがわかる。

線形性が強い関係

線形関数の値に対して、乱数でばらつきを与えた場合の相関係数の違いを示す。ばらつきが大きい方が相関係数は小さくなる。

負の線形性が強い場合は、相関係数がマイナスになる。

放物線(強い関係があるのに相関が低いケース)

以下のような放物線では、XとYにきちんとした数学的関係があるのに、相関係数がゼロに近くなる。

相関係数はXとYが単調増加/単調減少の度合いが強いほど、また線形関係に近いほど1に近くなるが、それ以外の関係が強い場合にはそれを補足できない場合がある。

反比例(負の線形性に見えてしまう場合)

以下は反比例関数の場合。

関数の形状や範囲に寄るが、この場合は相関係数の絶対値が0.8以上と1に近く、これだけ見ると負の線形性が強そうに見える。

対数関数(正の線形性に見えてしまう場合)

対数関数の場合。この場合は0.9以上とかなり強い線形性を示唆している。

相関係数に関する注意

本来の関係との乖離

先にみたように、線形関係ではないが数学的な関係を持つ場合に、相関係数からは全く関係がない、元の関係とは異なり線形関係を持つ、といった解釈になることがある。

相関係数が高い場合に、線形回帰式などで物事を予測する際には注意が必要。

変数間に解析的な関係が見いだせるならそれを重視すべきであり、よしんばそれがわからないにしても、定義域の範囲で「ある程度は当たる」程度に考えておくべきか。

因果関係

堂々と間違えられるケースが、「科学的な」記事やマスメディアなどでよくみられる。

相関係数は「変数間の単調な増加/減少傾向が強いかどうか」だけを示すもので、必ずしも因果関係を示唆しない。

気温が高いとビールはよく売れるがおでんは売れない。その二つに負の相関があるからといって、「ビール好きはおでんが嫌い」と言えないが、形を変えてこのような解釈がなされる恐れがある。

もともとのメカニズムで因果関係が示唆されていて、その上で相関係数の大きさを論じるなら意義もあるが、その場合でも、事象に対する寄与度などをよく考えておかないと「それだけが原因」と考えるような間違いを犯すことになる。

 

Python3 – 浮動小数点の誤差

概要

Pythonで浮動小数点誤差に行き当たった。

コンソール上にはErrorが表示され計算結果はnanになるが、配列を使っていたので、問題がどこで生じているかわかるまでにちょっと手間取った。

解析上の考え

中心が(cx, cy)、半径がrの円を考え、x座標を与えてそれに対する円上のy座標を計算しようとしていたとき。

このような円の方程式は、陰関数では以下のようになる。

     \begin{equation*} (x - c_x)^2 + (y - c_y)^2 = r^2 \end{equation*}

これをy ≥ 0で解くと、

     \begin{equation*} y = c_y + \sqrt{r^2 - (x - c_x)^2} \end{equation*}

問題の発生

中心の座標が(0.5, 0.5)、半径が0.3となるような円が欲しかったので、xの定義域を0.3~0.8として次のようなコードをPythonで書いた。実際はy < 0の部分も必要で、円周上だけでなく内部の座標をランダムに計算し、それらの相関係数を計算するというものだったが、本質的には変わらない。

その結果、cor=nanとなってしまう。実行時にsqrtの計算で警告が出ている。

そこで配列x、yを表示させてみると、yの最後の要素がnanになっている。

さらに、sqrtの中を表示させてみる。

最後のところで差し引かれる方の値がごくわずかに大きくなっていて、その結果根号の中が負の値となってnanが発生したらしい。

解決策

linspaceの範囲を円周上の上限・下限から少し内側にすることも考えられるが、今回はnumpy.around()を使った。

numpy.around(a, decimals=0)

これは値aを小数点以下decimalsの桁数で丸めてくれるもので、今回は小数点以下10桁目で丸めてみた。

この結果、正常に計算。

注意点

  • 解析上値が等しくなる(差し引きゼロになる)ような場合でも、浮動小数点ではごくわずかな誤差のため、エラーとなることがある
  • numpy.linspace()を使うと上限・下限の間を等分してくれるのでありがたいが、上限・下限自体が他の値と厳密に等しくなければならない場合には、浮動小数点の誤差が生じるため注意が必要

 

Python3 – numpy

全般

numpyのインストール

randomモジュール

配列・ndarray

配列の統計計算(一次元配列二次元配列

 

numpy – 統計量と組み込み関数の注意点

概要

Pythonのnumpyで統計量を扱う場合、分散・共分散で注意を要する点がある。Excelの関数でも同様だが、分散・共分散が標本値からそのまま計算した値か、不偏推定量として計算した値か、ということを意識する必要がある。

ここでは、ndarrayから律義に計算した値と、numpyの組み込み関数から計算した値を比較してみる。

標本数・総和・平均

これはどうしても変わらず、標本として与えた配列の要素数、全要素の和、総和を要素数で除した値。

分散・標準偏差の食い違い

var()、std()は標本の分散・標準偏差

分散・共分散については、提供されたメソッドvar()std()は標本数で除した標本分散・共分散となる。

cov()で与えられるのは不偏推定量

numpy.cov()で二つの一次元配列の分散・共分散行列を得ることができる。分散共分散行列は

     \begin{equation*} \left[ \begin{array}{ll} {\rm Var}(X) & {\rm Cov}(XY) \\ {\rm Cov}(XY) & {\rm Var}(Y) \end{array} \right] \end{equation*}

で定義されるが、以下の計算結果の対角要素は、先の分散の計算結果(2.0)と異なっている。

これらの値は、分散を標本分散ではなく不偏推定量で計算しているためで、試しに先の計算を不偏分散で計算してみると結果は整合する。

ちなみに標本分散は偏差の二乗和を標本数nで割り、不偏分散はn-1で割る。

共分散にも注意

上記では分散について確認したが、共分散についても注意が必要。

共分散は一般に以下で定義される。

     \begin{equation*} {\rm Cov}(X, Y) = E((X - \overline{X}) (Y - \overline{Y})) = \frac{1}{n} \sum_{i=1}^{n} (x_i - \overline{x}) (y_i - \overline{y}) \end{equation*}

これに従って計算した結果。

この結果は、numpy.cov()で計算した非対角要素2.5と食い違う。

ここでn→n-1として計算しなおすと結果は2.5となり同じ値となる。

共分散用は、たとえば回帰分析の残差変動を計算するときに利用するが、その時には標本の共分散を使うので、注意が必要。

numpyで統計計算をする際、特に共分散については、律義に定義式から計算するのがよさそう。

 

numpy – 配列の統計計算(二次元配列)

概要

二次元配列についても、一次元配列と同様な統計関係のメソッド群がある。

二次元についても、同じ機能のメソッドがnumpy、ndarrayのメソッドとして準備されている。

二次元配列の場合は、最小・総和・平均などの計算を全要素/列単位/行単位のいずれで行うかを区別する。具体例として、以下の3×3配列の最小値min()を全要素、列単位、行単位で計算した例を示す。

引数axisを指定しないと全要素が対象となり、axis=0を指定すると列ごとの計算結果を一次元の配列に、axis=1を指定すると行ごとの計算結果を一次元の配列にして返す。

axis=1の場合、計算は行ごとに行われるが、結果は列ベクトルではなく行ベクトルの配列となる。

また、numpyのメソッドではなくndarrayのメソッドでも同じように使える。

最小値・最大値

最小値/最大値を返す。

numpy.min(m, axis=0/1)
numpy.max(m, axis=0/1)
m.min(axis=0/1)
m.max(axis=0/1)

結果は同じなので、numpyのメソッドについて実行例を示す。

総和・総積

一次元配列の全要素の和・積を返す。

numpy.sum(m, axis=0/1)
numpy.prod(m, axis=0/1)
m.sum(axis=0/1)
m.prod(axis=0/1)

平均・分散・標準偏差

一次元配列の要素の平均、分散、標準偏差を返す。分散は標本分散なので、不偏分散が必要な場合はvar()*n/(n-1)とする(ただしn=len(v))。

numpy.mean(m, axis=0/1)
numpy.var(m, axis=0/1)
numpy.std(m, axis=0/1)
m.mean(axis=0/1)
m.var(axis=0/1)
m.std(axis=0/1)

累積和・累積積

一次元配列の要素について先頭から累積して積・和を計算し、それらを要素とする配列を返す。

二次元配列でaxis=0を指定した場合、列ごとに行方向に累積した値が並べられた二次元配列となり、axis=1を指定した場合は行ごとに列方向に累積した値が並べられた二次元配列となる。

numpy.cumsum(m)
numpy.cumprod(m)
m.cumsum()
m.cumprod()

 

numpy – 配列の統計計算(一次元配列)

概要

一次元配列について、様々な統計関係の計算をするメソッド群。

同じ機能のメソッドがnumpy、ndarrayのメソッドとして準備されている。例えば一次元配列vについて最小値を求めるメソッドは、numpy.min(v)v.min()のいずれも同じ結果を返す。

最小値・最大値

一次元配列の要素のうち最小値/最大値を返す。

numpy.min(v)
numpy.max(v)
v.min()
v.max()

結果は同じなので、numpyのメソッドについて実行例を示す。

総和・総積

一次元配列の全要素の和・積を返す。

numpy.sum(v)
numpy.prod(v)
v.sum()
v.prod()

平均・分散・標準偏差

一次元配列の要素の平均、分散、標準偏差を返す。分散は標本分散なので、不偏分散が必要な場合はvar()*n/(n-1)とする(ただしn=len(v))。

numpy.mean(v)
numpy.var(v)
numpy.std(v)
v.mean()
v.var()
v.std()

     \begin{eqnarray*} \overline{V} &=& \frac{1 + 2 + 3 + 4}{4} = 2.5 \\ {\rm Var}(V) &=& \frac{(1-2.5)^2+(2-2.5)^2+(3-2.5)^2+(4-2.5)^2}{4} \\ &=& \frac{2.25+0.25+0.25+2.25}{4} = \frac{5}{4} \\ &=& 1.25 \\ \sigma_V &=& \sqrt{1.25} = 1.1180 \ldots \end{eqnarray*}

累積和・累積積

一次元配列の要素について先頭から累積して積・和を計算し、それらを要素とする配列を返す。

numpy.cumsum(v)
numpy.cumprod(v)
v.cumsum()
v.cumprod()

 

numpy – 配列要素の演算

要素に対するスカラー演算

配列の要素に対するスカラー演算は、それぞれの要素に作用。

二次元配列も同じ。

なお、配列が後ろに来ても要素ごとに演算。

要素ごとの演算

同じ形状の配列同士の演算は、それぞれの要素ごとの演算。

要素ごとの比較

配列同士の比較は、要素ごとの比較結果の配列。

関数演算

numpyの関数の引数が配列の場合、要素ごとの演算結果となる。

 

numpy – ファイルの保存と読み込み

バイナリファイル

バイナリファイルとしての保存にはnp.save()、その読み込みにはnp.load()を使う。基本的なオプションは以下の通り。

numpy.dave("ファイル名", 配列)
配列 = numpy.load("ファイル名")

注意点として、ファイルの拡張子は.npy固定。

  • saveのファイル名に拡張子をつけない場合、自動的に拡張子.npyが付加される
  • saveのファイル名に別の拡張子を書いても、その後ろに.npyが付加される
  • loadのファイル名の拡張子は.npyでなければならない(違う場合はFileNotFoundError)

テキストファイル

概要

テキストファイルとしての保存にはnp.savetxt()、その読み込みにはnp.loadtxt()を使う。扱える配列の次元は1次元か2次元のみ。基本的なオプションは以下の通り。

numpy.savetxt(
    fname,         # ファイル名を任意の拡張子まで指定
    X,             # 書き込む配列
    fmt='%.18e',   # 書き込む書式
    delimiter=' ', # 行内の数値の区切り
    newline='n',   # 改行文字
    header='',     # ヘッダー文字列
    footer='',     # フッター文字列
    comments='# ', # コメントの開始文字
    encoding=None  # エンコーディング
)

numpy.loadtxt(
    fname,                 # ファイル名(拡張子まで)
    dtype=<class 'float'>, # データの型
    comments='#',          # コメント開始文字
    delimiter=None,        # 区切り文字(デフォルトはスペース)
    converters=None,
    skiprows=0,            # ファイル先頭から指定した行数だけ読み飛ばす
    usecols=None,          # タプルで指定した列のみ読み込み
    unpack=False,
    ndmin=0,
    encoding='bytes'max_rows=None)

なお、numpy.ndarrayはすべての要素の型が同じでなければならず、リストやタプルのように要素の型を混在させることはできない

数値データ

数値データのみを扱う場合、整数でも自動的に実数型に変換されて読み込まれる。

実行例

このときのファイルの内容は以下のようになっている。

3次元以上の配列を書き込もうとするとエラー。

文字列

文字列データについては、書き込む場合はfmt="%s"、読み込む場合にはdtype="unicode"を指定。

このとき、text.txtファイルの中は以下のようになっている。

CSVファイルの読み込み

以下のようなCSVファイルを配列に読み込む。

numpy.arrayはすべての要素の型が同じである必要があるため、上記のように文字列と数値が混在したファイルを一つの配列に読み込もうとするとエラーになる。

都道府県名と2次元のデータを別々に読み込むことは可能。

読み込み時にUnicodeエラーが出た場合

ファイル読み込み時に以下のようなエラーが出た場合。

引数にencoding='utf8'を指定して回避。

その他のメソッド

savez()/load()
複数のファイルを非圧縮で扱う。
savez_compressed()/load()
複数のファイルを圧縮して扱う
ndarray.tofile()
配列のストレージへの書き込み。