4．p％値の算出

4．ｐ％値の算出
１．パラメータが既知の場合（母集団扱い） ................................................................................. 1
1.1
1 確率変数 ............................................................................................................................ 1
1.2
相関のある 2 確率変数......................................................................................................... 2
1.2.1
X , Y が共に正規分布の場合 ......................................................................................... 2
1.2.2
非正規分布の場合.......................................................................................................... 6
1.2.3
算出例............................................................................................................................ 9
２．パラメータが未知の場合（標本扱い）................................................................................... 19
2.1
正規分布 ............................................................................................................................. 19
2.2
非正規分布 ......................................................................................................................... 24
2.2.1
対数正規分布 ............................................................................................................... 24
2.2.2
ワイブル分布 ............................................................................................................... 25
３．計算例...................................................................................................................................... 28
P％値の算出
１．パラメータが既知の場合（母集団扱い）
1.1 1 確率変数
パーセント点（percentile）とは確率 p の値のことで、本書では p % 値と呼ぶことにします。
文献1のように基準強度特性値として 5％値とする。2）基準弾性特性値は 50％値とする。など
と使用されます。50％値は分布形が対称の場合には平均値と同じ意味を持ちます。
5％値あるいは 50％値は、確率分布が既知であれば、クォンタイル関数（逆確率分布関数）か
ら算出できます。
例えば、平均値   10 、標準偏差   2.5 の母集団がある場合、5％値は下のようになります。
正規分布
5％値
xp  qnorm 0.05     
xp  5.888
対数正規分布
パラメータ

  ln 


 
2

ln 1  CV

2
 1  CV 
5％値
  2.272
  0.246
xp  qlnorm  0.05     
xp  6.471
2P ワイブル分布
パラメータ
2




  m  2    m  1  






  m 
 m    m  1  15
m  root CV  


m1
 







 m 


  10.952
 

m1
 

 m 
m  4.542
1
5％値
1
xp  ( ln ( 1  0.05) )
「木質構造設計規準・同解説
m

