実測データを使い最小二乗法で直線の信頼区間と予測区間を求めよう

# ■■■　最小二乗法　■■■
import numpy as np
import random
import matplotlib.pyplot as plt
from scipy.stats import t

smh = 158.4 # 平均（母平均）:２１才、女子身長
ssh = 4.5 # 標準偏差（母標準偏差）:２１才、女子身長

smw = 52.6 # 平均（母平均）:２１才、女子体重
ssw = 7.5 # 標準偏差（母平均）:２１才、女子体重

bmimw = 21 # 平均（母平均）:２０才台、女子BMI
bmisw = 2.9 # 標準偏差（母平均）:２０才台、女子BMI

sn = 2000 # 母数
en = 20 # 標本数
r = 0.5 #相関係数

# 17	158.4	4.5	16	52.6	7.5

h = np.random.normal(loc=smh, scale=ssh, size=sn)
h2w = [(p*p/10000, random.gauss(bmimw, bmisw) * p*p /10000) for p in h]
#BMI = w/h^2 = 22
#https://www.e-stat.go.jp/dbview?sid=0003224178

fig, ax = plt.subplots() # Figureオブジェクトとそれに属する一つのAxesオブジェクトを同時に作成

u,v = zip(*h2w)
A = np.vstack((np.array(u), np.ones(len(u)))).T
vv = np.array(v).T
uu, residuals, rank, s = np.linalg.lstsq(A, vv, rcond=None)
au = uu[0]
bu = uu[1]

#au = np.cov(u, v)[0][1] / np.var(u)
#bu = np.mean(v) - au * np.mean(u)
u1 = np.array([np.mean(u) - np.std(u)*3 , np.mean(u) + np.std(u)*3])
v1 = bu + au * u1  # 回帰曲線
ax.annotate("$w=$"+ "{:.1f}".format(au) + "$h^2+$" + "{:.1f}".format(bu)\
            ,xy = (np.mean(u) , np.mean(v))\
            ,xytext = (np.mean(u)+ np.std(u), np.mean(v)- 3*np.std(v))\
            ,arrowprops=dict(arrowstyle="->"))

x2y = random.sample(h2w, en) #無作為抽出
x,y = zip(*x2y)

A = np.vstack((np.array(x), np.ones(len(x)))).T
vv = np.array(y).T
uu, residuals, rank, s = np.linalg.lstsq(A, vv, rcond=None)
a = uu[0]
b = uu[1]
#a = np.cov(x, y)[0][1] / np.var(x)
#b = np.mean(y) - a * np.mean(x)
ax.annotate("$w=$"+ "{:.1f}".format(a) + "$h^2+$" + "{:.1f}".format(b)\
            ,xy = (np.mean(x) , np.mean(y))\
            ,xytext = (np.mean(x)+ np.std(x),np.mean(y)+ 3*np.std(y))\
            ,color="m",arrowprops=dict(arrowstyle="->"))

x1 = np.array([np.mean(x) - np.std(x)*3 , np.mean(x) + np.std(x)*3])
y1 = b + a * x1  # 回帰曲線
y_hat = [b + a * p for p in x]

alpha = 0.05 # 有意水準
n = len(x)
tvalue = t.ppf(1 - alpha / 2, n - 2) # 自由度n-2の時のt値(alphaを1/2することに注意)
ve = (np.sum(np.array(y) ** 2) - np.sum(np.array(y_hat) ** 2)) / (n - 2)
x_bar = np.mean(x)

conf = np.arange(np.mean(u) - np.std(u)*3 , np.mean(u) + np.std(u)*3, np.std(u)/16)
conf_hi = [b + a * p + tvalue * np.sqrt(1 / n + (p - x_bar) ** 2 / np.sum((x - x_bar) ** 2)) * ve for p in conf]
conf_lo = [b + a * p - tvalue * np.sqrt(1 / n + (p - x_bar) ** 2 / np.sum((x - x_bar) ** 2)) * ve for p in conf]

pred = np.arange(np.mean(x) - np.std(x)*3 , np.mean(x) + np.std(x)*3, np.std(x)/16)
pred_hi = [b + a * p + tvalue * np.sqrt(1 + 1 / n + (p - x_bar) ** 2 / np.sum((x - x_bar) ** 2)) * ve for p in pred]
pred_lo = [b + a * p - tvalue * np.sqrt(1 + 1 / n + (p - x_bar) ** 2 / np.sum((x - x_bar) ** 2)) * ve for p in pred]

ax.scatter(u,v,alpha=0.2)
ax.plot(u1,v1)
ax.plot(conf, conf_hi, 'k:', label='conf_95%')
ax.plot(conf, conf_lo, 'k:')
ax.fill_between(conf, conf_hi, conf_lo, facecolor='lightcoral', alpha=0.3)

ax.plot(pred, pred_hi, 'k:', label='pred_95%')
ax.plot(pred, pred_lo, 'k:')
ax.fill_between(pred, pred_hi, pred_lo, facecolor='lightgreen', alpha=0.2)

ax.scatter(x, y,color = "m")
ax.plot(x1, y1,color = "m")

ax.set_xlabel("Height$^2$ $h^2$ / m$^2$")
ax.set_ylabel("Weight $w$ / kg")

plt.show()

# ■■■　最小二乗法　■■■