「入門 信頼性工学」第５章
福井泰好「入門 信頼性工学（第２版）」森北出版 の読後メモ。
設例について、Python ［Google Colaboratory］による演算処理例を示す。
5章 信頼性データの統計的解析
5.0 概要
p75 観測データの態様
p75 分布の推定方法
- 線形回帰分析
- 最尤法
p75 分布の種類
\[ \left\{ \begin{array}{l} 正規分布\left\{ \begin{array}{l} 正規分布 \\ 対数正規分布 \end{array} \right. \\ 指数分布\\ポアソン分布\\ワイブル分布 \left\{ \begin{array}{l} ２母数 \\ ３母数 \end{array} \right.\end{array} \right. \]
5.1 回帰分析
p78 例題5.1 線形回帰分析（度数なしデータ、正規分布）
（処理手順）
- \(n\) 個の時刻データ \(x_i\) を設定
- ミーンランク（累積確率）\(F(x_i)=\cfrac{i}{n+1}\) を計算
- 正規分布を前提として、上側確率パーセント点 \(y=\Phi^{-1}(F)\)を計算
- 回帰直線を求めるに必要な係数を計算
- 一覧表を作成（表5.1、p166 表A1）
（使用モジュール）
- 正規分布の累積確率に対応するパーセント点：scipy.stats.norm.ppf
- リストの昇順並び替え：sorted
from scipy.stats import norm import pandas as pd # 故障時間の観測データ f_time = [33.5,21.5,19,24,35,27,23,29.5,31,26,28.5,25,27.5,32] f_time = sorted(f_time) # 昇順に並び替え # 累積確率F F_list = [item/(len(f_time)+1) for item in range(len(f_time)+1)] F_list.pop(0) # 上側確率パーセント点 y_value = [norm.ppf(f_value) for f_value in F_list] # 表５．１ 回帰直線を求めるに必要な係数 myu_x = sum(f_time)/len(f_time) myu_y = 0 DF = pd.DataFrame([f_time, F_list, y_value], index =['x','F','y']).T DF['xy'] = DF['x']*DF['y'] DF['x^2'] = DF['x']**2 DF['(x-μ_x)^2'] = (DF['x']-myu_x)**2 DF['(y-μ_y)^2'] = (DF['y']-myu_y)**2 DF.loc['Mean'] = DF.mean() pd.options.display.precision = 3 DF.at['Mean', 'F'] = None DF.at['Mean', 'y'] = myu_y DF
x F y xy x^2 (x-μ_x)^2 (y-μ_y)^2 0 19.000 0.067 -1.501 -28.521 361.000 69.246 2.253 1 21.500 0.133 -1.111 -23.882 462.250 33.889 1.234 2 23.000 0.200 -0.842 -19.357 529.000 18.675 0.708 3 24.000 0.267 -0.623 -14.950 576.000 11.032 0.388 4 25.000 0.333 -0.431 -10.768 625.000 5.389 0.186 5 26.000 0.400 -0.253 -6.587 676.000 1.746 0.064 6 27.000 0.467 -0.084 -2.259 729.000 0.103 0.007 7 27.500 0.533 0.084 2.300 756.250 0.032 0.007 8 28.500 0.600 0.253 7.220 812.250 1.389 0.064 9 29.500 0.667 0.431 12.706 870.250 4.746 0.186 10 31.000 0.733 0.623 19.311 961.000 13.532 0.388 11 32.000 0.800 0.842 26.932 1024.000 21.889 0.708 12 33.500 0.867 1.111 37.211 1122.250 38.175 1.234 13 35.000 0.933 1.501 52.538 1225.000 58.960 2.253 Mean 27.321 NaN 0.000 3.707 766.375 19.915 0.691
p78 例題5.1 回帰係数、期待値μ、標準偏差σ
# コードは上記の続き
# A1 (式5.3)
A1_value = (DF['xy']['Mean']-DF['x']['Mean']*DF['y']['Mean'])/(DF['x^2']['Mean']-DF['x']['Mean']**2)
print(f'回帰係数A1：{round(A1_value,4)}')
# A0 (式5.4)
A0_value = DF['y']['Mean']- A1_value * DF['x']['Mean']
print(f'回帰係数A0：{round(A0_value,3)}')
# 期待値μ (式5.6)
print(f'期待値μ：{round(-A0_value/A1_value,2)}時間')
# 標準偏差σ (式5.6)
print(f'標準偏差σ：{round(1/A1_value,3)}時間')
回帰係数A1：0.1861 回帰係数A0：-5.085 期待値μ：27.32時間 標準偏差σ：5.372時間
p78 例題5.1 図5.3 正規確率分布のプロット
from matplotlib import pyplot as plt
import numpy as np
# ｘとｙの範囲
xmin = 10
xmax = 40
ymin = np.log10(1/6000)
ymax = np.log10(1200)
# ２つのグラフを設定
fig = plt.figure()
ax1 = fig.add_subplot()
ax2 = ax1.twinx()
# 右ｙ軸に対数目盛を設定
ax2.set_yticks([])
_dy = np.array([0.001,0.01,0.05,0.1,0.2,0.3,0.5,0.7,0.8,0.9,0.95,0.99,0.999])
dy = norm.ppf(_dy, loc=0, scale=1) # 対数値に変換
for i in range(0,len(_dy)):
plt.text(xmax+0.5, dy[i], str(_dy[i]*100), ha = 'left', va = 'center')
plt.hlines(dy,xmin,xmax, color='grey',linestyle=":")
# xy値をプロット
ax1.scatter(f_time,y_value,color ='black')
# 回帰直線を描画
x_line = np.arange(xmin,xmax,0.1)
y_line = [A0_value + A1_value*xs for xs in x_line]
ax1.plot(x_line,y_line,label=f'$y={round(A0_value,3)}+{round(A1_value,4)}x$')
# 描画領域
plt.xlim(xmin,xmax)
ax1.set_ylim(ymin, ymax)
ax2.set_ylim(ymin, ymax)
# 凡例等
ax1.legend()
ax1.grid()
ax1.set_xlabel("failure time: $x$") #x軸ラベル
ax1.set_ylabel("percent point: $y$") #y1軸ラベル
ax2.set_ylabel("mean rank: $F$ [%]", labelpad = 40) #y2軸ラベル
plt.title('Normal probability plot')
plt.show()
p80 例題5.2 線形回帰分析（度数なしデータ、対数正規分布）
（処理手順）
- \(n\) 個の時刻データ \(z_i\) を設定
- 対数 \(x=\log_e z\) を計算
- ミーンランク（累積確率）\(F(i)=\cfrac{i}{n+1}\) を計算
- 正規分布を前提として、上側確率パーセント点 \(y=\Phi^{-1}(F)\)を計算
- 回帰直線を求めるに必要な係数を計算
- 一覧表を作成（表5.2、p166 表A2）
import math from scipy.stats import norm import pandas as pd # 故障時間の観測データ z_time = [170,430,230,140,210,285,340,190,395,240,310,260] z_time = sorted(z_time) # 昇順に並び替え x_time = [math.log(z_value) for z_value in z_time] # 対数化 # 累積確率F F_list = [item/(len(x_time)+1) for item in range(len(x_time)+1)] F_list.pop(0) # 上側確率パーセント点 y_value = [norm.ppf(F_value) for F_value in F_list] # 表５．２ pd.options.display.precision = 3 DF = pd.DataFrame([z_time, x_time, F_list, y_value], index =['t','x','F','y']).T # 回帰直線を求めるに必要な係数 myu_x = sum(x_time)/len(x_time) myu_y = 0 DF['xy'] = DF['x']*DF['y'] DF['x^2'] = DF['x']**2 DF['(x-μ_x)^2'] = (DF['x']-myu_x)**2 DF['(y-μ_y)^2'] = (DF['y']-myu_y)**2 DF.loc['Mean'] = DF.mean() pd.options.display.precision = 3 DF.at['Mean', 'F'] = None DF.at['Mean', 'y'] = myu_y DF
t x F y xy x^2 (x-μ_x)^2 (y-μ_y)^2 0 140.000 4.942 0.077 -1.426 -7.047 24.420 3.512e-01 2.034 1 170.000 5.136 0.154 -1.020 -5.239 26.376 1.588e-01 1.041 2 190.000 5.247 0.231 -0.736 -3.863 27.531 8.253e-02 0.542 3 210.000 5.347 0.308 -0.502 -2.686 28.592 3.504e-02 0.252 4 230.000 5.438 0.385 -0.293 -1.595 29.573 9.259e-03 0.086 5 240.000 5.481 0.462 -0.097 -0.529 30.037 2.880e-03 0.009 6 260.000 5.561 0.538 0.097 0.537 30.921 6.958e-04 0.009 7 285.000 5.652 0.615 0.293 1.658 31.951 1.397e-02 0.086 8 310.000 5.737 0.692 0.502 2.882 32.908 4.091e-02 0.252 9 340.000 5.829 0.769 0.736 4.292 33.977 8.681e-02 0.542 10 395.000 5.979 0.846 1.020 6.099 35.747 1.977e-01 1.041 11 430.000 6.064 0.923 1.426 8.647 36.769 2.804e-01 2.034 Mean 266.667 5.534 NaN 0.000 0.263 30.734 1.050e-01 0.661
p80 例題5.2 回帰係数、期待値μ、標準偏差σ
# コードは上記の続き
# A1 (式5.3)
A1_value = (DF['xy']['Mean']-DF['x']['Mean']*DF['y']['Mean'])/(DF['x^2']['Mean']-DF['x']['Mean']**2)
# A0 (式5.4)
A0_value = DF['y']['Mean']- A1_value * DF['x']['Mean']
print(f'回帰係数A0：{round(A0_value,3)}')
print(f'回帰係数A1：{round(A1_value,4)}')
# 期待値μ、標準偏差σ (式5.8)
myu_L = -A0_value/A1_value
sigma_L = 1/A1_value
myu = math.exp(myu_L + sigma_L**2/2)
sigma = math.sqrt(math.exp(2*myu_L+sigma_L**2)*(math.exp(sigma_L**2)-1))
print(f'期待値μ：{round(myu,2)}時間')
print(f'標準偏差σ：{round(sigma,2)}時間')
回帰係数A0：-13.856 回帰係数A1：2.5037 期待値μ：274.26時間 標準偏差σ：114.06時間
p80 例題5.2 図5.4 対数正規確率分布のプロット
from matplotlib import pyplot as plt
import numpy as np
# ｘとｙの範囲
xmin = 4.5
xmax = 7
ymin = np.log10(1/6000)
ymax = np.log10(1200)
# ２つのグラフを設定
fig = plt.figure()
ax1 = fig.add_subplot()
ax2 = ax1.twinx()
ax3 = ax1.twiny()
# 右ｙ軸に対数目盛を設定
ax2.set_yticks([])
_dy = np.array([0.001,0.01,0.05,0.1,0.2,0.3,0.5,0.7,0.8,0.9,0.95,0.99,0.999])
dy = norm.ppf(_dy, loc=0, scale=1) # パーセント値に変換
for i in range(0,len(_dy)):
plt.text(xmax+0.05, dy[i], str(_dy[i]*100), ha = 'left', va = 'center')
plt.hlines(dy,xmin,xmax, color='grey',linestyle=":")
# 上ｘ軸に対数目盛を設定
ax3.set_xticks([])
_dz = np.array([100,200,300,500,700,1000])
dz = [math.log(z_value) for z_value in _dz] # 対数値に変換
for j in range(0,len(dz)):
plt.text(dz[j],ymax+0.2, str(_dz[j]), ha = 'left', va = 'center')
# xy値をプロット
ax1.scatter(x_time,y_value,color ='black')
# 回帰直線を描画
x_line = np.arange(xmin,xmax,0.1)
y_line = [A0_value + A1_value*xs for xs in x_line]
ax1.plot(x_line,y_line,label=f'$y={round(A0_value,2)}+{round(A1_value,3)}x$')
# 描画領域
plt.xlim(xmin,xmax)
ax1.set_ylim(ymin, ymax)
ax2.set_ylim(ymin, ymax)
# 目盛り
ax1.grid(axis="x")
ax1.minorticks_on()
# 凡例等
ax1.legend()
ax1.set_xlabel("failure time: $x=\log_e z$") #x軸ラベル
ax1.set_ylabel("percent point: $y$") #y1軸ラベル
ax2.set_ylabel("mean rank: $F$ [%]", labelpad = 40) #y2軸ラベル
ax3.set_xlabel("Observed Time: z [hours]", labelpad = 20) #x2軸ラベル
plt.title('Log-Normal probability plot',size=15)
plt.show()
p81 例題5.3 線形回帰分析（累積度数データ、２母数ワイブル分布）
（処理手順）
- \(n\) 個の時刻データ \(t_i\) 、累積個数データを設定
- 累積個数データ →頻度 \(n_i\) を計算
- 累積個数データ →ミーンランク（累積確率）\(F(i)=\cfrac{i}{n+1}\) を計算
- ２母数ワイブル分布を前提に、上側確率パーセント点 \(y=\log_e \left( \log_e(1-F)^{-1} \right) \) を計算
- 対数 \(x=\log_e t\) を計算
- 回帰直線を求めるに必要な係数を計算
- 一覧表を作成（表5.4、p167 表A3）
import math
from scipy.stats import norm
import pandas as pd
# 観測データ
t_time = [2,3,7,14,25,40,48] # 時刻
cums = [2,4,14,36,64,86,96] # 累積個数
s_num = 100 # 全個数
# 累積度数 →頻度
cumse = [0]+ cums # 冒頭に0を追加
number = [i-j for i,j in zip(cumse[1:],cumse[:-1])] # 頻度
# 累積確率F
F_list = [item/(s_num+1) for item in cumse]
F_list.pop(0)
# 上側確率パーセント点
y_value = [math.log(math.log(1/(1-F_value))) for F_value in F_list]
# x=log t
x_time = [math.log(tvalue) for tvalue in t_time]
# 表５．３～５．４ 回帰直線を求めるに必要な係数
DF = pd.DataFrame([t_time, number,F_list, x_time, y_value], index =['t','n','F','x','y']).T
DF['nx'] = DF['n']*DF['x']
DF['ny'] = DF['n']*DF['y']
DF['nxy'] = DF['nx']*DF['y']
DF['nx^2'] = DF['nx']*DF['x']
nums = DF['n'].sum()
myu_x = DF['nx'].sum()/nums
myu_y = DF['ny'].sum()/nums
DF['n(x-μ_x)^2'] = DF['n']*(DF['x']-myu_x)**2
DF['n(y-μ_y)^2'] = DF['n']*(DF['y']-myu_y)**2
myu_xy = DF['nxy'].sum()/nums
myu_x2 = DF['nx^2'].sum()/nums
# 表示
print(f'μ_x= {round(myu_x,4)}')
print(f'μ_y= {round(myu_y,4)}')
print(f'Σpxy= {round(myu_xy,4)}')
print(f'Σμx^2= {round(myu_x2,4)}')
pd.options.display.precision = 4
DF
μ_x= 3.0323 μ_y= -0.2705 Σpxy= -0.0176 Σμx^2= 9.7272 t n F x y nx ny nxy nx^2 n(x-μ_x)^2 n(y-μ_y)^2 0 2.0 2.0 0.0198 0.6931 -3.9120 1.3863 -7.8240 -5.4232 0.9609 10.9430 26.5213 1 3.0 2.0 0.0396 1.0986 -3.2087 2.1972 -6.4174 -7.0502 2.4139 7.4781 17.2663 2 7.0 10.0 0.1386 1.9459 -1.9024 19.4591 -19.0238 -37.0187 37.8657 11.8017 26.6314 3 14.0 22.0 0.3564 2.6391 -0.8193 58.0593 -18.0249 -47.5689 153.2217 3.4015 6.6271 4 25.0 28.0 0.6337 3.2189 0.0042 90.1285 0.1174 0.3780 290.1125 0.9750 2.1123 5 40.0 22.0 0.8515 3.6889 0.6456 81.1553 14.2025 52.3913 299.3723 9.4850 18.4608 6 48.0 10.0 0.9505 3.8712 1.1005 38.7120 11.0050 42.6027 149.8620 7.0381 18.7957
p81 例題5.3 回帰係数、母数αβ、期待値μ、標準偏差σ
（処理手順）
- 回帰係数 \(A_0, A_1\) を計算
- 回帰係数 →母数 \(\alpha, \beta\) を計算
- 母数 →期待値 \(\mu\)、標準偏差\(\sigma\) を計算
# コードは上記の続き
from scipy.special import gamma
# A1(式5.3)_A0 (式5.4)
A1_value = (myu_xy -myu_x *myu_y)/(myu_x2 -myu_x**2)
A0_value = myu_y -A1_value * myu_x
print(f'回帰係数A0：{round(A0_value,3)}')
print(f'回帰係数A1：{round(A1_value,4)}')
# α、β(式5.10)
alpha = A1_value
beta = math.exp(-A0_value/alpha)
print(f'α：{round(alpha,3)}')
print(f'β：{round(beta,3)}')
# 期待値μ、標準偏差σ (式4.47_4.49)
myu_W = beta *gamma(1+1/alpha)
simga_W = beta *math.sqrt(gamma(1+2/alpha)-(gamma(1+1/alpha))**2)
print(f'期待値μ：{round(myu_W,2)}時間')
print(f'標準偏差σ：{round(simga_W,2)}時間')
回帰係数A0：-4.84 回帰係数A1：1.5071 α：1.507 β：24.822 期待値μ：22.4時間 標準偏差σ：15.14時間
p81 例題5.3 図5.5 ２母数ワイブル確率分布のプロット
from matplotlib import pyplot as plt
import numpy as np
# ｘとｙの範囲
xmin = -2.5
xmax = 4.6
ymin = -7
ymax = 2
# ３つのグラフを設定
fig = plt.figure()
ax1 = fig.add_subplot()
ax2 = ax1.twinx()
ax3 = ax1.twiny()
# 右ｙ軸に対数目盛を設定
ax2.set_yticks([])
_dy = np.array([0.001,0.01,0.02,0.05,0.1,0.2,0.4,0.5,0.7,0.9,0.95,0.99,0.999])
dy = [math.log(math.log(1/(1-dy))) for dy in _dy] # ln ln に変換
for i in range(0,len(_dy)):
plt.text(xmax+0.05, dy[i], str(_dy[i]*100), ha = 'left', va = 'center',size=9)
plt.hlines(dy,xmin,xmax, color='grey',linestyle=":") # 目盛り線を描画
# 上ｘ軸に対数目盛を設定
ax3.set_xticks([])
_dz = np.array([0.1,0.2,0.5,1,2,5,10,20,50,100])
dz = [math.log(z_value) for z_value in _dz] # 対数値に変換
for j in range(0,len(dz)):
plt.text(dz[j],ymax+0.3, str(_dz[j]), ha = 'center', va = 'center', size=7)
plt.vlines(dz,ymin,ymax, color='grey',linestyle=":") # 目盛り線を描画
# xy値を下ｘ、左軸平面にプロット
ax1.scatter(x_time,y_value,color ='black')
# 回帰直線を描画
x_line = np.arange(xmin,xmax,0.1)
y_line = [A0_value + A1_value*xs for xs in x_line]
ax1.plot(x_line,y_line,label=f'$y={round(A0_value,2)}+{round(A1_value,3)}x$')
# 描画領域
plt.xlim(xmin,xmax)
ax1.set_ylim(ymin, ymax)
ax2.set_ylim(ymin, ymax)
# 凡例等
ax1.legend(loc='lower right', fontsize=12)
ax1.set_xlabel(r"Elapsed time: $x=\ln t$") # 下ｘ軸ラベル
ax1.set_ylabel(r"$y=\ln\left(\ln\frac{1}{1-F}\right)$") # 左ｙ軸ラベル
ax2.set_ylabel("mean rank: $F$ [%]", labelpad = 30) # 右ｙ軸ラベル
ax3.set_xlabel("Observed Elapsed Time: t [hours]", labelpad = 20) # 上ｘ軸ラベル
plt.title('Weibull probability plot',size=15)
plt.show()
5.2 最尤法
p85 例題5.4 最尤法（度数データなし、指数分布）
（処理手順）
- 時刻データを設定
- 指数分布の母数（故障率λ）を計算
f_time = [380,60,140,430,230,310,100,160,270,290,190,250,210,340]
lambda_value = len(f_time)/sum(f_time)
print(f'指数分布の故障率λ：{round(lambda_value*100,4)}[%/hour]')
指数分布の故障率λ：0.4167[%/hour]
p86 例題5.5 最尤法（度数データなし、正規分布）
（処理手順）
- 時刻データを設定
- 正規分布の母数（平均\(\mu\)、標準偏差\(\sigma\)を計算）
import math
o_time = [380,60,140,430,230,310,100,160,270,290,190]
myu = sum(o_time)/len(o_time)
sigma = math.sqrt(sum([(o_value-myu)**2 for o_value in o_time])/len(o_time))
print(f'正規分布の平均μ：{round(myu,1)}[hour]')
print(f'正規分布の標準偏差σ：{round(sigma,1)}[hour]')
正規分布の平均μ：232.7[hour] 正規分布の標準偏差σ：110.5[hour]
p87 例題5.6 最尤法（度数データなし、ワイブル分布）
（処理手順）
- 疲労寿命データを設定
- 関数（αから次のαに更新）を設定
- 同関数にて、更新差異が小さくなるまで再帰演算 →\(\alpha,\beta\) を計算
# 例題５．６
import math
N_cycle = [3.66,4.67,5.43,5.98,6.52,7.08,7.37,7.61,8.22,8.63,9.19,9.48,10.35,10.88,11.21]
alpha = 3 # 暫定初期値
# 関数 αから次のα候補を求める
def next_alpha(alpha,num_list):
ln_x = [math.log(num) for num in num_list] # log x のリスト
x_a = [num**alpha for num in num_list] # xのa乗のリスト
term1 = sum([xa*lnx for xa, lnx in zip(x_a,ln_x)]) /sum(x_a)
term2 = 1/len(num_list) *sum(ln_x)
g_a = 1/(term1-term2) # p165 A13式
return (alpha + g_a)/2 # p164 A14式
# 関数 更新差分が収束するまで再帰処理
def recu_func(seed,num_list):
new_alpha = next_alpha(seed,num_list) # α
new_beta = (sum([num**new_alpha for num in num_list])/len(num_list))**(1/new_alpha) # β
diff = new_alpha - seed # α差分
if diff ＆lt;0.000001: # α差分の精度
return new_alpha, new_beta
else:
return recu_func(new_alpha,num_list)
# α、β演算
approx = recu_func(alpha,N_cycle)
print(f'最尤推定値α：{round(approx[0],3)}')
print(f'最尤推定値β：{round(approx[1],3)}')
最尤推定値α：4.079 最尤推定値β：8.563
5.3 分布のχ二乗適合度検定
p90 例題5.7 χ二乗適合度検定（度数データあり、ポアソン分布）
（処理手順）
- 有意水準\(\alpha\)、度数データ（保全数、頻度）を設定
- 度数データ →一覧表
- 度数データ →グラフ化
import pandas as pd
from matplotlib import pyplot as plt
alpha = 0.05 # 有意水準
DF = pd.DataFrame([[1,2,3,4,5,6,7,8]
,[8,11,20,17,17,13,8,6]]
,index=['maintenance','frequency'])
display(DF)
# グラフ化
plt.bar(DF.loc['maintenance'],DF.loc['frequency'])
plt.title('Maintenace Distribution')
plt.xlabel('Number of Maintenance')
plt.ylabel('Frequency')
plt.show()
0 1 2 3 4 5 6 7 maintenance 1 2 3 4 5 6 7 8 frequency 8 11 20 17 17 13 8 6
p90 例題5.7 表5.6 χ二乗検定の計算結果（ポアソン分布）
（処理手順）
- 平均値\(\mu\)を計算
- ポアソン分布を前提として、\(P=\cfrac{\mu^i e^{-\mu}}{i!}\) を計算
- 実測 \(\chi^2\) を計算
# コードは上記の続き
import math
DFT = DF.T
DFT['xy'] = DFT['maintenance']*DFT['frequency']
# 平均値
myu = DFT['xy'].sum()/DFT['frequency'].sum()
print(f'平均値：{myu}')
# P, nP
DFTP =pd.Series([myu **x *math.e**(-myu)/math.factorial(x) for x in DFT['maintenance']],name='P')
DFT = pd.concat([DFT, DFTP],axis=1)
DFT['nP'] = DFT['P']*100
DFT.iat[len(DFT['nP'])-1,4] += 100 -DFT['nP'].sum() # nP 縦合計を100%化
# χ二乗値
chi2 =pd.Series([(y-nP)**2 /nP for y, nP in zip(DFT['frequency'],DFT['nP'])],name='χ^2')
DFT = pd.concat([DFT, chi2],axis=1)
sigma_chi2 = chi2.sum()
print(f'χ二乗値：{round(sigma_chi2,4)}')
DFT
平均値：4.25 χ二乗値：2.217 maintenance frequency xy P nP χ^2 0 1 8 8 0.060623 6.062299 0.619350 1 2 11 22 0.128824 12.882386 0.275056 2 3 20 60 0.182500 18.250047 0.167799 3 4 17 68 0.193907 19.390675 0.294746 4 5 17 85 0.164821 16.482074 0.016275 5 6 13 78 0.116748 11.674802 0.150422 6 7 8 56 0.070883 7.088273 0.117271 7 8 6 48 0.037656 8.169443 0.576108
ｐ90 例題5.7 χ二乗検定、グラフ
（処理手順）
- 自由度\(\nu\)を設定
- 有意水準\(\alpha\)下での理論値 \(\chi^2\) を計算
- 理論値 \(\chi^2\) ＜ 実測 \(\chi^2\) のとき、仮説（ポアソン分布）を棄却
（使用モジュール）
- \(\chi^2\) 分布の残存確率から\(\chi^2\) 値を求める：scipy.stats.chi2.isf
- \(\chi^2\) 値から、確率密度を求める：scipy.stats.chi2.pdf
# コードは上記の続き
import numpy as np
# 自由度ν
m_value = 1 # ポアソン分布の母数は、μの１個のみ
nu = len(DFT['maintenance']) -1 -m_value
# 有意水準α下でのχ二乗値
from scipy.stats import chi2
chi2_theory = chi2.isf(alpha, nu)
print(f'有意水準{alpha*100}%でのχ二乗値：{round(chi2_theory,4)}')
# 判定結果
if chi2_theory ＆lt;= sigma_chi2:
print(f'有意水準{alpha*100}%で、仮説（ポアソン分布に適合）を棄却')
else:
print(f'有意水準{alpha*100}%で、仮説（ポアソン分布に適合）は未棄却')
# # χ二乗分布のグラフ描画用データ
chi2_list = np.linspace(0,20,100)
fx_list = [chi2.pdf(chi2_value, nu) for chi2_value in chi2_list]
# グラフ化
plt.plot(chi2_list,fx_list, label=f'ν={nu}')
plt.fill_between(chi2_list, fx_list, 0, where=(chi2_list>=chi2_theory), facecolor='blue', alpha=0.5)
# 目盛り、凡例
plt.xticks([0, sigma_chi2, chi2_theory])
plt.title(f'Probability Density')
plt.xlabel(f'$ \chi^2 $')
plt.ylabel(f'$ f\ (\chi^2) $')
plt.legend()
plt.xlim(0,)
plt.ylim(0,)
plt.show()
有意水準5.0%でのχ二乗値：12.5916 有意水準5.0%で、仮説（ポアソン分布に適合）は未棄却
演習問題５
演習問題5.1（線形回帰分析、度数データ、正規分布）
（１）前提数値表の作成 p178 解表１（作業手順）
- 保全時間、保全数のデータを設定
- 保全数 →累積保全数
- 累積保全数 →累積確率\(F\)
- 正規分布を前提に、累積確率\(F\) →上側確率パーセント点\(y\)
- 一覧表を作成
（使用モジュール）
- 配列から累積度数配列を作成：itertools.accumulate
import itertools from scipy.stats import norm import pandas as pd # 保全時間の観測データ m_time = [7,7.5,8,8.5,9,9.5,10,10.5,11,11.5] m_num = [4,6,6,15,15,3,14,8,4,5] # DFデータ生成 cnum = list(itertools.accumulate(m_num)) # 累積個数 alln = cnum[-1] # 全個数 p_list = [num/alln for num in m_num] # 相対度数 F_list = [item/(alln+1) for item in cnum] # 累積確率F y_value = [norm.ppf(f_value) for f_value in F_list] # 上側確率パーセント点 DF = pd.DataFrame([m_time,F_list,y_value],index=['保全時間x','累積確率F','y']) pd.options.display.precision = 4 display(DF)
0 1 2 3 4 5 6 7 8 9 保全時間x 7.0000 7.5000 8.0000 8.5000 9.0000 9.5000 10.0000 10.5000 11.0000 11.5000 累積確率F 0.0494 0.1235 0.1975 0.3827 0.5679 0.6049 0.7778 0.8765 0.9259 0.9877 y -1.6509 -1.1579 -0.8505 -0.2984 0.1710 0.2662 0.7647 1.1579 1.4461 2.2462
（２）回帰直線の導出、母数（期待値、標準偏差）
（作業手順）
- 回帰直線の係数\(A_0,A_1\) を算出
- \(A_0,A_1\) →期待値\(\mu\)、標準偏差\(\sigma\)を計算
# コードは上記の続き
# 回帰直線を求めるに必要な係数
myu_x = sum([iv*jv for iv,jv in zip(m_time,m_num)])/alln
myu_y = sum([kv*lv for kv,lv in zip(y_value,m_num)])/alln
sigma_pxy = sum([pv*xv*yv for pv,xv,yv in zip(p_list,m_time,y_value)])
sigma_px2 = sum([pv*xv**2 for pv,xv in zip(p_list,m_time)])
# A1 (式5.3)
A1_value = (sigma_pxy -myu_x *myu_y)/(sigma_px2 -myu_x**2)
print(f'回帰係数A1：{round(A1_value,3)}')
# A0 (式5.4)
A0_value = myu_y - A1_value * myu_x
print(f'回帰係数A0：{round(A0_value,3)}')
# 期待値μ (式5.6)
myu = -A0_value/A1_value
print(f'期待値μ：{round(myu,3)} 時間')
# 標準偏差σ (式5.6)
sigma = 1/A1_value
print(f'標準偏差σ：{round(sigma,3)} 時間')
回帰係数A1：0.796 回帰係数A0：-7.121 期待値μ：8.948 時間 標準偏差σ：1.257 時間
（３）観測データと正規分布のグラフ対比
（作業手順）
- 観測データをグラフ表示
- 正規分布をグラフ表示
# コードは上記の続き
import numpy as np
from matplotlib import pyplot as plt
# 観測データのグラフ化
plt.bar(m_time,m_num,width=0.4, label='Observed Data')
# 正規分布のグラフ化
norm_x = np.linspace(0, 15, 101) # 区間を100等分する101点
norm_y = 50 *np.exp(-(norm_x-myu)**2 / (2*sigma**2)) / (np.sqrt(2 * np.pi)*sigma)
plt.plot(norm_x, norm_y, color='black', linestyle ='--', label='Normal Distribution')
# 凡例
plt.xlim(0,)
plt.ylim(0,18)
plt.xlabel('hours')
plt.ylabel('frequency')
plt.legend()
plt.show()
演習問題5.2（線形回帰分析、度数データ、対数正規分布）
（１）前提数値表の作成 p178 解表２
（作業手順）
- 故障発生時間、故障数のデータを設定
- 故障発生時間 →対数時間 \(x\)
- 故障数 →累積個数 →累積確率 \(F\) →上側確率パーセント点 \(y\)
- 一覧表を作成
# 演習問題５．２ import numpy as np import math import itertools from scipy.stats import norm import pandas as pd # 故障発生時間の観測データ z_time = np.linspace(10, 190, 10) # 区間を10等分 f_num = [4,21,17,11,8,4,4,2,3,6] # DFデータ生成 t_time = [math.log(zv) for zv in z_time] cnum = list(itertools.accumulate(f_num)) # 累積個数 alln = cnum[-1] # 全個数 p_list = [num/alln for num in f_num] # 相対度数 F_list = [item/(alln+1) for item in cnum] # 累積確率F y_value = [norm.ppf(f_value) for f_value in F_list] # 上側確率パーセント点 DF = pd.DataFrame([t_time,F_list,y_value],index=['対数時間x','累積確率F','y']) pd.options.display.precision = 4 DF
0 1 2 3 4 5 6 7 8 9 対数時間x 2.3026 3.4012 3.9120 4.2485 4.4998 4.7005 4.8675 5.0106 5.1358 5.2470 累積確率F 0.0494 0.3086 0.5185 0.6543 0.7531 0.8025 0.8519 0.8765 0.9136 0.9877 y -1.6509 -0.4997 0.0464 0.3970 0.6842 0.8505 1.0444 1.1579 1.3631 2.2462
（２）回帰直線の導出、母数（期待値、標準偏差）
（作業手順）
- 回帰直線の係数\(A_0,A_1\) を算出
- \(A_0,A_1\) →期待値\(\mu\)、標準偏差\(\sigma\)
# コードは上記の続き
# 回帰直線を求めるに必要な係数
myu_x = sum([iv*jv for iv,jv in zip(t_time,f_num)])/alln
myu_y = sum([kv*lv for kv,lv in zip(y_value,f_num)])/alln
sigma_pxy = sum([pv*xv*yv for pv,xv,yv in zip(p_list,t_time,y_value)])
sigma_px2 = sum([pv*xv**2 for pv,xv in zip(p_list,t_time)])
# A1 (式5.3)
A1_value = (sigma_pxy -myu_x *myu_y)/(sigma_px2 -myu_x**2)
# A0 (式5.4)
A0_value = myu_y - A1_value * myu_x
print(f'回帰係数A0：{round(A0_value,3)}')
print(f'回帰係数A1：{round(A1_value,4)}')
# 期待値μ、標準偏差σ (式5.8)
myu_L = -A0_value/A1_value
sigma_L = 1/A1_value
myu = math.exp(myu_L + sigma_L**2/2)
sigma = math.sqrt(math.exp(2*myu_L+sigma_L**2)*(math.exp(sigma_L**2)-1))
print(f'期待値μ：{round(myu,2)} 分')
print(f'標準偏差σ：{round(sigma,2)} 分')
回帰係数A0：-4.575 回帰係数A1：1.1905 期待値μ：66.39 分 標準偏差σ：67.21 分
（３）観測データと対数正規分布のグラフ対比
（作業手順）
- 観測データをグラフ表示
- 関数（対数正規分布の確率密度）を設定
- 同関数に基づき、対数正規分布をグラフ作成
# コードは上記の続き
import numpy as np
from matplotlib import pyplot as plt
# 観測データのグラフ化
plt.bar(z_time,f_num,width=10, label='Observed Data')
# 関数 対数正規分布の確率密度 f(z)
def log_norm(z_value, myu, std):
return 1/(std * z_value) * norm.pdf((math.log(z_value)-myu)/std)
# 対数正規分布のグラフ化
norm_x = np.linspace(0.1, 200, 1001) # 区間を1000等分
norm_y = []
for z_value in norm_x:
norm_y.append(1500*log_norm(z_value, myu_L, sigma_L))
plt.plot(norm_x, norm_y, color='black', linestyle ='--', label='Log-Normal Distribution')
# 凡例
plt.xlim(0,)
plt.ylim(0,30)
plt.xlabel('minutes')
plt.ylabel('frequency')
plt.legend()
plt.show()
演習問題5.3（線形回帰分析、度数データなし、ワイブル分布）
（１）前提数値表の作成 p179 解表３
（作業手順）
- 疲労寿命データを設定
- 疲労寿命 →累積確率\(F\) →上側確率パーセント点\(y\)
- 一覧表を作成
import math
import pandas as pd
# 疲労寿命データ
N_cycle = [10.35,10.88,11.21,12.66,12.67,13.43,13.98,14.52,15.08,16.37,17.61,18.22,18.63,19.19,19.48]
N_cycle = sorted(N_cycle) # 昇順に並び替え
x_time = [math.log(N_value) for N_value in N_cycle] # 対数化
# 累積確率F(ミーンランク)
F_list = [item/(len(N_cycle)+1) for item in range(len(N_cycle)+1)]
F_list.pop(0)
# 上側確率パーセント点
y_value = [math.log(math.log(1/(1-F_value))) for F_value in F_list]
# 解表３
pd.options.display.precision = 3
DF = pd.DataFrame([N_cycle, x_time, F_list, y_value], index =['実測値','対数値x','累積確率F','y=ln{ln(1-F)^(-1)}']).T
DF
実測値 対数値x 累積確率F y=ln{ln(1-F)^(-1)} 0 10.35 2.337 0.062 -2.740 1 10.88 2.387 0.125 -2.013 2 11.21 2.417 0.188 -1.572 3 12.66 2.538 0.250 -1.246 4 12.67 2.539 0.312 -0.982 5 13.43 2.597 0.375 -0.755 6 13.98 2.638 0.438 -0.553 7 14.52 2.676 0.500 -0.367 8 15.08 2.713 0.562 -0.190 9 16.37 2.795 0.625 -0.019 10 17.61 2.868 0.688 0.151 11 18.22 2.903 0.750 0.327 12 18.63 2.925 0.812 0.515 13 19.19 2.954 0.875 0.732 14 19.48 2.969 0.938 1.020
（２）回帰直線の導出、母数（期待値、標準偏差）
（作業手順）
- 回帰直線を導出
- 母数（\(\alpha,\beta\)）を計算
- 期待値\(\theta\)、標準偏差\(\sigma\)を計算
# コードは上記の続き
from scipy.special import gamma
# 回帰直線を求めるに必要な係数
myu_x = DF['対数値x'].mean()
myu_y = DF['y=ln{ln(1-F)^(-1)}'].mean()
sum_pxy = sum(DF['対数値x']*DF['y=ln{ln(1-F)^(-1)}']/len(N_cycle))
sum_px2 = sum(DF['対数値x']**2 /len(N_cycle))
print(f'μ_x= {round(myu_x,4)}')
print(f'μ_y= {round(myu_y,4)}')
print(f'Σ_pxy= {round(sum_pxy,4)}')
print(f'Σ_px^2= {round(sum_px2,4)}')
# A1(式5.3)_A0 (式5.4)
A1_value = (sum_pxy -myu_x *myu_y)/(sum_px2 -myu_x**2)
A0_value = myu_y -A1_value * myu_x
print(f'回帰直線：y = {round(A0_value,4)} + {round(A1_value,4)} x')
# α、β(式5.10)
alpha = A1_value
beta = math.exp(-A0_value/alpha)
print(f'α：{round(alpha,4)}')
print(f'β：{round(beta,4)}')
# 期待値Θ、標準偏差σ (式4.47_4.49)
myu_W = beta *gamma(1+1/alpha)
simga_W = beta *math.sqrt(gamma(1+2/alpha)-(gamma(1+1/alpha))**2)
print(f'期待値Θ：{round(myu_W,2)}[kcycle]')
print(f'標準偏差σ：{round(simga_W,2)}[kcycle]')
μ_x= 2.6838 μ_y= -0.5128 Σ_pxy= -1.1704 Σ_px^2= 7.2456 回帰直線：y = -13.4801 + 4.8316 x α：4.8316 β：16.2805 期待値Θ：14.92[kcycle] 標準偏差σ：3.53[kcycle]
（３）観測データとワイブル分布のグラフ対比
（作業手順）
- 観測データをグラフ表示
- 関数（ワイブル分布の確率密度）を設定
- 同関数に基づき、ワイブル分布をグラフ表示
# コードは上記の続き
import numpy as np
import math
from matplotlib import pyplot as plt
# プロット図
level_line = [1/20,1/20,1/20,1/20,1/20,1/20,1/20,1/20,1/20,1/20,1/20,1/20,1/20,1/20,1/20]
plt.scatter(N_cycle, level_line,color='orange')
# ワイブル
t_list = np.arange(0.01,max(N_cycle),0.01) # 設定データ
def Weibull(t_value, alpha, beta): # 関数 ワイブル分布の確率密度 f(t)
return alpha/beta * (t_value/beta)**(alpha-1) * math.exp(-(t_value/beta)**alpha)
ft_list = []
for t_value in t_list:
ft_list.append(Weibull(t_value,alpha,beta))
plt.plot(t_list,ft_list,label=f'α={round(alpha,4)}, β={round(beta,4)}')
# 期待値Θ
plt.vlines(myu_W,0,Weibull(myu_W,alpha,beta),color='black',linestyles='--',label=f'Expected value = {round(myu_W,2)}')
# 凡例
plt.xlim(0,)
plt.ylim(0,)
plt.legend()
plt.xlabel('Fatigue Life [kcycle]')
plt.ylabel('f (kcycle)')
plt.title('Probability Density [Weibull Dist.]')
plt.show()
演習問題5.4（最尤法、度数データなし、正規分布）
（作業手順）- 故障時刻データを設定
- 平均\(\mu\) を計算
- 標準偏差\(\sigma\) を計算
import math
f_time = [33.5,21.5,19,24,35,27,23,29.5,31,26,28.5,25,27.5,32]
myu = sum(f_time)/len(f_time)
sigma = math.sqrt(sum([(o_value-myu)**2 for o_value in f_time])/len(f_time))
print(f'正規分布の平均μ：{round(myu,2)}[hour]')
print(f'正規分布の標準偏差σ：{round(sigma,3)}[hour]')
正規分布の平均μ：27.32[hour] 正規分布の標準偏差σ：4.463[hour]
演習問題5.5（最尤法、累積度数データ、対数正規分布）
（１）対数の最尤推定値\(\mu_L,\sigma_L\)、真数の母数\(\mu,\sigma\)（作業手順）
- 経過時間と累積故障数のデータを設定
- 累積故障数データ →頻度\(n\)
- 経過時間\(z\) →対数時間x
- 対数平均の最尤推定値\(\mu_L\)、対数標準偏差の最尤推定値\(\sigma_L\)を計算
- 真数の平均\(\mu\)、真数の標準偏差\(\sigma\)を計算
import numpy as np
import math
import pandas as pd
z_time = [16,18,22,26,29,32,34,37]
f_cums = [1,2,7,18,30,42,48,50]
# 度数
cumse = [0]+ f_cums # 冒頭に0を追加
f_num = [i-j for i,j in zip(cumse[1:],cumse[:-1])]
# DFデータ生成
t_time = [math.log(zv) for zv in z_time]
alln = f_cums[-1] # 全個数
DF = pd.DataFrame([z_time,t_time,f_num],index=['実測値z','対数時間x','度数n']).T
DF['xn'] = DF['対数時間x']*DF['度数n']
# 対数平均の最尤推定値
myu_L = DF['xn'].sum()/DF['度数n'].sum() # 式5.18の度数n反映版
print(f'対数の平均の最尤推定値：{round(myu_L,3)}')
# 対数標準偏差の最尤推定値
DF['x-μ'] = DF['対数時間x'] -myu_L
DF['n(x-μ)^2'] = DF['度数n'] *DF['x-μ']**2
sigma_L = math.sqrt(1/alln *sum(DF['n(x-μ)^2'])) # 式5.19の度数n反映版
print(f'対数の標準偏差の最尤推定値：{round(sigma_L,4)}')
# 真数の平均μ
myu = math.exp(myu_L + sigma_L**2 /2) # 式4.38
print(f'真数の平均：{round(myu,2)}[時間]')
# 真数の標準偏差σ
myu = math.sqrt(math.exp(2*myu_L +sigma_L**2 ) *(math.exp(sigma_L**2)-1)) # 式4.40
print(f'真数の標準偏差：{round(myu,3)}[時間]')
pd.options.display.precision = 4
DF
対数の平均の最尤推定値：3.347 対数の標準偏差の最尤推定値：0.1708 真数の平均：28.83[時間] 真数の標準偏差：4.961[時間] 実測値z 対数時間x 度数n xn x-μ n(x-μ)^2 0 16.0 2.7726 1.0 2.7726 -0.5741 0.3296 1 18.0 2.8904 1.0 2.8904 -0.4563 0.2082 2 22.0 3.0910 5.0 15.4552 -0.2556 0.3267 3 26.0 3.2581 11.0 35.8391 -0.0886 0.0863 4 29.0 3.3673 12.0 40.4075 0.0206 0.0051 5 32.0 3.4657 12.0 41.5888 0.1191 0.1701 6 34.0 3.5264 6.0 21.1582 0.1797 0.1937 7 37.0 3.6109 2.0 7.2218 0.2642 0.1397
（２）観測データと対数正規分布のグラフ対比
（作業手順）
- 観測データをグラフ表示
- 関数（対数正規分布の確率密度）を設定
- 同関数に基づき、対数正規分布をグラフ表示
# コードは上記の続き
import numpy as np
from scipy.stats import norm
from matplotlib import pyplot as plt
# 描画範囲
x_up = 60
# 観測データのグラフ化
plt.bar(z_time,f_num, width =3, label='Observed Data')
# 関数 対数正規分布の確率密度 f(z)
def log_norm(z_value, myu, std):
return 1/(std * z_value) * norm.pdf((math.log(z_value)-myu)/std)
# 対数正規分布のグラフ化
norm_x = np.linspace(0.1, x_up, 101) # 区間を100等分
norm_y = []
for z_value in norm_x:
norm_y.append(200*log_norm(z_value, myu_L, sigma_L))
plt.plot(norm_x, norm_y, color='black', linestyle ='--', label='Log-Normal Distribution')
# 凡例
plt.xlim(0, x_up)
plt.ylim(0,)
plt.xlabel('hours')
plt.ylabel('frequency')
plt.legend()
plt.show()
演習問題5.6（最尤法、度数データなし、２母数ワイブル分布）
（作業手順）- 故障時刻データを設定
- 関数（\(\alpha\)から次の\(\alpha\)に更新）を設定
- 同関数にて、更新差異が小さくなるまで再帰演算 →\(\alpha\) 値
- \(\alpha\) 値 →\(\beta\)値
- \(\alpha,\beta\) 値 →平均\(\mu\)、標準偏差\(\sigma\)
# 演習問題５．６
import math
from scipy.special import gamma
# 設定データ
t_time = [27,24.5,29,31.5,20,25.5,30,28,33,23]
t_time.sort()
alpha = 8 # 暫定初期値
# 関数 αから次のα候補を求める
def next_alpha(alpha,num_list):
ln_x = [math.log(num) for num in num_list] # log x のリスト
x_a = [num**alpha for num in num_list] # xのa乗のリスト
term1 = sum([xa*lnx for xa, lnx in zip(x_a,ln_x)]) /sum(x_a)
term2 = 1/len(num_list) *sum(ln_x)
g_a = 1/(term1 -term2) # p165 A13式
return (alpha + g_a)/2 # p164 A10式
# 関数 更新差分が収束するまで再帰処理
def recu_func(seed,num_list):
new_alpha = next_alpha(seed,num_list) # α
diff = new_alpha - seed # α差分
if diff <0.0001: # α差分の精度
return new_alpha
else:
return recu_func(new_alpha,num_list)
# α、β演算
h_alpha = recu_func(alpha,t_time)
h_beta = (sum([tv**h_alpha for tv in t_time])/len(t_time))**(1/h_alpha)
print(f'最尤推定値α：{round(h_alpha,3)}')
print(f'最尤推定値β：{round(h_beta,3)}')
# 平均μ
myu = h_beta * gamma(1+1/h_alpha) # 式4.47
print(f'期待値μ＝{round(myu,1)}[時間]')
# 標準偏差σ
std = h_beta * math.sqrt(gamma(1+2/h_alpha)-(gamma(1+1/h_alpha))**2) # 式4.49
print(f'標準偏差σ＝{round(std,2)}[時間]')
最尤推定値α：8.436 最尤推定値β：28.777 期待値μ＝27.2[時間] 標準偏差σ＝3.84[時間]
演習問題5.7（χ二乗適合度検定、度数データ、正規分布）
（１）表5.10 保全時間分布
（作業手順）
- 保全時間の階級境界値、保全数\(y\)を設定
- 階級境界値 →階級値\(x\)を計算
- 一覧表
import pandas as pd import numpy as np # 設定データ m_time = np.arange(10, 15.5, 0.5) # 階級 m_num = [3,6,7,12,15,5,13,8,6,5] alpha = 0.05 # 有意水準 # DF mid_time = [(i+j)/2 for i,j in zip(m_time[1:],m_time[:-1])] # 階級値 DF = pd.DataFrame([mid_time,m_num],index=['保全時間x','保全数y']) DF
0 1 2 3 4 5 6 7 8 9 保全時間x 10.25 10.75 11.25 11.75 12.25 12.75 13.25 13.75 14.25 14.75 保全数y 3.00 6.00 7.00 12.00 15.00 5.00 13.00 8.00 6.00 5.00
（作業手順）
- 平均、分散を計算
- 階級の確率\(P, nP_i \) を計算
- \(\chi^2\) 値を計算
- 一覧表
# コードは上記の続き
import math
from scipy.stats import chi2
from scipy.stats import norm
DFT = DF.T
# 平均値
DFT['xy'] = DFT['保全時間x']*DFT['保全数y']
myu = DFT['xy'].sum()/DFT['保全数y'].sum()
print(f'平均値：{myu}')
# 分散
DFT['y(x-μ)^2'] = DFT['保全数y'] *(DFT['保全時間x']-myu)**2
var = DFT['y(x-μ)^2'].sum()/DFT['保全数y'].sum()
print(f'分散：{round(var,4)}')
print(f'標準偏差：{round(math.sqrt(var),4)}')
# 階級の確率 P_i
t_list = [(mv-myu)/math.sqrt(var) for mv in m_time] # t値
P_list = [norm.cdf(tv) for tv in t_list] # 階級下限、上限のt値に対応する累積確率
Pi_list = [k-l for k,l in zip(P_list[1:],P_list[:-1])] # 階級幅に対応する確率
Pi_list[0] =P_list[1] # 初端は下限までの全確率
Pi_list[-1] = 1-P_list[-2] # 終端は上限までの全確率
DFT['P_i'] = Pi_list
# nP_i
DFT['nP_i'] = sum(m_num) *DFT['P_i']
# χ二乗値
DFT['χ^2'] = (DFT['保全数y'] -DFT['nP_i'])**2 /DFT['nP_i']
chi2_sum = DFT['χ^2'].sum()
print(f'データのχ二乗値：{round(chi2_sum,4)}')
# p180 解表４
DFT
平均値：12.55 分散：1.4475 標準偏差：1.2031 データのχ二乗値：6.9001 保全時間x 保全数y xy y(x-μ)^2 P_i nP_i χ^2 0 10.25 3.0 30.75 15.87 0.044200 3.535977 0.081242 1 10.75 6.0 64.50 19.44 0.054618 4.369462 0.608463 2 11.25 7.0 78.75 11.83 0.092587 7.406967 0.022360 3 11.75 12.0 141.00 7.68 0.132379 10.590285 0.187653 4 12.25 15.0 183.75 1.35 0.159642 12.771333 0.388914 5 12.75 5.0 63.75 0.20 0.162383 12.990614 4.915081 6 13.25 13.0 172.25 6.37 0.139315 11.145211 0.308675 7 13.75 8.0 110.00 11.52 0.100814 8.065090 0.000525 8 14.25 6.0 85.50 17.34 0.061531 4.922509 0.235853 9 14.75 5.0 73.75 24.20 0.052532 4.202553 0.151318
（処理手順）
- 自由度\(\nu\)を設定
- 有意水準\(\alpha\) 下での理論値 \(\chi^2\) を計算
- 理論値 \(\chi^2\) ＜ 実測 \(\chi^2\) のとき、仮説（正規分布）を棄却
# コードは上記の続き
# 有意水準α下でのχ二乗値
m_value = 2 # 正規分布の母数は、２個（μ、σ）
nu_value = len(m_num) -1 -m_value # 自由度
chi2_theory = chi2.isf(alpha, nu_value)
print(f'有意水準{alpha}のχ二乗値：{round(chi2_theory,4)}')
# 判定結果
if chi2_theory <= chi2_sum:
print(f'有意水準{alpha*100}%で、仮説（正規分布に適合）を棄却')
else:
print(f'有意水準{alpha*100}%で、仮説（正規分布に適合）は未棄却')
# データ値のグラフ
from matplotlib import pyplot as plt
# 正規分布のグラフ描画用データ
norm_list = np.linspace(0,20,100)
fx_list = [50*norm.pdf(x_value,myu,var) for x_value in norm_list]
# 正規分布のグラフ
plt.plot(norm_list,fx_list, color='black', linestyle ='--' ,label=f'μ={myu}\nσ={round(math.sqrt(var),2)}')
plt.bar(DF.loc['保全時間x'],DF.loc['保全数y'], width =0.4)
plt.title('Maintenace Distribution')
plt.xlabel('Time')
plt.ylabel('Numbers')
plt.legend()
plt.show()
有意水準0.05のχ二乗値：14.0671 有意水準5.0%で、仮説（正規分布に適合）は未棄却