xp  5.695
許容応力度・許容耐力設計法」
4-1
日本建築学会
P％値の算出
このように、対数正規分布（6.471）＞正規分布（5.888）＞2P ワイブル分布（5.695）となっ
て、平均、標準偏差が同じでも仮定する分布形により 5％値は異なる値となります。
表-1 に変動係数の変化させた場合（0.1～0.5)の対数正規分布、正規分布、2P ワイブル分布の
5%値の比を示します。変動係数によらず対数正規分布の 5％値が最も高くなり、変動係数が 0.1
～0.3 程度ならば正規分布と 2P ワイブル分布の 5％値は近い値となることが分かります。
表- 1 対数正規分布、正規分布、2P ワイブル分布の 5％値の比
変動係数
対数正規分布
正規分布
2P ワイブル分布
0.1
101
100
98
0.2
106
100
96
0.3
117
100
98
0.4
144
100
109
0.5
232
100
155
“R”で p % 値を求めるための関数例は以下の通りです。
pcal<-function(μ,CV,p){#p%値を平均値、変動係数から求める関数
σ<-μ*CV#標準偏差
λ<-log(μ/(1+CV^2)^0.5)#対数正規分布パラメータ
ξ<-(log(1+CV^2))^0.5#対数正規パラメータ
f2p<-function(m){#2P ワイブルパラメータ算出関数
atai<-CV-(gamma(1+2/m)-gamma(1+1/m)^2)^0.5/gamma(1+1/m)
return(atai)
}
m<-uniroot(f2p,lower=1,upper=15)$root#2P ワイブルパラメータ m
η<-μ/gamma(1+1/m)#2P ワイブル分布パラメータη
normp<-qnorm(p,μ,σ)#正規分布 P%値
logp<-qlnorm(p,λ,ξ)#対数正規分布 P%値
weip<-qweibull(p,m,η)#2P ワイブル分布 P%値
c(" 正 規 "=normp," 対 数 正 規 "=logp,"2P ワ イ ブ ル "=weip," 対 数 正 規 / 正 規 "=logp/normp,"2P ワ イ ブ ル / 正 規
"=weip/normp)
}
1.2
相関のある 2 確率変数
相関のある 2 つの確率変数 X , Y の場合です。 X を MOE、 Y を MOR とすると、MOE による
等級区分材の下限値算出などに用いられます。
1.2.1
X , Y が共に正規分布の場合
2 つの確率変数 X , Y を正規分布、相関係数を  とし、 X , Y の平均および標準偏差をそれぞれ
4-2
P％値の算出
 x ,  x ,  y ,  y とすれば、同時確率密度関数(joint probability density function)は 1）式で表され
ます。
f xy ( x, y,  xy ) 
1
2     x   y  1   xy
2
 x   2
 x  x

x
 
  2   xy  
 x
  x 

1
 exp  
2
 2  (1   xy )
  y  x
  
  y
  y  y
  
  y



2
・・・ 1）




 x  10,  x  2.5,  y  40,  y  10 の各確率密度関数を描画すると、図－１のようになります。
（周辺密度関数(marginal density function)と呼ばれます。）
0.2
 y  10,  y  2.5
0.15
0.1
 x  40,  x  10
0.05
0
0
20
40
60
80
図－ 1 周辺密度関数
各確率変数の 5％値をクォンタイル関数で算出すると下記のようになります。
x05  qnorm 0.05   x   x
x05  5.888
y05  qnorm  0.05   y   y
y05  23.551
同時確率密度は図－2 のようになり、4  x  5 の範囲の同時確率密度は図-3 のようになります。
（X＝MOE、Y＝MOR  xy  0.7 、  x  10,  x  2.5,  y  40,  y  10 ）
すなわち、図－3 は MOE による区分材の確率密度を表すことになります。
4-3
P％値の算出
図－ 2 同時確率密度
図－ 3
4  x  5 の同時確率密度
X  x に お け る Y の p% 値 を 求 め る た め に は 、 X  x で の Y の 条 件 付 確 率 密 度 関 数


（conditional density function f x y y x ）を知る必要があります。


同時確率密度関数（ f xy  x, y  ）、Y の確率密度関数（ f y  y  ）、条件付確率密度関数（ f x y y x ）
には 2）式の関係が成り立ちます。
f xy  x, y   f x y  x y   f x  x 
・・・ 2）
4-4
P％値の算出
1）式を変形すると 3）式となりますから 2）式を参照して、 X  x での Y の条件付確率密度関


数（ f x y y x ）は 4）式となります。
f xy  x, y  
1
2     y  1   xy
2
2


y

 

 y    y   xy    x   x    
x
 1 

 
exp    
 
2
2 
 y  1   xy
 


 



・・・ 3）
 1  x   2 
1
x

 exp    
 
2   y
 2   x  
fx y  y x 
1
2     y  1   xy
4）式は平均値が  y   xy 
2
2


y

 

 y    y   xy    x   x    
x
 1 

 
exp    
 
2 
 y  1   xy2
 


 




・・・ 4）
y
2
  x   x  、標準偏差が  y  1   xy の正規分布になっており、Y
x
軸に平行に切り取った分布は正規分布となることを示しています。また、重みなしの直線回帰式
を算出する仮定と同じになります。逆に言えば、重みなしの直線回帰は確率変数 X , Y を共に正規
分布と仮定していることになります。
x  5,10,15 における Y の条件付確率密度は図-4 となります。
0.06
x  10
x5
0.04
x  15
0.02
0
0
20
40
60
80
図－ 4 条件付確率密度関数（ x  5,10,15  x  10,  x  2.5,  y  40,  y  10 ）
条件付確率密度関数（この場合は正規分布）が既知であるため、5％値は正規分布のクォンタイ
ルタイル関数を使って算出でき、 x  5,10,15 について算出すると以下のようになります。
4-5
P％値の算出

y

x
yp  x0  p  qnorm p   y  xy 
yp ( 5  0.05)  14.253
2
  x0   x   y  1  xy 

yp ( 15  0.05)  42.253
yp ( 20  0.05)  56.253
区分を行わない場合（周辺確率密度関数）の 5％値は y p  23.551 でしたが、この値と等しくな
る x を求めると x  8.321 となります。このため、 x  8.321 より小さい区分とした場合には y の
5％値が下がり、逆に x  8.321 より大きい区分とすると y の 5％値は上がることになります。
x0  0  20
60
5％値
40
20
0
0
5
10
15
20
x
図－ 5 y の 5％値（ y p ）
1.2.2 非正規分布の場合
確率変数 X , Y が非正規分布であっても、何らか方法で正規分布扱いができれば 3）式に準じて
同時確率密度関数が定義できそうです。
このために等価正規分布パラメータという概念を導入します。等価正規分布パラメータとは非
正規分布確率変数 X , Y の任意の値において、確率と確率密度が等しい正規分布のパラメータのこ
とを言います。
（理論的には Rosenblatt 変換）
例えば、確率変数 X の x における正規分布のパラメータを  zx ,  zx とすると z x 
x   zx
は標準
 zx
正規分布となりますから 5）、6）式を満たすパラメータを見出すことができるということです。
 x   zx
  zx
確率が等しい： FX  x    

    zx 

4-6
・・・ 5）
P％値の算出
 x   zx 
    zx 
  zx 
確率密度が等しい： f X  x    
ここに
・・・ 6）
Fx  x  , f x ( x) ：累積確率分布関数、確率密度関数

 ,    ：標準正規分布の累積確率分布関数、確率密度関数
6）式、7）式を解いて  zx ,  zx は 8）式となります。
 zx  x   zx  
1
F  x 
x
,  zx 
  1  Fx  x  
fx  x
・・・ 7）
すなわち、非正規分布の同時確率密度関数は 3）式における  x ,  x ,  y ,  y を  zx ,  zx ,  zy ,  zy（等
価正規分布パラメータ）に置き換えることで定義できます。しかし、相関係数（  xy ）と等価正
規分布空間における相関係数（  zxy ）とは同じ値にはなりませんので、  zxy を  xy から算出する
必要があります。
相関係数（  xy ）は平均値  x ,  y のまわりの結合二次積率（共分散 Cov  X , Y  ）を標準偏差  x ,  y
で除したものとして 8）式で定義されます。
（  x ,  y ,  x ,  y は等価正規分布パラメータではありま
せん。
）
 xy 
ここに
   x      y     f  x, y,   dxdy
x
y
xy
zxy
 x  y
・・・ 8）
f xy  x, y,  zxy  ：同時確率密度関数
8）式での定義域は仮定する確率分布形より異なりますので、標準正規化空間（ Z 空間）で書
き換えると 9）式となります。これを解いて換算相関係数  zxy を算出することになりますが、非
線形の方程式を解くことになります。
4-7
P％値の算出
 FX 1   z x    x
 xy    

x
  
  FY

 
 
 
ここに
Fx
1
  z   
y
y
y

  f ZXY  zx , z y ,  zxy  dzx dz y ・・・ 9）


  , Fy 1   ：確率変数 X , Y のクォンタイル関数
f zxy ( z x , z y ,  zxy ) 
zx 
1
1
2    1   zxy
2

1
2
2
 exp  
 zx  2   zxy  zx  z y  z y
2
 2  (1   zxy )




y   zy
x   zx
, zy 
、  zxy ：換算相関係数
 zx
 zy
 zx ,  zx ,  zy ,  zy ：等価正規分布パラメータ
以上のことから、非正規分布の同時確率密度関数は 10）式となります。
f xy  x, y  
1
2     zy  1   zxy
2
2


 zy

 

  x   zx    
 y    zy   zxy 
 zx
 1 

 
exp    
 
2
2 


1


 
zy
zxy


 




・・・ 10）
 1  x   2 
1
zx

 exp    
 
2

2     zy

zx

 
ただし、  zx ,  zx ,  zy ,  zy （等価正規分布パラメータ）、  zxy ：換算相関係数


同様に X  x での Y の条件付確率密度関数（ f x y y x ）は 11）式となります。
fx y  y x 
1
2     zy  1   zxy
2
2


 zy

 

  x   zx    
 y    zy   zxy 
 zx
 1 

   ・・・ 11）
exp    
 
2
2 
 zy  1   zxy
 


 



すなわち、平均値が 0   zy   zxy 
 zy
2
  x   zx  、標準偏差が  0   zy  1   zxy の正規分布
 zx
になっていることが分かります。
4-8
P％値の算出
このため、 X  x での Y の 5％値（ y p ）は、クォンタイル関数（ qnorm  p, 0 ,  0  ）から求めら
れるように感じますが、 0 ,  0 が y の関数ですので、12）式を満たす y p を求める収束計算が必要
となります。（ pnorm 
 ：正規分布の確率分布関数）
pnorm ( y p , 0 ,  0 )  p
1.2.3
・・・ 12）
算出例
X ：平均値  x  100 、変動係数 cv x  0.1 、 Y ：平均値  y  35 、変動係数 cv y  0.2 、相関
係数  xy  0.7 の X の 25％値、50％値、75％値に対する Y の 5％値算出例を示します。
① X ＝正規分布、 Y ＝正規分布
パラメータ
 x  100
cv x  0.1
 x   x  cv x
 x  10
 y  35
cv y  0.2
 y   y  cv y
y  7
xy  0.7
相関係数
周辺密度関数(marginal density function)
f x ( x)  dnorm x   x   x
0.04
fx( x) 0.02
0
60
80
100
120
x
4-9
140
P％値の算出
Y（正規分布）
Y周辺密度関数(marginal density function)
f y ( y)  dnorm  y   y   y
Y周辺密度関数のp％値
ypm ( p)  qnorm  p   y   y
Y条件付確率密度関数（conditional density function)

