信頼性モデルの統計解析　共立出版　真壁肇・宮村鐵夫・鈴木和幸先生

#信頼性モデルの統計解析 

#p.42 指数分布に関する推定と検定

n=20;r=5

t=sort(sample(10:100,r,replace=T))

#MTBF

mu_hat=(sum(t)+(n-r)*t[length(t)])/r

#区間推定

alpha=0.05

mu_und=2*sum(t)/qchisq(1-alpha/2,df=2*r)

mu_up=2*sum(t)/qchisq(alpha/2,df=2*r)

#故障率の信頼区間

lam_up=1/mu_und

lam_und=1/mu_up

#時刻tにおける信頼度R(t)

R_L=exp(-t/mu_und);R_U=exp(-t/mu_up)

plot(t,R_L,xlim=c(min(t),max(t)),ylim=c(min(c(R_L,R_U)),max(c(R_L,R_U))),xlab="t",ylab="R",col=2,type="o")

par(new=T)

plot(t,R_U,xlim=c(min(t),max(t)),ylim=c(min(c(R_L,R_U)),max(c(R_L,R_U))),xlab="t",ylab="R",col=3,type="o")



#p.49 観測打ち切りの場合：

n=20;

t=sort(sample(10:100,n,replace=T))

X=sort(sample(10:100,n,replace=T))

delta=ifelse(t<=X,1,0)

#MTBFの最尤推定量

mu_hat=(sum((1-delta)*X)+sum(delta*t))/r

alpha=0.05;r=sum(delta)

#MTBFの信頼区間

p_L=r*qf(alpha/2,2*r,2*(n-r+1))/((n-r+1)+r*qf(alpha/2,2*r,2*(n-r+1)))

p_U=(r+1)*qf(1-alpha/2,2*(r+1),2*(n-r))/(n-r+(r+1)*qf(1-alpha/2,2*(r+1),2*(n-r)))

mu_L=-t[1]/log(1-p_U);mu_U=-t[1]/log(1-p_L)



#p.90 B.位置母数ｒが未知の場合の母数r,mu,sigmaの推定

n=10

t=sample(10:100,n,replace=T)

#Broyden Quasi-Newton Algorithm p.213

f=function(x){

f1=sum(log(t-x[1])-x[2])/(x[3]^2)

f2=n/x[3]+sum((log(t-x[1])-x[2])^2)/(x[3]^2)

f3=sum((log(t-x[1])-x[2])/(t-x[1]))/(x[3]^2)+sum(1/(t-x[1]))

return(c(f1,f2,f3))

}

X=rep(10,3)

ite=100000;eta=0.001

H=diag(f(X))

pre_f_val=10;f_val=pre_f_val-1

for(l in 1:ite){

if(abs(pre_f_val-f_val)>eta*10^(-3)){  

pre_f_val=f_val

X_pre=X  

X=X-eta*H%*%f(X)  

X[1]=ifelse(X[1]>min(t),min(t)-eta*10^(-3),X[1])

s=X-X_pre

y=f(X)-f(X_pre)

H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H

f_val=sum(abs(f(X)))

print(sum(abs(f(X))))

}

}  

#p.155

library(dplyr)

z=c(20.1,8.5,4.2,2,13.2,8.2,6.1,13.5,8.4,6.8,8.2,11.3,3.9,6.2,13.5,18.2,10.1,3.8,12.6,5.2)

N=50

r=length(z)

#Broyden Quasi-Newton Algorithm p.213

f=function(x){

f1=-(N-r)/(x[1]+x[2])+r/x[1]-sum(z)

f2=(N-r)*(1/x[2]-1/sum(x))-sum(z)

return(c(f1,f2))

}

X=rep(1,2)

ite=10^6

H=diag(f(X))

eta=10^(-5)

for(l in 1:ite){

X_pre=X  

X=X-eta*H%*%f(X)  

s=X-X_pre

y=f(X)-f(X_pre)

H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H

loglik=(N-r)*log(X[2]/sum(X))+r*log(X[1])-sum(X)*sum(z)

print(loglik)

}


#信頼性モデルの統計解析 

#ベイズ法　p.172

library(dplyr)

n=20;r=5

t=sample(c(1:100),n,replace=T)

lam=n/sum(t)

g_lam=1/lam

beta=0.1;alpha=0.2

#E[lam|t1,t2,...,tn]

(beta+n)/(alpha+sum(t))

#V[lam|t1,t2,...,tn]

(beta+n)/((alpha+sum(t))^2)

rgamma(100,beta+n,(alpha+sum(t)))

#信頼性のベイズ推定量

#二項サンプリング

r0=1;n0=2

n=6;r=5

#E[P|r]

(r+r0)/(n+n0)

#V[P|r]

(r+r0)*(n+n0-r-r0)/((n+n0+1)*(n+n0)^2)

alpha=0.05

und=(r+r0)/(r+r0+(n+n0-r-r0)*qf(1-alpha/2,2*n+2*n0-2*r-2*r0,2*r+2*r0))

up=(r+r0)*qf(1-alpha/2,2*r+2*r0,2*n+2*n0-2*r-2*r0)/(n+n0-r-r0+(r+r0)*qf(1-alpha/2,2*r+2*r0,2*n+2*n0-2*r-2*r0))

#パスカルサンプリング

n=10;s=3

r0=1;n0=2

#E[P|n]

p_hat=(n+r0-s)/(n+n0)

#V[P|n]

(n+r0-2)*(n0-r0+s)/((n+n0+1)*(n+n0)^2)


alpha=0.05

und=(n+r0-s)/((n+r0-s)+(n0-r0+s)*qf(1-alpha/2,2*(n0-r0+s),2*(n+r0-2)))

up=(n+r0-s)*qf(1-alpha/2,2*(n+r0-s),2*(n0-r0+s))/(n0-r0+s+(n+r0-s)*qf(1-alpha/2,2*(n+r0-s),2*(n0-r0+s)))



#B.計量値の場合：

#a.故障率のベイズ推定

r=20;n=30

t=sample(10:20,r,replace=T)

total_time=sum(t)+t[length(t)]*(n-r)

beta=2;alpha=1

#E[lam|total_time](ベイズ推定量)

lam_hat=(beta+r)/(alpha+total_time)


a=0.05

und=qchisq(a/2,df=2*(beta+r))/(2*(alpha+total_time))

up=qchisq(1-a/2,df=2*(beta+r))/(2*(alpha+total_time))

#MTBF(mean time between failure)

#E[MTBF|total_time]

(alpha+total_time)/(beta+r-1)



#8.3　ベイズ法による信頼性実証試験

t0=10

N=1000

lam0=0.1;lam1=0.3

alpha=10;beta=0.5

r_star=8

vec1=c();vec2=c()

for(j in (r_star+1):N){

fun=function(z){

value=exp(-z)*z^(beta+j-1)

return(value)

}

vec1=c(vec1,prod(rep(t0,j)/(c(1:j)*(alpha+t0)))*((1/(alpha+t0))^(beta-1))*integrate(fun,0,lam0*(alpha+t0))$value)

vec2=c(vec2,prod(rep(t0,j)/(c(1:j)*(alpha+t0)))*((1/(alpha+t0))^(beta-1))*gamma(beta+j))

}

vec_data=data.frame(vec1=vec1,vec2=vec2)%>%filter(vec1>0)

gamma_star=sum(vec_data$vec1)/sum(vec_data$vec2)