修理する７ 「修理あり故障累積回数確率」の実例

• イベント・逆イベントの起きやすさ・閾値回数をパラメタにして「修理あり故障累積回数確率」がどうなるかをみてみる

# 修復ありポアッソン


matrixA<-function(p,q,n,open=FALSE){
	k<-p
	h<-q
	N<-n
	M<-matrix(0,N,N)

	M[1,1]<--k
	M[1,2]<-h
	M[2,1]<-k
	M[2,2]<--(k+h)
	M[2,3]<-h

	for(i in 2:(N-1)){
		M[i,i-1]<-k
		M[i,i]<--(k+h)
		M[i,i+1]<-h
	}

	M[N,N-1]<-k
	if(open){
		M[N,N]<--(k+h)
	}else{
		M[N,N]<--h
	}
	
	M
}


my.ShokiChi<-function(A,v0,t){
	n<-length(v0)
	# 固有値分解
	eigen.out<-eigen(A)
	V<-eigen.out[[2]]
	U<-solve(V)
	B<-diag(eigen.out[[1]])
	# 計算する
	x<-t
	Ys<-matrix(0,length(x),n)
	# 初期値も適当に
	Ys0<-v0
	for(i in 1:(length(x))){
		Bex<-diag(exp(eigen.out[[1]]*x[i]))
		#Aex<-V%*%Bex%*%U
		#Ys[i,]<-Aex%*%Ys0
		Ys[i,]<-V%*%(Bex%*%(U%*%Ys0))
		#print(Ys[i,])
	}
	ReYs<-Re(Ys)
	ReYs[which(ReYs<0)]<-0
	#print(range(ReYs))
	ReYs
}


# pは単位時間あたりの平均イベント生起回数
# qは単位時間当たりの平均修復イベント生起回数

# pとqとを与えて、「得点」が閾値以下である確率を求める

Syufuku.Poisson<-function(p,q,n,maxk=NULL,t=1,open=FALSE){
	k<-0:n
	if(!is.null(maxk)){
		k<-0:maxk
	}
	ret<-rep(0,length(k))
	if(p==0){
		if(q==0){
			ret<-rep(-1,length(k))
			
		}else{
			ret[which(k>=0)]<-1
		}
	}else{
		if(q==0){
			ret<-ppois(k,p*t)
			#print((ret))
		}else{
			A<-matrixA(p=p,q=q,n=n,open=open)
			v0<-c(1,rep(0,n-1))
			my.out<-my.ShokiChi(A,v0,t)
			ret<-cumsum(my.out)
			#print(range(ret))
		}
	}
	ret
}

################

t<-1
#kq<-2^(seq(from=-1,to=1,length=20))
kq<-seq(from=1,to=20,by=0.1)
kp<-kq
q<-kq
p<-kp
n<-20
maxk<-10

pqt<-expand.grid(p,q,t)
#out<-rep(0,length(pqt[,1]))
out<-matrix(0,length(pqt[,1]),maxk)
for(i in 1:length(pqt[,1])){
	#A<-matrixA(pqt[i,1],pqt[i,2],n)
	#print(A)
	out[i,]<-Syufuku.Poisson(p=pqt[i,1],q=pqt[i,2],t=pqt[i,3],n=n)[1:maxk]
}

#plot(pqt)


library(rgl)
plot3d(log(pqt[,1]),log(pqt[,2]),out[,1],zlim=c(0,1))
plot3d(pqt[,1],pqt[,2],out[,1],zlim=c(0,1))
points3d(pqt[,1],pqt[,2],out[,2],col=2)

plot3d(log(pqt[,1]),log(pqt[,2]),out[,1],type="l",zlim=c(0,1))

plot3d(pqt[,1],pqt[,2],out[,1],zlim=c(0,1))

for(i in 2:maxk){
	points3d(pqt[,1],pqt[,2],out[,i],col=i)

}

pp<-4
qq<-4
t<-1
n<-100
plot(Syufuku.Poisson(p=pp,q=qq,t=t,n=n),ylim=c(0,1))