y

x

y

x
f yxc  y  x0   dnorm y   y  xy 
2
  x0   x   y  1  xy 

Y条件付確率密度関数のp％値
ypc  x0  p  qnorm p   y  xy 
2
  x0   x   y  1  xy 

ypm ( 0.05)  23.486
Y周辺密度関数の5％値
Y条件付確率密度関数の5％値
x25  93.255
ypc  x25  0.05  23.472
x50  100
ypc  x50  0.05  26.777
x75  106.745
ypc  x75  0.05  30.082
確率密度
0.08
fy( y)
0.06
fyxc y  x25
fyxc y  x50
fyxc y  x75
0.04
0.02
0
10
20
30
40
y
4-10
50
60
P％値の算出
② X ＝正規分布、 Y ＝2P ワイブル分布
X（正規分布）
パラメータ
 x  100
 x  10
周辺(marginal)
f x ( x)  dnorm x   x   x
確率密度関数
x   x  3   x   x  3   x
0.04
fx( x) 0.02
0
60
80
100
120
140
x
確率分布関数
Fx ( x)  pnorm x   x   x
クォンタイル関数
qFx  px 
px  pmin if px  0
px  pmax if px  1
atai  qnorm px   x   x
atai
x25  qFx ( 0.25 )
x25  93.255
x50  qFx ( 0.5 )
x50  100
x75  qFx ( 0.75 )
x75  106.745
等価正規変換
 zx ( x) 
  q Fx ( x)  
f x ( x)
 zx ( x)  x   zx ( x)  q F x ( x) 
4-11
P％値の算出
Y（2Pワイブル）
パラメータ
 y  35
 y  35
y  7
y  7
Yの2Pワイブル分布パラメータ算出
mcal      
CV 


m1



m  root  CV 



 

m2
m1
  

m 
 m 
 
m1

 m 
2



 m



m1
 

 m 
 
 
m
 y  mcal   y   y 1
 y  37.801
my  mcal   y   y 2
my  5.792
周辺(marginal)
確率密度関数
my
 y my
my
y 

 
f y ( y)  
  y  exp    

 y
  y 
確率分布関数
 y my
 
F y ( y)  1  exp 
 y  
   
クォンタイル関数
qF y  py 
py  pmin if py  0
py  pmax if py  1
1
atai   ln  1  py 
my
y
atai
等価正規変換
 zy ( y) 
  q F y ( y)  
f y ( y)
 zy ( y)  y   zy ( y) q F y ( y) 
4-12
P％値の算出
Y周辺密度関数のp％値
Y周辺密度関数の5％値
ypm ( p)  qFy ( p)
ypm ( 0.05)  22.635
Y条件付確率密度関数（conditional density function)

f yxc  y  x0  zxy  dnorm y   zy ( y)  zxy 

 zy ( y)
 zx  x0 
2
  x0   zx  x0    zy ( y)  1  zxy 

Y条件付確率密度関数のp％値
ypccal  xyp  zxy  p 

y

x
yp  qnorm p   y  zxy 
2
  xyp   x   y  1  zxy 

atai  root  p  f yxc  yp  xyp  zxy  yp
atai
換算相関係数の算出
関数定義
ｘ、ｙ、ｚ
ｘ、ｚ ｙの定義
x  zx  qFx    zx 
zx ( x)  q F x ( x) 
y  zy  qFy    zy 
zy ( y)  q Fy ( y) 
ｚ空間同時確率密度関数
f z  zx  zy  zxy 
  zx2  2 zxy  zx  zy  zy2  


 exp  
2



2
2  1  zxy
2   1  zxy
 

1
a

a  10
XY空間相関係数定義

xycal  zxy  







a
a a
 x  zx   x   y  zy   y 

 
  f z  zx  zy  zxy d zx d zy
x
y



XY空間同時確率密度関数
f xy  x  y  zxy  f z  zx ( x)  zy ( y)  zxy 
4-13
f x ( x)  f y ( y)
  zx ( x)     zy ( y) 
