統計学入門のお供 Julia (1) - 標準得点、偏差値得点 - Qiita
統計学入門のお供 Julia (2) - ピアソンの積率相関係数 - Qiita
統計学入門のお供 Julia (3) - 順位相関係数と自己相関係数 - Qiita
統計学入門のお供 Julia (4) - 最小二乗法と決定係数 - Qiita
統計学入門のお供 Julia (5) - ベイズの定理 - Qiita
統計学入門のお供 Julia (6) - ド・メレの問題 = Oiita
統計学を勉強するため「統計学入門」（東京大学出版会）を読み始めました。練習問題があるのですが、これも勉強中のJuliaを使って計算していきたいと思います。
#１．最小二乗法 least squares method
(x,y)の２次元データの散布図を描き、そのグラフが直線として見れるときに、１次関数の y= bx + a の b, a を求める方法を最小二乗法といいます。
具体的には、以下のLを最小にするa,bを求めていきます。
\begin{align}
&L=\sum_{i=1}^n\{y_i-(bx_i+a)\}^2
\tag{1}\\
\\
&最小値（極小値）では、それぞれの偏微分が0になる。\\
&a,bで偏微分する\\
\\
\\
\\
&na + (\sum x_i)b = \sum y_i \tag{2}\\
&(\sum x_i)a + (\sum x_i^2)b = \sum x_i y_i \tag{3}\\
\\
\\
&これを解くと\\
\\
\\
&b = \frac{\sum x_i y_i - n\bar{x}\bar{y}}{\sum x_i^2 -n\bar{x}^2} \tag{4}\\
&a = \bar{y} - b\bar{x} \tag{5}\\
&(何故ならxとyの平均値より： \quad \bar{x}=\sum x_i/n, \quad \bar{y}=\sum y_i/n \quad)
\\
\\
\\
&\sum(x_i-\bar{x})(y_i-\bar{y})=\sum x_i y_i - n\bar{x}\bar{y} \tag{6}\\
&\sum(x_i - \bar{x})^2 = \sum x_i^2 - n\bar{x}^2 \tag{7}\\
\\
\\
\\
&b = \frac{ \sum(x_i - \bar{x})(y_i - \bar{y})}{\sum(x_i - \bar{x})^2} \tag{8}\\
\\
&つまり
\\
\\
&b = \frac{ ( \sum(x_i - \bar{x})(y_i - \bar{y}) )/n}{ ( \sum(x_i - \bar{x})^2 ) /n}
= \frac{C_{xy}}{(S_x)^2}
= \frac{共分散}{xの標準偏差^2} \tag{9}
\\
\\
\\
&またピアソンの積率相関係数の定義は以下の通りだから\\
&r = \frac {\sum (x_i - \bar{x}) (y_i - \bar{y})} {\sqrt {\sum(x_i-\bar{x})^2} \sqrt {\sum(y_i-\bar{y})^2}} \tag{10}
\\
\\
&式(8)と式(10)を比べて、bとrの関係は以下の通りになる。
\\
\\
&b = r \frac{S_y}{S_x} \tag{11}
\\
\\
&\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \\
\end{align}
式(9)によれば、bはjuliaの以下の式で簡単に計算できます。
b = cov(x,y) / std(x)^2
統計学入門のお供 Julia (2) - ピアソンの積率相関係数 - Qiitaで計算した「都道府県別の自民党支持率と持ち家比率」の例をここでも考えてみます。最小二乗法で１次関数 y= bx + a の b と aを求めます。
using Statistics
using Plots
x=[41.4,76.3,59.2,51.8,52.5,53.2,62.4,55,57.7,63.2,37.5,48.5,32.4,20.5,
47.9,68.9,68.5,52.5,63.3,58.8,59.7,48.4,40.7,51,50.9,34.3,25.8,32.1,
34.4,55.1,60.3,57,45.6,54.2,55.1,55.7,70.3,61.8,47.6,42.5,71.3,55.2,
65.2,42.9,54.7,62,48.2]
y=[52.8,71.2,72.6,63.7,81.3,81.8,70.9,74,73.2,72.9,66.7,65.7,43.7,55.5,
79.6,85.7,75.3,80.5,73,77,77.5,69.2,60,78.2,79.5,61.8,49.6,59.6,72.1,
71,76.3,72.8,71.8,60.7,67,71.8,71.2,68.3,68.5,54.8,76,65.8,69.4,66.9,
69.7,71.2,59]
b=cov(x,y) / std(x)^2
a=mean(y)-b*mean(x)
println(b)
println(a)
x1=0
y1=a
push!(x,x1)
push!(y,y1)
x2=80
y2=b*x2+a
push!(x,x2)
push!(y,y2)
println(x)
println(y)
plot(x,y,
seriestype=:scatter,
title="Jimin-House(1983)",
lab="House",
xlabel = "a (Jimin)",
ylabel = "b (House)",
legend=:topleft
)
最小二乗法で散布図に最も適した直線を求めました。つまり１次関数の傾きbと切片を求めたわけです。その上で、この一次関数の直線が以下の２点を通ること計算させて、散布図に一緒にプロットさせました。赤の直線はこの２点を通るように手動で描いたものです。散布図と直線を一緒のグラフに描く方法がわかりませんでしたので。
(0, 45.0703), (80, 82.1926)
一応printlnの出力は以下の通りです。b と aの値が出力されています。
0.4640285720454987 # bの値
45.07033826566744 # aの値
[41.4, 76.3, 59.2, ..., 62.0, 48.2, 0.0, 80.0]
[52.8, 71.2, 72.6, ..., 71.2, 59.0, 45.0703, 82.1926]
最小二乗法については以下のサイトもわかりやすいです。
最小二乗法（直線）の簡単な説明
https://mathtrain.jp/leastsquares
二変数の二次関数
https://mathtrain.jp/quadratic
方針：最小→極小→微分が0という数2の考え方を使います
#２．決定係数
\begin{align}
&\hat{y_i} = bx_i + a \\
&最小二乗法で求めた直線上の点 \hat{y_i} と実データ y_i とがどの程度違うのか、\\
&その距離を計算することにする。\\
\\
\\
&\sum d_i^2 \equiv \sum(y_i - \hat{y_i})^2 = (1-r^2) \sum (y_i-\bar{y})^2 \tag{12}\\
\\
&rは式(10)のもの。\\
\\
\\
&r^2を決定係数という。\\
&r^2が１に近づくほど、 y_i は \hat{y_i} に近づく。\\
&つまりr^2が１に近づくほど、従属変数yは独立変数xによって決定されている。\\
&図で言えば、r^2が１に近づくほど、散布図は直線に集まっている。\\
\\
&\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad\\
\end{align}
#３．ブートストラップ
P65の問3.4を解きます。
この問題をここに掲載するのが良いかわからないけど成り行きでここに掲載します。
ブートストラップとは、得られたサンプルから何回もリサンプリングを行い、それらのリサンプルから統計量を求めることで、推測を行うような方法のことである。
現段階ではブートストラップの意味を理解していないけど、問題にある手順をプログラム化し、ヒストグラムを描いてみました。
x=[71, 68, 66, 67, 70, 71, 70, 73, 72, 65, 66]
y=[69, 64, 65, 63, 65, 62, 65, 64, 66, 59, 62]
plot(x,y,seriestype=:scatter)
function myrandom11(x,y)
n = 11
idx=rand(1:n, 11)
xx = zeros(n)
yy = zeros(n)
for i in 1:n
j=idx[i]
xx[i] = x[j]
yy[i] = y[j]
end
cor(xx,yy)
end
function do200(x,y)
n = 200
rr = zeros(n)
for i in 1:n
rr[i] = myrandom11(x,y)
end
rr
end
rr = do200(x,y)
histogram(rr , bins=10)
histogramはPlotsの関数です。度数とか自動的に計算してくれるようで優れモノですね。ちょっと感動しました。
今回は以上です。