「入門 信頼性工学」第１章～第２章
福井泰好「入門 信頼性工学（第２版）」森北出版 の読後メモ。
設例について、Python ［Google Colaboratory］による演算処理例を示す。
１章 信頼性工学の概要
p5 図1.4 部品数と部品の信頼度Rがアイテムの信頼度におよぼす影響
（作業手順）
- 部品数を 1～\(10^6\) 個まで設定（グラフ横軸）
- 部品の信頼度\(R\)を複数設定
- アイテム（直列システム全体）の信頼度を計算（グラフ縦軸）
- グラフ化
import numpy as np
import matplotlib.pyplot as plt
start = 0
end = 10**6
X_list = np.arange(start, end, 1) #部品数
R_list = [0.9, 0.99, 0.999, 0.9999, 0.99999] # 部品の信頼度
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.grid(color='gray')
for rel in range(len(R_list)):
Y_list = R_list[rel]**X_list # システム全体の信頼度
ax.plot(X_list, Y_list, label=fr'R={R_list[rel]}')
plt.xscale('log') #横軸を対数目盛り
plt.xlabel('Number of parts')
plt.ylabel('Overall system reliability')
plt.xlim(1,) # 部品点数は１以上
plt.ylim(0,)
plt.legend()
plt.show()
- 部品点数が増加するにつれて、直列システムの信頼度は減少する。
- 同減少度合いは、部品の信頼度が小さいほど、強い。
p11 例題1.1 信頼度R
（作業手順）
- 故障時刻と信頼度を求めたい時刻（運用時刻）のリストを設定
- 特定の運用時刻における残存数を計算する関数を設定
- 各残存数から各信頼度を計算
- 一覧表化
import pandas as pd
Num = 3 # アイテム数
f_time = [6,12,20] # 故障時刻
r_time = [5,10,25] # 運用時刻
# 関数 残存数
def residual(f_time,snap):
available = Num
for failure in f_time:
if failure <= snap:
available -= 1
return available
# 残存数と信頼度の計算
res_list = [] # 残存数のリスト
R_list = [] # 信頼度R のリスト
for snap in r_time:
res = residual(f_time,snap)
res_list.append(res)
R_list.append(res/Num*100)
# DF化
DF = pd.DataFrame([r_time, '', res_list, R_list]).T
DF.columns = ['hours','start','available','Reliability[％]']
DF['start'] = Num
pd.options.display.precision = 1
DF
hours start available Reliability[％] 0 5.0 3 3.0 100.0 1 10.0 3 2.0 66.7 2 25.0 3 0.0 0.0
- p14 演習問題1.1 も同様のプログラムで計算できる。
p12 例題1.2 保全度M
（作業手順）
- 保守時間のリストを設定
- 保守時間が基準値未満のものを抽出
- 基準値未満のものの個数と全体個数から保全度を計算
m_time = [30,39,42,44,46,47,50,53,56,58,61,63,64,65,66,67,69,73,76,80,82,85,88,93]
threshold = 70 # 基準指定
don_list = [item for item in m_time if item < threshold] # 基準値未満の抽出
Maintenance = len(don_list)/len(m_time) # 保全度の計算
print(f'保全度：{Maintenance:.2%}')
保全度：70.83%
- p14 演習問題1.2 も同様のプログラムで計算できる。
p13 例題1.3 故障率λ
（作業手順）
- 故障時刻（start）と故障率を求めたい時刻（運用時刻）のリストを設定
- 特定の運用時刻における残存数（available）を計算
- 残存数リストの要素差分を取得し新規故障数（failure）とする
- 新規故障数（failure）を運用時刻の残存数（available）で除して新規故障割合を求める
- 新規故障割合を運用時刻差（end - start）で除して故障率（\(\lambda\)）を求める
- 一覧表化
import pandas as pd
Num = 4 # アイテム数
f_time = [3,6,9,13] # 故障時刻
r_time = [5,10,15,20] # 運用時刻
# 関数 残存数
def residual(f_time,snap):
available = Num
for failure in f_time:
if failure <= snap:
available -= 1
return available
# 残存数の計算
res_list = [] # 残存数のリスト
for snap in r_time:
res = residual(f_time,snap)
res_list.append(res)
# 新規故障リスト
res_list.insert(0,Num) # 冒頭に当初アイテム数を追加
dif = [i-j for i, j in zip(res_list[:-1], res_list[1:])]
# 開始時刻リスト
s_time = [rs-5 for rs in r_time]
# DF化
DF = pd.DataFrame([s_time, r_time, res_list[:-1], dif]).T
DF.columns = ['start','end', 'available', 'failure']
# 故障率λを計算
DF['λ[%/time]'] = (DF['failure']/DF['available'])/(DF['end']-DF['start'])*100
pd.options.display.precision = 1
DF
start end available failure λ[%/time] 0 0 5 4 1 5.0 1 5 10 3 2 13.3 2 10 15 1 1 20.0 3 15 20 0 0 NaN
- p14 演習問題1.3 も同様のプログラムで計算できる。
p14 演習問題1.4 アベイラビリティA
（作業手順）
- 稼働時間のリストを設定
- 非可動時間のリストを計算
- 稼働時間を全体時間（稼働＋非可動）で除してアベイラビリティを求める
up_hours = [20,21,20,19,21,21,22] # アップ時間
down_hours = [24-time for time in up_hours]
available = sum(up_hours)/(sum(up_hours)+sum(down_hours))
print(f'{available:.2%}')
85.71%
- p14 例題1.4 も同様のプログラムで計算できる。
２章 確率と統計の基礎
p18 例題2.2 順列と組み合わせ
全m個中r個故障の集団から、ｎ個を任意抽出したとき、s個故障の確率
\(\cfrac{{}_r C_s \times {}_{m-r} C_{n-s}}{{}_m C_n}\)
import math
all = 10
fail_p = 3
choice = 3
fail_c = 2
prob = math.comb(fail_p,fail_c)*math.comb(all-fail_p,choice-fail_c)/math.comb(all,choice)
print(f'確率：{prob*100}%')
確率：17.5%
p22 例題2.3 独立事象
(作業工程)
- A工場品である確率 ×A工場品の合格率 ＝A工場の合格品である確率 （１）の解
- B工場品である確率 ×B工場品の合格率 ＝B工場の合格品である確率
- A工場の合格品である確率 +B工場の合格品である確率 ＝合格品である確率 （２）の解
- ２個目がB工場品である確率 ×B工場品の合格率 ＝２個目がB工場の合格品である確率
- (1個目が)A工場の合格品である確率 +２個目がB工場の合格品である確率 （３）の解
import pandas as pd
# 個数と不良確率のデータ
a_factory = [100, 0.02]
b_factory = [150, 0.05]
DF = pd.DataFrame([a_factory, b_factory]
,index=['A_factory', 'B_factory']
, columns = ['個数','不良率'] )
# A工場の合格品である確率
prob_a = DF['個数']['A_factory']/DF['個数'].sum() # A工場品である確率
pass_a = 1-DF['不良率']['A_factory'] # A工場品の合格率
answer1 = prob_a * pass_a
print(f'(1)A工場の合格品である確率：{answer1*100}％')
# 合格品である確率
prob_b = DF['個数']['B_factory']/DF['個数'].sum() # B工場品である確率
pass_b = 1-DF['不良率']['B_factory'] # B工場品の合格率
answer2b = prob_b * pass_b # B工場の合格品である確率
answer2 = answer1 + answer2b
print(f'(2)合格品である確率：{answer2*100}％')
# A合格品 →B合格品である確率
prob_b_2 = DF['個数']['B_factory']/(DF['個数'].sum()-1) # 2個目がB工場品である確率
answer3b = prob_b_2 * pass_b # 2個目がB工場の合格品である確率
answer3 = answer1 * answer3b
print(f'(3)A合格品→B合格品である確率：{round(answer3*100,2)}％')
DF
(1)A工場の合格品である確率：39.2％ (2)合格品である確率：96.2％ (3)A合格品→B合格品である確率：22.43％ 個数 不良率 A_factory 100 0.02 B_factory 150 0.05
p24 例題2.4(1) 反復試行の確率
（作業工程）
import math
# 分数を扱う
from decimal import Decimal
from fractions import Fraction
# 関数 x個中y個が不良品である確率
def prob_func(xn,yn,ratio):
nume, deno = Decimal(str(ratio)).as_integer_ratio() # 分子、分母に分解
fault = Fraction(nume,deno) **yn
normal = Fraction(deno-nume,deno)**(xn-yn)
prob = math.comb(xn,yn) *fault *normal
return float(prob)*100
# 例題2.4(1) の設定値
s_num = 5
def_num = 2
def_per = 0.1
# 関数を用いて確率を演算
probability = prob_func(s_num, def_num, def_per)
print(f'(1) 5個中2個が不良品である確率：{round(probability,2)}％')
(1) 5個中2個が不良品である確率：7.29％
p24 例題2.4(2) 反復試行の確率
（作業工程）
- 関数（ｘ個中ｙ個が不良品である確率）を設定
- ｙ個の各個数について上記関数で確率を計算
- 上記確率を合計
# 例題2.4(2) 5個中3個以上が不良品である確率
over_dnum = 3
# 各個について確率を計算して集計
sum_prob = 0
for num in range(over_dnum, s_num+1):
sum_prob += prob_func(s_num, num, def_per)
# 表示
print(f'(2) 5個中3個以上が不良品である確率：{round(sum_prob,3)}％')
(2) 5個中3個以上が不良品である確率：0.856％
p26 表2.1 度数分布
（作業工程）
- 階級と度数を設定
- 累積度数を計算（１階級下の累積度数 ＋度数 ＝現階級の累積度数）
- 一覧表化
# 表２．１
import pandas as pd
from matplotlib import pyplot as plt
import numpy as np
bins_list = np.arange(0, 91, 10) # 階級
freq_list = [0,1,2,3,9,13,12,4,1,0] # 度数
cumu_list = [0] # 累積度数の空リスト
index_list = ['Class','Frequency','Cumulative_Freq']
DF = pd.DataFrame([bins_list, freq_list, cumu_list], index =index_list)
# 累積度数の行を生成
for freq in DF.loc['Frequency']: # 度数を取り出す
cumu = cumu_list[-1]+freq # 累積度数を計算
cumu_list.append(cumu)
cumu_list.pop(0) # 先頭要素を削除
cumu_series = pd.DataFrame(cumu_list).T
# 累積度数の行を結合
DF = pd.concat([DF,cumu_series],axis=0) # 結合
DF.drop('Cumulative_Freq', axis=0, inplace=True) # 不要行の削除
DF.index = index_list # インデックスの再付与
DF
p27 図2.7～2.8 ヒストグラム、累積度数多角形
# コードは上記の続き
plt.bar(DF.loc['Class'], DF.loc['Frequency'], width=9, align="edge", color="lightskyblue", edgecolor="black", label='Frequency')
plt.plot(DF.loc['Class'],DF.loc['Cumulative_Freq'],label='Cumulative_Frequency', color='orange', marker='o')
plt.xlim(0,)
plt.xticks(bins_list)
plt.xlabel('Inspection Time')
plt.ylabel('Frequency')
plt.legend()
plt.grid()
plt.show()
p28 例題2.5 平均、メジアン、モード
（作業工程）
- 平均：numpy.mean で演算
- メジアン：numpy.median で演算
- モード：度数分布表を作成 →最大値を持つ行の階級値
import numpy as np
import pandas as pd
# 設定値
scores = [7230,6560,10800,10990,3200,5520,8680,4650,7810,5320,10080,7260,3600,5640,7050,6650,7180,5840,7500,10310,4630,9020,7620,10040,8860,9200,7420,9350,10400,9530]
width = 1000 # 階級の幅
print(f'平均：{int(np.mean(scores))}')
print(f'メジアン：{int(np.median(scores))}')
# 関数 点数、階級幅 → 階級、度数
def FDDF(score_list, width):
data = np.asarray(score_list)
bins = np.arange(0, data.max()+width+1, width) # 最大値の上限を含む階級値まで階級値を作成
hist = np.histogram(data, bins)[0] # 度数を計算
return bins, hist
# 関数を使用
bins, hist = FDDF(scores,width)
# DF化
indexs = pd.Index([f'{bins[i]}～{bins[i+1]}' for i in range(hist.size)], # インデックス列を生成
name='階級（以上～未満）') # インデックス名称
DF = pd.DataFrame({'階級値': (bins[1:] + bins[:-1]) / 2, # ズレ数列の平均 →中間値
'度数': hist},
index = indexs)
# モード
max_row = DF.loc[DF['度数'].idxmax()] # 最大値を持つ行を取得
print(f'モード：{int(max_row[0])}') # 該当行の階級値
DF
平均：7598 メジアン：7460 モード：7500
p30 例題2.6 散布度（範囲、平均、分散、変動係数）
（作業工程）
- 範囲：最大値numpy.max－最小値numpy.min
- 平均：numpy.mean で演算
- 標準偏差：numpy.std で計算
- 変動係数：標準偏差／平均
import numpy as np
numbers = [38,61,52,14,69,84,55,62,65,27,52,53,76,45,71,56,60,40,46,74]
range = np.max(numbers)-np.min(numbers)
myu = np.mean(numbers)
sigma = np.std(numbers)
eta = sigma/myu
print(f'範囲：{range}\n平均：{myu}\n分散：{round(sigma,3)}\n変動係数：{round(eta,4)}')
範囲：70 平均：55.0 分散：16.565 変動係数：0.3012
p32 例題2.7 ２つのデータを混合
（作業工程）
- 全データ →混合データの平均、標準偏差
- ２データのデータ数・平均・標準偏差 →混合データの平均、標準偏差
import pandas as pd
import numpy as np
# データ値
a_worker = [50.07, 49.96, 50.05, 49.88, 50.09, 50.07, 49.94, 50.03, 49.96, 50.05]
b_worker = [49.98, 50.06, 49.99, 50.01, 50.02, 50.04, 50.03, 49.98, 49.97, 50.02]
all_worker = a_worker + b_worker
# データ数
a_len = len(a_worker)
b_len = len(b_worker)
all_len = len(all_worker)
# 平均
a_mu = np.mean(a_worker)
b_mu = np.mean(b_worker)
all_mu = np.mean(all_worker)
# 標準偏差
a_std = np.std(a_worker)
b_std = np.std(b_worker)
all_std = np.std(all_worker)
# DF化
index_list = ['Num','Mean','Std']
DF = pd.DataFrame([[a_len,b_len,all_len],[a_mu,b_mu,all_mu],[a_std,b_std,all_std]]
,index = index_list
,columns = ['Worker_A','Worker_B','all(definition)'])
# 公式 全平均,全標準偏差
def myu_std(DataFrame):
mux = DF['Worker_A']['Mean']
muy = DF['Worker_B']['Mean']
sx = DF['Worker_A']['Std']
sy = DF['Worker_B']['Std']
m = DF['Worker_A']['Num']
n = DF['Worker_B']['Num']
N = DF['all(definition)']['Num']
mu_all = (m *mux + n* muy)/(m +n) # 全体の平均
s_all = np.sqrt((m *sx**2 +n *sy**2 + (m*n/N)*(mux-muy)**2)/N) # 全体の標準偏差
return N, mu_all, s_all
# 公式での演算
fcalc = pd.DataFrame([item for item in myu_std(DF)], index=index_list, columns= ['all (By_formula)'])
DF = pd.concat([DF,fcalc],axis=1)
DF
Worker_A Worker_B all(definition) all (By_formula) Num 10.000000 10.000000 20.000000 20.000000 Mean 50.010000 50.010000 50.010000 50.010000 Std 0.066332 0.027928 0.050892 0.050892
- 全データから定義に基づき計算した値 all(definition) と公式で求めた値 all (By_formula) は一致している。
p35 例題2.8(1) ２変数の和
（作業工程）
- 確率変数の和が指定値を満足するような確率変数を求める
- 当該確率変数が生起する確率を求める
（使用モジュール）
# 例題２．８（１）
import pandas as pd
# 設定値
x_list = [0,1,2,3]
x_prob = [0.2, 0.4, 0.2, 0.2]
y_list = [0,1,2,3,4]
y_prob = [0.4, 0.2, 0.1, 0.2, 0.1]
# 辞書化
x_dic = dict(zip(x_list, x_prob))
y_dic = dict(zip(y_list, y_prob))
# 関数 足し算
def addition(a_value,b_value):
return a_value + b_value
# 関数 掛け算
def multiple(a_value,b_value):
return a_value * b_value
# 関数 解を与える要素（X,Yの値リスト）を求める
def target(dic_a, dic_b, func, answer): # 引数：辞書２つ、関数、解
target_a = []
target_b = []
for a_ele in dic_a:
for b_ele in dic_b:
if func(a_ele,b_ele) == answer:
target_a.append(a_ele)
target_b.append(b_ele)
return target_a,target_b
# (1) X+Y=7 となる場合 のX,Yを求める
function = addition # 加法を指定
answer = 7 # 解を指定
x_answer = target(x_dic,y_dic,function, answer)[0]
y_answer = target(x_dic,y_dic,function, answer)[1]
# 表示
print(f'X+Y={answer} となる組合せは、{len(x_answer)}通り。')
total_p = 0
for answer in zip(x_answer,y_answer):
prob = round(x_dic[answer[0]]*y_dic[answer[1]],3)
print(f'(X,Y)={answer}の確率：{prob}')
total_p += prob
print(f'合計確率：{total_p}')
X+Y=7 となる組合せは、1通り。 (X,Y)=(3, 4)の確率：0.02 合計確率：0.02
p35 例題2.8(2) ２変数の積
（処理工程）
- 確率変数の積が指定値を満足するような確率変数を求める
- 当該確率変数が生起する確率を求める
# コードは上記の続き
# (2) X*Y=4 となる場合 のX,Yを求める
function = multiple # 乗法を指定
answer = 4 # 解を指定
x_answer = target(x_dic,y_dic,function, answer)[0]
y_answer = target(x_dic,y_dic,function, answer)[1]
# 表示
print(f'X＊Y={answer} となる組合せは、{len(x_answer)}通り。')
total_p = 0
for answer in zip(x_answer,y_answer):
prob = round(x_dic[answer[0]]*y_dic[answer[1]],3)
print(f'(X,Y)={answer}の確率：{prob}')
total_p += prob
print(f'合計確率：{total_p}')
X＊Y=4 となる組合せは、2通り。 (X,Y)=(1, 4)の確率：0.04 (X,Y)=(2, 2)の確率：0.02 合計確率：0.06
p35 例題2.8(3) 確率変数の平均と分散
（処理工程）
- 関数（確率変数の平均）で計算
- 関数（確率変数の分散）で計算
# コードは上記の続き
# 関数 確率変数の平均
def myu(dic):
myu_value = 0
for variable in dic:
myu_value += variable * dic[variable]
return myu_value
# 関数 確率変数の分散
def variance(dic,myu):
var_value = 0
for variable in dic:
var_value += (variable - myu)**2 * dic[variable]
return var_value
# 関数で演算
y_myu = myu(y_dic) # 平均
y_var = variance(y_dic, y_myu) # 分散
# 表示
print(f'Ｙの平均：{y_myu}\nＹの分散：{y_var}')
Ｙの平均：1.4 Ｙの分散：2.04
p38 演習問題2.1 事象確率の積算
（作業工程）
- 特定の故障台数が生ずる確率を計算
- 故障台数を０～３台について合計
- 判定基準との大小比較に従い、結果を表示
import math
# 設定値
total =1000
max_failure = 3
prob = 1/100
rare = 1/100
# 故障台数０～３台の事象確率を積算
probability = 0
for failure in range(max_failure+1):
probability += math.comb(total,failure) *prob **failure *(1-prob)**(total-failure)
# 表示
print(f'生起確率：{round(probability,7)}')
if probability <= rare:
print(f'故障台数{max_failure}台以下との事象は、珍しい')
else:
print(f'故障台数{max_failure}台以下との事象は、珍しくはない')
生起確率：0.0100727 故障台数3台以下との事象は、珍しくはない
p38 演習問題2.2 生起確率
# 演習問題２．２
import math
# 設定値
A_factory = 4
B_factory = 3
select = 3
A_select = 2
# 計算
total = A_factory + B_factory
prob = (math.comb(A_factory,A_select)*math.comb(B_factory, select-A_select))/math.comb(total, select)
# 表示
print(f'確率：{round(prob *100,2)}％')
確率：51.43％
p38 演習問題2.4
import math
# 設定値
A_factory = 150
B_factory = 100
select = 2
at_least = 1
# A工場分が１～２個となる事象確率を積算
total = A_factory + B_factory
probability = 0
for a_num in range(at_least, select+1):
prob = (math.comb(A_factory,a_num)*math.comb(B_factory, select-a_num))/math.comb(total, select)
probability += prob
# 表示
print(f'確率：{round(probability *100,2)}％')
確率：84.1％
p38 演習問題2.5
import pandas as pd
import numpy as np
# データ数
a_len = 45
b_len = 50
# 平均
a_mu = 50.022
b_mu = 49.995
# 標準偏差
a_std = 0.121
b_std = 0.047
# DF化
all_len = a_len + b_len
index_list = ['Num','Mean','Std']
DF = pd.DataFrame([[a_len,b_len,all_len],[a_mu,b_mu,],[a_std,b_std,]]
,index = index_list
,columns = ['旋盤A','旋盤B','All'])
# 公式 全平均,全標準偏差
def myu_std(DataFrame):
mux = DF['旋盤A']['Mean']
muy = DF['旋盤B']['Mean']
sx = DF['旋盤A']['Std']
sy = DF['旋盤B']['Std']
m = DF['旋盤A']['Num']
n = DF['旋盤B']['Num']
N = DF['All']['Num']
mu_all = (m *mux + n* muy)/(m +n) # 全体の平均
s_all = np.sqrt((m *sx**2 +n *sy**2 + (m*n/N)*(mux-muy)**2)/N) # 全体の標準偏差
return N, mu_all, s_all
# 公式での演算
fcalc = pd.DataFrame([item for item in myu_std(DF)], index=index_list, columns= ['all(formula)'])
# 表示
DF = pd.concat([DF,fcalc],axis=1)
DFA = DF.drop('All', axis=1)
DFA.round({'旋盤A':3, '旋盤B':3, 'all(formula)':5})
p38 演習問題2.6(1) ２つのデータの混合
# 設定値
x_list = [0,1,2,3,4]
x_prob = [0.1, 0.2, 0.4, 0.2, 0.1]
y_list = [0,1,2,3]
y_prob = [0.6, 0.2, 0.1, 0.1]
# 辞書化
x_dic = dict(zip(x_list, x_prob))
y_dic = dict(zip(y_list, y_prob))
# 関数 足し算
def addition(a_value,b_value):
return a_value + b_value
# 関数 解を与える要素（X,Yの値リスト）を求める
def target(dic_a, dic_b, func, answer): # 引数：辞書２つ、関数、解
target_a = []
target_b = []
for a_ele in dic_a:
for b_ele in dic_b:
if func(a_ele,b_ele) == answer:
target_a.append(a_ele)
target_b.append(b_ele)
return target_a,target_b
# X+Y=4 となる場合 のX,Yを求める
function = addition # 加法を指定
answer = 4 # 解を指定
x_answer = target(x_dic,y_dic,function, answer)[0]
y_answer = target(x_dic,y_dic,function, answer)[1]
# 表示
print(f'X+Y={answer} となる組合せは、{len(x_answer)}通り。')
total_p = 0
for answer in zip(x_answer,y_answer):
prob = round(x_dic[answer[0]]*y_dic[answer[1]],3)
print(f'(X,Y)={answer}の確率：{prob}')
total_p += prob
print(f'合計確率：{total_p}')
X+Y=4 となる組合せは、4通り。 (X,Y)=(1, 3)の確率：0.02 (X,Y)=(2, 2)の確率：0.04 (X,Y)=(3, 1)の確率：0.04 (X,Y)=(4, 0)の確率：0.06 合計確率：0.16
p38 演習問題2.6(2) ２つのデータの混合（平均）
# コードは上記の続き
# 関数 確率変数の平均
def myu(dic):
myu_value = 0
for variable in dic:
myu_value += variable * dic[variable]
return myu_value
# 関数で各平均を演算
x_myu = myu(x_dic)
y_myu = myu(y_dic)
xy_myu = x_myu * y_myu
# 表示
print(f'Xの平均：{round(x_myu,2)}\nYの平均：{round(y_myu,2)}\nXYの平均：{round(xy_myu,2)}')
Xの平均：2.0 Yの平均：0.7 XYの平均：1.4