P％値の算出
算出
精度
初期値
TOL  10
 10
zxy  xy
zxy  0.7
zxy  root  xy  xycal  zxy  zxy
zxy  0.70204
確認
xycal  zxy  0.7
xy  xycal  zxy  0
Y条件付確率密度関数（conditional density function)

f yxc  y  x0  zxy  dnorm y   zy ( y)  zxy 

 zy ( y)
 zx  x0 
2
  x0   zx  x0    zy ( y)  1  zxy 

Y条件付確率密度関数のp％値
ypccal  xyp  zxy  p 

y

x
yp  qnorm p   y  zxy 
2
  xyp   x   y  1  zxy 

atai  root  p  f yxc  yp  xyp  zxy  yp
atai
Y条件付確率密度関数の5％値
ypccal  x25  zxy  0.05  27.501
x25  93.255
ypccal  x50  zxy  0.05  30.929
x50  100
ypccal  x75  zxy  0.05  34.123
x75  106.745
確率密度
0.1
fy( y)


0.05
fyxc y  x50   zxy
fyxc y  x75   zxy
fyxc y  x25   zxy
0
10
20
30
40
y
4-14
50
60
P％値の算出
③ X ＝2P ワイブル分布、 Y ＝2P ワイブル分布
X（2Pワイブル分布 ）
パラメータ
 x  100
 x  10
Xの2Pワイブル分布パラメータ算出
 x  mcal   x   x 1
 x  104.323
mx  mcal   x   x 2
mx  12.093
周辺(marginal)
mx
 x mx
mx
x 

 
f x ( x)  
  x exp     

 x
  x 
確率密度関数
x   x  3   x   x  3   x
0.06
0.04
fx( x)
0.02
0
60
80
100
120
140
x
確率分布関数
 x mx
 
Fx ( x)  1  exp 
 x  
   
クォンタイル関数
qFx  px 
px  pmin if px  0
px  pmax if px  1
1
atai   ln  1  px 
atai
x25  qFx ( 0.25 )
x25  94.11
x50  qFx ( 0.5 )
x50  101.209
x75  qFx ( 0.75 )
x75  107.179
4-15
mx
 x
P％値の算出
等価正規変換
 zx ( x) 
  q Fx ( x)  
 zx ( x)  x   zx ( x)  q F x ( x) 
f x ( x)
Y（2Pワイブル）
パラメータ
 y  35
y  7
Yの2Pワイブル分布パラメータ算出
 y  mcal   y   y 1
 y  37.801
my  mcal   y   y 2
my  5.792
周辺(marginal)
確率密度関数
my
 y my
my
y 

 
f y ( y)  
  y exp    

 y
  y 
確率分布関数
 y my
 
Fy ( y)  1  exp 
 y  
   
クォンタイル関数
qFy  py 
py  pmin if py  0
py  pmax if py  1
1
atai   ln  1  py 
my
 y
atai
等価正規変換
 zy ( y) 
  q F y ( y)  
f y ( y)
Y周辺密度関数のp％値
 zy ( y)  y   zy ( y)  q Fy ( y) 
Y周辺密度関数の5％値
ypm ( p)  qFy ( p)
ypm ( 0.05)  22.635
Y条件付確率密度関数（conditional density function)

f yxc  y  x0  zxy  dnorm y   zy ( y)  zxy 

 zy ( y)
 zx  x0 
4-16
2
  x0   zx  x0    zy ( y)  1  zxy 

P％値の算出
Y条件付確率密度関数のp％値
ypccal  xyp  zxy  p 

y

x
yp  qnorm p   y  zxy 
2
  xyp   x   y  1  zxy 

atai  root  p  f yxc  yp  xyp  zxy   yp 
atai
換算相関係数の算出
算出
精度
初期値
TOL  10
 10
zxy  xy
zxy  0.7
zxy  root  xy  xycal  zxy  zxy
zxy  0.70147
確認
xycal  zxy  0.699999999999315
xy  xycal  zxy   6.846  10
 13
Y条件付確率密度関数の5％値
ypccal  x25  zxy  0.05  27.506
x25  94.11
ypccal  x50  zxy  0.05  30.93
x50  101.209
ypccal  x75  zxy  0.05  34.121
x75  107.179
確率密度
0.1
fy( y)


0.05
fyxc y  x50   zxy
fyxc y  x75   zxy
fyxc y  x25   zxy
0
10
20
30
40
y
4-17
50
60
P％値の算出
X  Y の組み合わせとして、④正規分布－対数正規分布、⑤対数正規分布－対数正規分布、⑥
対数正規分布－2Ｐワイブル分布を付け加えて、算出結果を表－1、図に示します。
①正規分布－正規分布の場合が最も低い値となります。当然ながら Y の分布形の違いが Y の
5％値に影響し、2Ｐワイブル分布＞対数正規分布＞正規分布の順になっています。なお、換算相
関係数の差は少数点以下 3 桁目が異なるだけの結果となりました。
表- 2 5％値の算出結果
項目
X
Y
相関係数
X
Y
相関係数
X
Y
相関係数
X
Y
相関係数
X
Y
相関係数
X
Y
相関係数
36
X
分布形
正規分布
正規分布
0.7000
正規分布
2Pワイブル分布
0.7020
2Pワイブル分布
2Pワイブル分布
0.7015
正規分布
対数正規分布
0.7069
対数正規分布
対数正規分布
0.7038
対数正規分布
2Pワイブル分布
0.7068
34
％
25％
50％
75％
25％
50％
75％
25％
50％
75％
25％
50％
75％
25％
50％
75％
25％
50％
75％
34
34
34
33
31
32
Y
5％値
23.47
26.78
30.08
27.50
30.93
34.12
27.51
30.93
34.12
26.23
29.19
32.52
26.23
29.18
32.50
27.46
30.92
34.14
％値
93.26
100.00
106.75
93.26
100.00
106.75
94.11
101.21
107.18
93.26
100.00
106.75
93.03
99.50
106.43
93.03
99.50
106.43
33
31
31
30
5％値
30
28
27
28
29
29
26
26
28
27
26
24
23
X25%
X50%
X75%
22
20
正規－正規
正規－2P
2P-2P
正規－対数
対数－対数
正規：正規分布、 2P：2Pワイブル分布、対数：対数正規分布
図－ 6 5％値の算出結果
4-18
対数-2P
P％値の算出
２．パラメータが未知の場合（標本扱い）
2.1
正規分布
標本平均値(  es 
n
1
2
   xi   x  )が得られてい
