ワイブル分布の確率密度関数への近似

R
Python
Author

Nobukuni Hyakutake

Published

2023-02-20

1 目的

寿命データを最尤推定法でワイブル分布の確率密度関数に近似します。

2 解析するデータの読み込み

寿命のデータ

library(tidyverse)
library(plotly)
library(gt)
dsu01<-read_csv("
  life
  26
  30
  32
  30
  25
  27
  24
  23
  29
  29
  29
  31
  38
  33
  32
  29
  30
  18
  28
  25
  24
  23
  28
  29
  27
  26
  28
  28
  28
  27
  ",show_col_types = FALSE
)
dsu11<-count(dsu01,life) |> 
  mutate(n=n/sum(n))
f1<-plot_ly(x=dsu11$life,y=dsu11$n,type='bar',name="Actual")
f1

3 ワイブル分布の確率密度関数

下記式のηとmを最尤推定します。

\[ f(t)=\frac{m}{\eta}\left(\frac{t}{\eta}\right)^{m-1}e^{-\left(\frac{t}{\eta}\right)^m} \]

4 Rでのワイブル近似

4.1 パッケージの読み込み

library(survival)

4.2 ワイブルのパラメーターを算出

survregを使う。

t<-survreg(Surv(dsu01$life)~1, dist="weibull")
t
Call:
survreg(formula = Surv(dsu01$life) ~ 1, dist = "weibull")

Coefficients:
(Intercept) 
   3.382502 

Scale= 0.1275018 

Loglik(model)= -82.8   Loglik(intercept only)= -82.8
n= 30

survregの結果は、上述の確率密度関数のパラメーターとして使うには換算する必要がある。

\(m\)

m_fit<-1/t$scale
m_fit
[1] 7.843027

\(\eta\)

eta_fit<-exp(unname(t$coefficients))
eta_fit
[1] 29.44436

4.3 グラフ化

tsr<-c(1:max(dsu01$life))
ttr<-str_c("Scale parameter η: ",as.character(round(eta_fit,3)),", Shape parameter m: ",as.character(round(m_fit,3)))
f2<-f1|> 
  add_trace(x=tsr,y=dweibull(tsr,shape=m_fit,scale=eta_fit),type="scatter",mode="line",name="Weibull fit") |>
  layout(title="Weibull fit with 'survival' package - R") |> 
  add_annotations(
    x= 0.03,
    y= 0.95,
    xref = "paper",
    yref = "paper",
    text = ttr,
    showarrow=FALSE)
f2

4.4 複数タイプから一気にηとmを計算

t3<-read_csv("
type,life
a,5
b,4
a,5
a,6
b,3
b,2
a,4
b,2
",show_col_types = FALSE)

t4<-split(t3,t3$type)

mm2<-function(i){
  t<-survreg(Surv(i$life)~1, dist="weibull")
  t1<-tibble(eta=exp(unname(t$coefficients)),m=1/t$scale)
  return(t1)
}

qq3<-map(t4,mm2) |> 
  bind_rows(.id="id") |> 
  rename(type=id) |> 
  arrange(type)
gt(qq3)
type eta m
a 5.308623 7.975578
b 3.058647 3.613229

4.5 R version

R.version.string
[1] "R version 4.2.2 (2022-10-31)"

5 Pythonでのワイブル近似

import numpy as np
import pandas as pd
from reliability.Fitters import Fit_Weibull_2P
import plotly.graph_objects as go

5.1 関数定義

pythonはワイブル関数の自作が必要

def weibull_d(m,t,eta):
  k=(m/eta)*(t/eta)**(m-1)*np.exp(-1*((t/eta)**m))
  return k

5.2 フィット

結果詳細。Alphaは尺度パラメーターη、Betaは形状パラメーターmに相当する。

dsu02=r.dsu01
life_max=int(max(dsu02['life']))
dsu03=dsu02['life'].values.tolist()
wb2 = Fit_Weibull_2P(failures=dsu03)
Results from Fit_Weibull_2P (95% CI):
Analysis method: Maximum Likelihood Estimation (MLE)
Optimizer: TNC
Failures / Right censored: 30/0 (0% right censored) 

Parameter  Point Estimate  Standard Error  Lower CI  Upper CI
    Alpha         29.4443        0.725964   28.0553   30.9021
     Beta         7.84309          1.0011   6.10715   10.0725 

Goodness of fit    Value
 Log-likelihood -82.7942
           AICc  170.033
            BIC  172.391
             AD   1.2344

5.3 グラフ化

tt="Scale parameter η: "+str(round(wb2.alpha,3))+", Shape parameter m: "+str(round(wb2.beta,3))
ts=np.arange(1,life_max+1)
ta=r.dsu11
fig0201=go.Bar(
        x=ta['life'], y=ta["n"],name='Actual',
        marker={"color": "blue"}
        )
fig0202=go.Scatter(
        x=ts, y=weibull_d(wb2.beta,ts,wb2.alpha), name='Weibull fit',
        line={"color": "orange"},
    )
layout=go.Layout(
        title={"text":"Weibull fit with 'leliability' package - Python "})
fig04=go.Figure(data=[fig0201,fig0202],layout=layout)
fig04=fig04.update_layout(annotations=[
  go.layout.Annotation(
    xref="paper",
    yref="paper",
    x=0.03,
    y=0.95,
    showarrow=False,
    text=tt
)
])
fig04

5.4 複数タイプから一気にηとmを計算

pt1=pd.DataFrame([
["a",5],
["b",4],
["a",5],
["a",6],
["b",3],
["b",2],
["a",4],
["b",2]
])
pt1.columns=["type","life"]
pt2=pt1.drop_duplicates("type")

5.4.1 タイプごとの結果

def mp(i):
  print("type: "+i)
  dpu01=pt1.loc[pt1["type"]==i]
  dpu02=dpu01['life'].values.tolist()
  wb21 = Fit_Weibull_2P(failures=dpu02)
  dpu03=pd.DataFrame({"type":[i],"eta":[wb21.alpha],"m":[wb21.beta]})
  return dpu03

dpu04=pd.concat((mp(r) for r in pt2["type"]),ignore_index=True)
type: a
Results from Fit_Weibull_2P (95% CI):
Analysis method: Maximum Likelihood Estimation (MLE)
Optimizer: TNC
Failures / Right censored: 4/0 (0% right censored) 

Parameter  Point Estimate  Standard Error  Lower CI  Upper CI
    Alpha         5.30868        0.351929   4.66185   6.04527
     Beta         7.97623          3.1118   3.71291   17.1349 

Goodness of fit    Value
 Log-likelihood -4.32774
           AICc  24.6555
            BIC  11.4281
             AD  3.03487 

type: b
Results from Fit_Weibull_2P (95% CI):
Analysis method: Maximum Likelihood Estimation (MLE)
Optimizer: TNC
Failures / Right censored: 4/0 (0% right censored) 

Parameter  Point Estimate  Standard Error  Lower CI  Upper CI
    Alpha         3.05865        0.449056   2.29382   4.07848
     Beta         3.61323         1.39996   1.69081   7.72139 

Goodness of fit    Value
 Log-likelihood -4.90322
           AICc  25.8064
            BIC   12.579
             AD  3.02246

5.4.2 まとめ結果

dpu04
  type     eta       m
0    a 5.30868 7.97623
1    b 3.05865 3.61323

5.5 参考サイト

Python package reliability - Fitting a specific distribution to data