システム信頼性解析法  日科技連　阿部俊一先生

```r

#データからの最良線形不偏推定量の求め方
#p.20

n=100

x=c(rnorm(n,2,10),rnorm(n,5,9))
sortx=sort(x)

#定数の打ち切り
x=sortx[-c(c(1:25),c(175:200))]

#F(x)を正規分布とした場合
mu=mean(x)
sd=sd(x)

z=(x-mu)/sd
rij=diag(sd^2,length(x))
ri=rep(mu,length(x))

#位置パラメータの計算
#Broyden Quasi-Newton Algorithm p.213

f=function(XX){
lam=XX[1]
XX=XX[-c(1)]
f1=sum(XX*ri)-mu
f2=t(rij)%*%XX-lam*ri
return(c(f1,f2))
}

X=sample(1:(2*n),length(x)+1)/(n*100)

H=diag(f(X))
eta=10^(-2)

while(sum(abs(f(X)))>10^(-4))({
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
print(sum(abs(f(X))))
})

lambda_a=X[1]
a=X[-c(1)]

#最良線形不変推定量(正規分布の場合は標本平均となるため)
mu_childer=sum(a*x)


#尺度パラメータの計算
#Broyden Quasi-Newton Algorithm p.213

f=function(XX){
lam=XX[1]
XX=XX[-c(1)]
f1=sum(XX*ri)-sd
f2=t(rij)%*%XX-lam*ri
return(c(f1,f2))
}

X=sample(1:(2*n),length(x)+1)/(n*100)

H=diag(f(X))
eta=10^(-2)

while(sum(abs(f(X)))>10^(-4))({
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
print(sum(abs(f(X))))
})


lambda_b=X[1]
b=X[-c(1)]

#最良線形不変推定量(正規分布の場合は標本平均となるため)
sig_childer=sum(b*x)

#生成された標本
x_p=mu_childer+sig_childer*qnorm(seq(0,1,0.01))

#V(mu_childer)
V_mu_childer=sum(t(t(c(a)))%*%t(c(a))*rij)


#V(sig_childer)
V_sig_childer=sum(t(t(c(b)))%*%t(c(b))*rij)


#C(mu_childer,sig_childer)
C_mu_sig=sum(t(t(c(a)))%*%t(c(b))*rij)


```

```r

#データからの最良線形不偏推定量の求め方
#p.34 指数分布

n=150

x=c(rexp(n,5))
x=sort(x)
r=100
x=x[1:r]

z=c(r-c(1:r)+1)*c(x[1],c(x[2:r]-x[1:(r-1)]))


#位置パラメータの計算
#Broyden Quasi-Newton Algorithm p.213
#sitaはc_iのiに依存しない

f=function(XX){
lam=XX[1]
XX=XX[-c(1)]
f1=sum(XX)-1
f2=2*XX+lam
return(c(f1,f2))
}

X=sample(1:(2*r+1),length(z)+1)/(r*100)

H=diag(f(X))
eta=10^(-2)

while(sum(abs(f(X)))>10^(-6))({
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
cost=sum(X[-1]^2)
print(sum(abs(f(X))))
})

lambda_sita=X[1]
c=X[-c(1)]

#最良線形不変推定量
sita_hat=sum(c*x)


```

```r

#データからの最良不偏推定量の求め方
#p.34 指数分布
#平均2乗誤差の場合

n=100

x=c(rexp(n*2,5))
x=sort(x)
r=100
x=x[1:r]

z=c(n-c(1:r)+1)*c(x[1],c(x[2:r]-x[1:(r-1)]))

#dQ/dc_j=0を解く
mat=diag(2,length(z))
mat[mat==0]=1
vec=rep(1,length(z))
c=solve(mat)%*%vec

sita_hathat=sum(c*z)

```
```r

#p.45
#3.2.7 muとsitaの推定　指数分布の場合

n=10
r=7
x=c(55,57,86,90,98,116,174)


loglik=function(param){
SITA=param[1]
MU=param[2]
return(-r*log(SITA)-(sum(x-MU)+(n-r)*(x[r]-MU))/(SITA))
}

mu=x[1]
sita=(sum(x)+(n-r)*x[r]-n*x[1])/r

#最尤推定での対数尤度
loglik(c(sita,mu))

#最良不変推定量
sita_hat=sita
mu_hat=mu-sita/n

#信頼区間
mu_confidence_inf=mu-(r*sita/(n*(r-1)))*qf(1-0.1/2,2,2*r-2)
mu_confidence_sup=mu-(r*sita/(n*(r-1)))*qf(0.05/2,2,2*r-2)
sita_confidence_inf=2*r*sita/(qchisq(1-0.1/2,2*r-2))
sita_confidence_sup=2*r*sita/(qchisq(0.1/2,2*r-2))

```
```r
#p.187 和の形のモデル


s=c(1,1,1,3,4,3,2,3,5)
t=c(15,11,10,18,12,11,11,13,12)
delta=c(rep(1,8),0)
z=matrix(c(rep(1,9),c(1,5,8,2,3,1,2,0,0),c(2,6,7,2,3,1,2,2,1)),nrow=9)
lam0=1

beta=rep(1,3)
H=diag(f(X))
eta=10^(-3)
ite=10^4

for(l in 1:ite){
dbeta=apply(z*c((1/z%*%beta)*delta)-z*c(lam0*(t-s)),2,sum)  
ddbeta=-t(z*c(delta/((z%*%beta)^2)))%*%z  
beta=beta-eta*solve(ddbeta)%*%dbeta
loglik=sum(delta*log(lam0*z%*%beta))-sum(lam0*z%*%beta*c(t-s))
print(loglik)
}
      

```

```r
#p.187 交互作用のモデル


s=c(1,1,1,3,4,3,2,3,5)
t=c(15,11,10,18,12,11,11,13,12)
delta=c(rep(1,8),0)
z=matrix(c(rep(1,9),c(1,5,8,2,3,1,2,0,0),c(2,6,7,2,3,1,2,2,1)),nrow=9)
lam0=1


beta=rep(1,3)
r=array(1,dim=c(3,3))
H=diag(f(X))
eta=10^(-7)
ite=50000
b=0.7
m1=10^(-1)
m2=10^(-1)

for(l in 1:ite){
value=z%*%beta
for(i in 1:nrow(z)){
value[i]=value[i]+sum(t(t(c(z[i,])))%*%t(c(z[i,]))*r)  
}
dbeta=apply(z*c(delta/value)-lam0*z*c(t-s),2,sum) 
dr=t(z*c(delta/value))%*%z-t(z*c(t-s))%*%z*lam0

m1=b*m1-eta*dbeta
m2=b*m2-eta*dr
beta=beta-m1
r=r-m2
loglik=sum(delta*log(lam0*value))-sum(lam0*value*c(t-s))
print(loglik)
}
      

```