n  1 i 1
1 n
  xi es )、標本不偏標準偏差(  es 
n i 1
るものとすれば、u p を標準正規分布の p % 点として、変量 Tes  es  u p   es は確率変数となりま
す。
例えば、平均値   100 、標準偏差   20 とする母集団から大きさ n1  50 の乱数を発生させ、
es 、  es から、変量 Tes  es  u p   es （5％値： u p  1.644854 ）を計算します。
そして、これを n  1000 回繰り返して真の 5％値（=67.103)より大きい個数（ U n ）と小さい個
数（ Ln )を計算してみます。
（ n  U n  Ln ）
2000
1500
Frequency
2500
3000
3500
Histogram of tl
0
500
1000
> #5%値の実験（正規分布）
> n<-1000#繰り返し回数
> n1<-50#乱数発生回数
> μ<-100#平均値（母集団)
> σ<-20#標準偏差（母集団）
> Ln<-0
> Un<-0
> for (i in 1:n){
+ x<-rnorm(n1,μ,σ)#乱数発生
+ tl[i]<-mean(x)+qnorm(0.05,0,1)*sd(x)#5%値算出
+ if (tl[i]<qnorm(0.05,μ,σ)){Ln<-Ln+1}
+ if (tl[i]>qnorm(0.05,μ,σ)){Un<-Un+1}
+}
> c(n,Ln,Un,Ln+Un)
[1] 1000 463 537 1000
>hist(tl)
40
50
60
70
80
90
tl
変量 Tes  es  u p   es はヒストグラムから分かるように何らかの分布に従い、真の 5％値より
小さくなる個数は 463、大きくなる個数は 537 となることが分かります。
4-19
P％値の算出
発生回数を 100、150、200 回とした場合について表-4 に示します。いずれの結果も標本から求
めた 5％値は真の 5％値より大きい・小さい確率が 50％程度となることが分かります。
表- 3 5％値算出の実験結果（繰り返し数 1000 回）
発生回数
真の 5％値の対して
小さくなる個数
大きいなる個数
50
463
537
100
516
484
150
503
497
200
477
523
強度の 5％値を算出する際に、真の p % 値より変量 Ta   a  u p   a が小さくなる確率（  % ）
を設定し算出する場合があります。
Pr 
 を確率を意味するものとすれば 13）式を満たす k p を使って、14）式から信頼水準が明
確な TL p が計算できることになります。
Pr  es  k p   es    u p     
・・・ 13）
ここに 、  ：母集団（真）の平均、  ：母集団（真）の標準偏差
n
1 n
1
2
 es ：標本平均、  es ：標本不偏標準偏差、 es    xi ,  es 
   xi  es 
n i 1
n  1 i 1
TL p  es  k p   es
・・・ 14）
TL p を正確に言うと、真の p % 値より下回る水準を信頼として、下側信頼水準  % の p % 値と
呼ぶことになります。逆の言い方をして、上側不信頼水準 1    % と呼んでも良いかも知れませ
ん。
この概念は図－7 となります。
4-20
P％値の算出
K・σ ａ
μ a －K・σ ａ
Up・σ
Ｎ（μ、σ）
１－α
μ a －K・σ ａ
図－ 7
μ
Low Tolerance Limits with One-side confidence of Normal Distribution
k p 値の算出は式 13）式を変形して、15）式となることを利用します。
     es  k p   es  

  p
Pr  


 


 Pr    ( es  k p   es )  u p     Pr    u p    es  k p   es 
 n   es    /   n  u p

 Pr  es    U p    k p   es   Pr 
 n  kes 
 es / 


15）式における
n  ( es   ) /   n  u p
 es / 
・・・ 15）
は、自由度  n  1 、非心度λ＝ n  u p の非心 t 分
布に従いますので、個数 n と下側信頼水準  % 、 p % 値を与え、式 16）となるように、収束計算
をして k p 値を算出します。
Pr t  ,    n  k p   
ここに
 k p 
Pr t  ,   
・・・ 16）
n
t  ,   ：非心 t 分布、 Pr t  ,    ：非心 t 分布の確率積分値（＝確率分布関数）
“R”には非心 t 分布の確率分布関数（pt(t 値、自由度、非心度））がありますので、16）式か
ら収束計算することができますが、非心 t 分布の確率分布関数がない場合には、自由度νが奇数の
4-21
P％値の算出
ときは A)式から、自由度νが偶数のときには B)式から求めることができます。
自由度νが奇数：
( 1) / 2
1 exp{ 2 B(1  x 2 ) / 2}
Pr [t ( ,  )]   ( B )  
dx  2   M 2i 1
 0
1  x2
i 1
A
A)
自由度νが偶数：
Pr [t ( ,  )]   ( )  2 
(  2) / 2

i 0
B)
M 2i
M i に関しては
i  0  M 0  A B   ( B )   ( A B )
A


i  1  M 1  B  AM 0 
  ( ) 
2


B
i  2  M 2    AM 1  M 0 
2
j 1
i jMj 
 B  a j  AM j 1  M j  2  ( j  3)
j
1
a2  1, a j 
( j  3, 4,・・・・  2)
( j  2)  a j 1
ここに

 ：標準正規分布の確率密度関数
A
 ：標準正規分布の確率分布関数
t

,B 
  t2

b  b2  a  c
t の初期値： t0 
a
ここで A 

t'

t ' の初期値は t ' 
B
2
1 
u
1


a  1    P , b   1 
 4  2 
 4

2
2
,c    up


であり、
  t '2
1
b  b 2  ac 


a
2
1  U
1


a  1    P , b   1 
 4  2 
 4
“R”で k p を算出する関数の一例は下記の通りです。
#ｋ値算出
4-22

2
2
,c    up

P％値の算出
kcalc<-function(n,p,α){#k 値の算出関数
if (p<=0.5){p<-1-p}
v<-(n-1)#自由度
up<-qnorm(p,0,1)#標準正規分布のｐ%点
λ<-up*n^0.5#非心度
a0<-(1-1/4/v)^2-up/2/v#初期値
b0<-λ*(1-1/4/v)#初期値
c0<-λ^2-up^2#初期値
t<-(b0+(b0^2-a0*c0)^0.5)/a0
f<-function(t){pt(t,v,λ)-α}#非心ｔ分布
t<-uniroot(f,c(0.1*t,2*t))$root#収束
c("n"=n,"k"=round(t/n^0.5,4))
}
この関数を使って p  0.05,   0.75 として n  50,100,150,200 について K 値を算出すると以
下のようになります。
> kcalc(50,0.05,0.75)
n
k
50.0000 1.8109
> kcalc(100,0.05,0.75)
n
k
100.0000 1.7576
> kcalc(150,0.05,0.75)
n
k
150.0000 1.7354
> kcalc(200,0.05,0.75)
n
k
200.0000 1.7225
この K 値を使って、上記と同様なシミュレーションを行うと下記のようになります。
> #5%値の実験(その 2)
> test2<-function(n,n1,μ,σ,p,α){
+ kcalc<-function(n,p,α){#k 値の算出関数
+ if (p<=0.5){p<-1-p}
+ v<-(n-1)#自由度
+ up<-qnorm(p,0,1)#標準正規分布のｐ%点
+ λ<-up*n^0.5#非心度
+ a0<-(1-1/4/v)^2-up/2/v#初期値
+ b0<-λ*(1-1/4/v)#初期値
+ c0<-λ^2-up^2#初期値
+ t<-(b0+(b0^2-a0*c0)^0.5)/a0
+ f<-function(t){pt(t,v,λ)-α}#非心ｔ分布
+ t<-uniroot(f,c(0.1*t,2*t))$root#収束
+ atai<-round(t/n^0.5,4)
+ return(atai)
+}
+ K<-kcalc(n1,p,α)
+ Ln<-0
+ Un<-0
+ for (i in 1:n){
+ x<-rnorm(n1,μ,σ)#乱数発生
+ tl[i]<-mean(x)-K*sd(x)#5%値算出
+ if (tl[i]<qnorm(p,μ,σ)){Ln<-Ln+1}
+ if (tl[i]>qnorm(p,μ,σ)){Un<-Un+1}
+}
+ c(n,Ln,Un,Ln+Un)
+}
4-23
P％値の算出
> test2(1000,50,100,200,0.05,0.75)
[1] 1000 737 263 1000
> test2(1000,100,100,200,0.05,0.75)
[1] 1000 747 253 1000
> test2(1000,150,100,200,0.05,0.75)
[1] 1000 758 242 1000
> test2(1000,200,100,200,0.05,0.75)
[1] 1000 743 257 1000
真の 5％値を下回る個数が 750 前後になり下側信頼水準 75％に近いことが分かると思います。
表- 4 K 値を用いた 5％値算出実験結果（繰り返し数 1000 回）
真の 5％値の対して
発生回数
2.2
小さくなる個数
K値
50
737
1.8109
100
747
1.7576
150
758
1.7354
200
743
1.7225
非正規分布
2.2.1
対数正規分布
対数正規分布とはデータの対数値が正規分布に従うことですので、17）式で計算できることに
なります。
TL p  Exp(es  k p   es )
ここに
・・・ 17）
es , es ：対数正規分布の平均値・標準偏差の推定値
なお、等価正規分布パラメータを使って 18）式としても表現できます。
TL p   es  k p    es
ここに
  es 
  1  p 
f log  x p 
x p  Flog
f log 
1
・・・ 18）
 es  x p    es   1  p 
 p, es , es 
 ：対数正規分布の確率密度関数
4-24
P％値の算出
1
Flog

2.2.2
 ：対数正規分布のクォンタイル関数
ワイブル分布
パラメータ  200, m  12 の 2P ワイブル母集団から大きさ n1  50,100,150,200 個の乱数を
発生させ、最尤法を用いてパラメータ  es , mes を推定し、そのパラメータから 5％値を計算します。
そして、これを n  1000 回繰り返して真の 5％値（= 156.148)より大きい個数（ U n ）と小さい個
数（ Ln )を計算してみます。
（ n  U n  Ln ）
> #5%値の実験(その 3)
> test3<-function(n,n1,η,m,p,α){
+ Ln<-0
+ Un<-0
+ tl<-0
+ for (i in 1:n){
+ x<-rweibull(n1,m,η)#乱数発生
+ #2P ワイブル分布のパラメータ推定
+ LL2p<-function(para){#対数尤度関数の定義（2P ワイブル、完全）
+ atai<-sum(log(dweibull(x,para[1],para[2],log=FALSE)))
+ return(atai)
+}
+ fm<-function(mes){#m 関数定義
+ s1<-gamma(1+1/mes)
+ s2<-gamma(1+2/mes)
+ atai<-((n1-1)/n1)^0.5*sd(x)/mean(x)-(s2-s1^2)^0.5/s1
+ return(atai)
+}
+ mes0<-uniroot(fm,c(m/2,m*2))$root#m 初期値
+ ηes0<-mean(x)/gamma(1+1/mes0)#η初期値
+ para<-c(mes0,ηes0)#パラメータ初期値
+ para2p<-optim(para, LL2p, control=list(fnscale=-1))$par
+ ηes<-para2p[2]
+ mes<-para2p[1]
+ tl[i]<-qweibull(p,mes,ηes)#5%値算出
+ if (tl[i]<qweibull(p,m,η)){Ln<-Ln+1}
+ if (tl[i]>qweibull(p,m,η)){Un<-Un+1}
+}
+ c(n,Ln,Un,Ln+Un)
+}
> test3(1000,50,200,12,0.05,0.75)
[1] 1000 452 548 1000
> test3(1000,100,200,12,0.05,0.75)
[1] 1000 473 527 1000
> test3(1000,150,200,12,0.05,0.75)
[1] 1000 465 535 1000
> test3(1000,200,200,12,0.05,0.75)
[1] 1000 466 534 1000
結果を表－5 に示しますが、標本から求めた 5％値は真の 5％値より大きい・小さい確率がほぼ
50％となることが分かります。
4-25
P％値の算出
表- 5 2P ワイブルの 5％値算出実験結果（繰り返し数 1000 回）
発生回数
真の 5％値の対して
小さくなる個数
大きいなる個数
50
452
548
100
473
527
150
465
535
200
466
534
正規分布と同様に、何らかの下側信頼水準  % を設定手法について考えてみます。
ワイブル分布の p % 値を x p とすると x p における等価正規分布パラメータを使って x p は 19）式
で表現されます。
x p   m  u p    m
・・・ 19）
ここに、 u p ：標準正規分布の p % 点
同様にして、標本データから求められたパラメータ（ es , mes ,  es ）からの推定 p % 値（ x pes ）
は 20）式で表現できることになります。
x pes   m es  u p    m es
ここに
・・・ 20）
 m es ,   m es ： x pes における等価正規分布パラメータ
すなわち、下側信頼水準  % を設定すると正規分布と同様な 21)式を満たす k p が存在すること
になります。
Pr   m es  k p    m es   m  u p    m   
・・・ 21）
k p は正規分布の場合には非心 t 分布に従うことから算出できましたが、ワイブル分布の場合に
どのような分布に従うかは分っていません。
（筆者だけかも？）
4-26
P％値の算出
ここでは、ワイブル分布の下側信頼水準  ％の p % 値を正規分布で求めた k p と等価正規パラ
メータを使って算出する 22）式について考えてみたいと思います。
TL p   m es  k p    m es
ここに
  m es 
  1  p 
f wei  x p 
・・・ 22）
1
 m es  x p    m es   1  p  、 x p   es  es   ln 1  p  mes
パラメータ  200, m  12 の 2P ワイブル母集団から大きさ n1  50,100,150,200 個の乱数を
発生させ、最尤法を用いてパラメータ  es , mes を推定し 22）式から 5％値を計算します。そして、
これを n  1000 回繰り返して真の 5％値（= 156.148)よりと小さい個数（ Ln )を計算してみます。
> #5%値の実験(その 4)
> test4<-function(n,n1,η,m,p,K){
+ Ln<-0
+ Un<-0
+ tl<-0
+ for (i in 1:n){
+ x<-rweibull(n1,m,η)#乱数発生
+ #2P ワイブル分布のパラメータ推定
+ LL2p<-function(para){#対数尤度関数の定義（2P ワイブル、完全）
+ atai<-sum(log(dweibull(x,para[1],para[2],log=FALSE)))
+ return(atai)
+}
+ fm<-function(mes){#m 関数定義
+ s1<-gamma(1+1/mes)
+ s2<-gamma(1+2/mes)
+ atai<-((n1-1)/n1)^0.5*sd(x)/mean(x)-(s2-s1^2)^0.5/s1
+ return(atai)
+}
+ mes0<-uniroot(fm,c(m/2,m*2))$root#m 初期値
+ ηes0<-mean(x)/gamma(1+1/mes0)#η初期値
+ para<-c(mes0,ηes0)#パラメータ初期値
+ para2p<-optim(para, LL2p, control=list(fnscale=-1))$par
+ ηes<-para2p[2]
+ mes<-para2p[1]
+ xp<-qweibull(p,mes,ηes)#推定 5%値算出
+ σn<-dnorm(qnorm(pweibull(xp,mes,ηes),0,1),0,1)/dweibull(xp,mes,ηes)#等価正規分布パラメータσ
+ μn<-xp-σn*qnorm(pweibull(xp,mes,ηes),0,1)#等価正規パラメータμ
+ tl[i]<-μn-K*σn#5%値算出
+ if (tl[i]<qweibull(p,m,η)){Ln<-Ln+1}
+ if (tl[i]>qweibull(p,m,η)){Un<-Un+1}
+}
+ c(n,Ln,Un,Ln+Un)
+}
> test4(1000,50,200,12,0.05,1.8109)
[1] 1000 745 255 1000
> test4(1000,100,200,12,0.05,1.7576)
[1] 1000 768 232 1000
4-27
P％値の算出
> test4(1000,150,200,12,0.05,1.7354)
[1] 1000 758 242 1000
> test4(1000,200,200,12,0.05,1.7225)
[1] 1000 770 230 1000
計算結果を表-6 に示しますが、22）式から算出される 5％値は真の 5％値を下回る個数が 750
前後になり下側信頼水準 75％に近いことが分かると思います。
表- 6 K 値を用いたワイブル分布の 5％値算出実験結果（繰り返し数 1000 回）
発生回数
真の 5％値の対して
小さくなる個数
K値
50
745
1.8109
100
768
1.7576
150
758
1.7354
200
770
1.7225
３．計算例
データが得られた場合、①母集団パラメータを推定し、その後② TL 値を算出することになりま
す。
完全データ扱いのパラメータ算出し、 TL 値を算出する関数例を下記に示します。
#パラメータ及び TL 値の算出（完全データ扱い）
paraTLcalfull<-function(x,p,α){#最尤法によるパラメータ推定（完全データ扱い）
#x:データ
#p:p%
#α:下側信頼水準
n<-length(x)#データ大きさ
#ｋ値算出
if (p<=0.5){p<-1-p}
v<-(n-1)#自由度
up<-qnorm(p,0,1)#標準正規分布のｐ%点
λ<-up*n^0.5#非心度
a0<-(1-1/4/v)^2-up/2/v#初期値
b0<-λ*(1-1/4/v)#初期値
c0<-λ^2-up^2#初期値
t<-(b0+(b0^2-a0*c0)^0.5)/a0
f<-function(t){pt(t,v,λ)-α}#非心ｔ分布
t<-uniroot(f,c(0.1*t,2*t))$root#収束
K<-round(t/n^0.5,4)#K 値
#パラメータ算出
#正規分布
LL.norm<-function(para){#対数尤度関数(完全）
atai<-sum(dnorm(x,mean=para[1],sd=para[2],log=TRUE))
return(atai)
4-28
P％値の算出
}
μes<-mean(x)
σes<-((n-1)/n)^0.5*sd(x)
LLmax.norm<-LL.norm(c(μes,σes))#最大対数尤度
AIC.norm<-2*2-2*LLmax.norm
#対数正規分布
LL.lnorm<-function(para){#対数尤度関数(完全)
atai<-sum(dlnorm(x,meanlog=para[1],sdlog=para[2],log=TRUE))
return(atai)
}
lx<-log(x)#データ log 化
λes<-mean(lx)
ξes<-((n-1)/n)^0.5*sd(lx)
LLmax.lnorm<-LL.lnorm(c(λes,ξes))#最大対数尤度
AIC.lnorm<-2*2-2*LLmax.lnorm
#2P ワイブル分布
LL.weibull<-function(para){#対数尤度関数(完全)
atai<-sum(dweibull(x,para[1],para[2],log=TRUE))
return(atai)
}
fm<-function(mes){#fm 関数(2P ワイブル用）
s1<-gamma(1+1/mes)
s2<-gamma(1+2/mes)
atai<-((n-1)/n)^0.5*sd(x)/mean(x)-(s2-s1^2)^0.5/s1
return(atai)
}
mes0<-1
mes0<-uniroot(fm,c(1,15))$root#m 初期値
ηes0<-mean(x)/gamma(1+1/mes0)#η初期値
para<-c(mes0,ηes0)#パラメータ初期値
para<-optim(para, LL.weibull, method="SANN",control=list(fnscale=-1))$par#2P ワイブル分布のパラメー
タ推定"シミュレーテッドアニーリング法"
para<-optim(para, LL.weibull, method="BFGS",control=list(fnscale=-1))$par#2P ワイブル分布のパラメータ
推定"準ニュートン法"
mes<-para[1]
ηes<-para[2]
LLmax.weibull<-LL.weibull(para)#最大対数尤度
AIC.weibull<-2*2-2*LLmax.weibull
#最小 AIC
minAIC<-min(AIC.norm,AIC.lnorm,AIC.weibull)
if(minAIC==AIC.norm){fname<-"正規分布"}
if(minAIC==AIC.lnorm){fname<-"対数正規分布"}
if(minAIC==AIC.weibull){fname<-"2P ワイブル分布"}
#図化
hist(x,main="最尤法に よる推定（完全 データ扱い )¥n（正規：黒）（対数正 規：青）（ワイブル：赤） ",
freq=F,ylim=c(0,1.2*dnorm(μes,μes,σes)))
curve(dnorm(x,μes,σes),add=T,col="black",lwd=ifelse(minAIC==AIC.norm,2,1))
curve(dlnorm(x,λes,ξes),add=T,col="blue",lwd=ifelse(minAIC==AIC.lnorm,2,1))
curve(dweibull(x,mes,ηes),add=T,col="red",lwd=ifelse(minAIC==AIC.weibull,2,1))
#TL 値の算出
TLnorm<-μes-K*σes#正規分布 TL 値
TLlnorm<-exp(λes-K*ξes)#対数正規分布 TL 値
4-29
P％値の算出
#2p ワイブル
xp<-qweibull(p,mes,ηes)#推定 p%値算出
σn<-dnorm(qnorm(pweibull(xp,mes,ηes),0,1),0,1)/dweibull(xp,mes,ηes)#等価正規分布パラメータσn
μn<-xp-σn*qnorm(pweibull(xp,mes,ηes),0,1)#等価正規パラメータμn
TLweibull<-μn-K*σn#2P ワイブル TL 値
p0<-p
if (p>=0.5){p0<-1-p}
kekka<-c(
c("n"=n,"p"=p0,"α"=α,"K"=K),
c("μes"=μes,"σes"=σes,"AICnorm"=AIC.norm),
c("λes"=λes,"ξes"=ξes,"AIClnorm"=AIC.lnorm),
c("mes"=mes,"ηes"=ηes,"AICweibull"=AIC.weibull),
c("AIC 最小分布:"=fname),
c("TLnorm"=TLnorm,"TLlnorm"=TLlnorm,"TLweibull"=TLweibull)
)
return(kekka)
}
正規分布乱数（   100,   20 ）を 1000 個発生させ、下側信頼水準   75% の 5％値の TL を
求めた例を下記に示します。
> x<-rnorm(1000,100,20)
> kekka<-paraTLcalfull(x,0.05,0.75)
> kekka
n
p
α
"1000"
"0.05"
"0.75"
K
μes
σes
"1.6784" "100.944079580483" "20.2931478713422"
AICnorm
λes
ξes
"8862.44363797644" "4.59292149201174" "0.213445695477769"
AIClnorm
mes
ηes
"8938.97438563091" "5.41592797533742" "109.17202743562"
AICweibull
AIC 最小分布:
TLnorm
"8895.33456294908"
"正規分布" "66.8840601932219"
TLlnorm
TLweibull
"69.0391846387337" "77.2049942552876"
0.010
0.000
Density
0.020
最尤法による推定（完全データ扱い)
（正規：黒）（対数正規：青）（ワイブル：赤）
50
100
x
4-30
150