Computational Codes PDF Print E-mail
Written by ahmadi-j   
Sunday, 28 August 2011 18:19

The R  codes used to obtain numerical results in the paper:

"Bayesian Reconstruction of the  Missing  Failure Times  in Exponential Distribution" JSCS

**************************

Example 1: GSCS-2010-0395
***************************
n<-12
Y<-array(n)
Y<-c(0.6,2.4,5.6,6.6,6.6,9,9.3,13.0,14.3,18.1,24.4,33.8)
f11=function(r){
xx=0
i=0
while (i<=(r-1)) {xx=c(xx,Y[i]);i=i+1}
sum(xx)
}
f12=function(s){
xx=0
i=(s+1)
while (i<=n) {xx=c(xx,Y[i]);i=i+1}
sum(xx)
}
***********************************************
Example1:  Point reconstruction
=======================
g11=function(r,l,s,k,m,a,b) {(Y[s]*(b+f11(r)+f12(s)+(l-r-k)*Y[r]+(s-l+k+1)*Y[s])^(1-(n-s+r+a+2))-
Y[r]*(b+f11(r)+f12(s)+(l-r+m+1)*Y[r]+(s-l-m)*Y[s])^(1-(n-s+r+a+2)))}
g12=function(r,l,s,k,m,a,b) {((b+f11(r)+f12(s)+(l-r-k)*Y[r]+(s-l+k+1)*Y[s])^(2-(n-s+r+a+2))-
(b+f11(r)+f12(s)+(l-r+m+1)*Y[r]+(s-l-m)*Y[s])^(2-(n-s+r+a+2)))/(k+m+1)/(2-(n-s+r+a+2))}
g1=function(r,l,s,k,a,b){
xx=0
m=0
while (m<=(s-l-1)) {xx=c(xx,choose(s-l-1,m)*(-1)^(s-l-1-m)*
(g11(r,l,s,k,m,a,b)-g12(r,l,s,k,m,a,b))/(k+m+1)/(1-(n-s+r+a+2)));m=m+1}
sum(xx)
}
g=function(r,l,s,a,b){
xx=0
k=0
while (k<=(l-r-1)) {xx=c(xx,choose(l-r-1,k)*(-1)^k*
g1(r,l,s,k,a,b));k=k+1}
sum(xx)
}
h11=function(r,l,s,k,m,a,b) {((b+f11(r)+f12(s)+(l-r-k)*Y[r]+(s-l+k+1)*Y[s])^(1-(n-s+r+a+2))-
(b+f11(r)+f12(s)+(l-r+m+1)*Y[r]+(s-l-m)*Y[s])^(1-(n-s+r+a+2)))/(k+m+1)/(1-(n-s+r+a+2))}
h1=function(r,l,s,k,a,b){
xx=0
m=0
while (m<=(s-l-1)) {xx=c(xx,choose(s-l-1,m)*(-1)^(s-l-1-m)*
((b+f11(r)+f12(s)+(l-r-k)*Y[r]+(s-l+k+1)*Y[s])^(1-(n-s+r+a+2))-(b+f11(r)+f12(s)+(l-r+m+1)*Y[r]+(s-l-m)*Y[s])^(1-(n-s+r+a+2)))/(k+m+1)/(1-(n-s+r+a+2)));m=m+1}
sum(xx)
}
h=function(r,l,s,a,b){
xx=0
k=0
while (k<=(l-r-1)) {xx=c(xx,choose(l-r-1,k)*(-1)^k*h1(r,l,s,k,a,b));k=k+1}
sum(xx)
}
for (r in 1:(n-2)){
for (s in (r+2):n){
for (l in (r+1):(s-1)){     
yhat=function(a,b)   g(r,l,s,a,b)/h(r,l,s,a,b)
}
}
}
A=array (,c(3*(n-2),(n-1),n))
for (r in 1:(n-2)){
for (s in (r+2):n){
for (l in (r+1):(s-1)){              
A[3*(r-1)+1,l,s]=yhat(0.5,0.5)
A[3*(r-1)+2,l,s]=yhat(1,1)
A[3*(r-1)+3,l,s]=yhat(2,2)
}
}
}
round(A,4)
B=array (,c(3*(n-2),n,n-1))
for (r in 1:(n-2)){
for (s in (r+2):n){
for (l in (r+1):(s-1)){     

B[3*(r-1)+1,s,l]=A[3*(r-1)+1,l,s]
B[3*(r-1)+2,s,l]=A[3*(r-1)+2,l,s]
B[3*(r-1)+3,s,l]=A[3*(r-1)+3,l,s]
}
}
}
round(B,4)
***********************************************
Example1:  Interval reconstruction
=======================
Intg1=function(r,l,s,x,k,a,b){
xx=0
m=0
while (m<=(s-l-1)) {xx=c(xx,choose(s-l-1,m)*(-1)^(s-l-1-m)*
((b+f11(r)+f12(s)+(l-r-k)*Y[r]+(s-l+k+1)*Y[s])^(1-(n-s+r+a+2))-
(b+f11(r)+f12(s)+(l-r-k)*Y[r]+(k+m+1)*x+(s-l-m)*Y[s])^(1-(n-s+r+a+2)))/(k+m+1)/(1-(n-s+r+a+2)));m=m+1}
sum(xx)
}
Intg=function(r,l,s,x,a,b){
xx=0
k=0
while (k<=(l-r-1)) {xx=c(xx,choose(l-r-1,k)*(-1)^k*
Intg1(r,l,s,x,k,a,b));k=k+1}
sum(xx)
}
Inth1=function(r,l,s,k,a,b){
xx=0
m=0
while (m<=(s-l-1)) {xx=c(xx,choose(s-l-1,m)*(-1)^(s-l-1-m)*
((b+f11(r)+f12(s)+(l-r-k)*Y[r]+(s-l+k+1)*Y[s])^(1-(n-s+r+a+2))-
(b+f11(r)+f12(s)+(l-r+m+1)*Y[r]+(s-l-m)*Y[s])^(1-(n-s+r+a+2)))/(k+m+1)/(1-(n-s+r+a+2)));m=m+1}
sum(xx)
}
Inth=function(r,l,s,a,b){
xx=0
k=0
while (k<=(l-r-1)) {xx=c(xx,choose(l-r-1,k)*(-1)^k*
Inth1(r,l,s,k,a,b));k=k+1}
sum(xx)
}
for (r in 1:(n-2)){
for (s in (r+2):n){
for (l in (r+1):(s-1)){     
Intlow=function(a,b,c) {
Intf1=function(x) { Intg(r,l,s,x,a,b)/Inth(r,l,s,a,b)-(1-c/2)}
uniroot(Intf1, lower=Y[r] , upper=Y[s])$root
}
Intup=function(a,b,c) {
Intf2=function(x) { Intg(r,l,s,x,a,b)/Inth(r,l,s,a,b)-(c/2)}
uniroot(Intf2, lower=Y[r] , upper=Y[s])$root
}
}
}
}
c=0.05
AInt=array (,c(3*(n-2),2*(n-1),n))
for (r in 1:(n-2)){
for (s in (r+2):n){
for (l in (r+1):(s-1)){     

AInt[3*(r-1)+1,2*(l-1),s]=Intlow(0.5,0.5,c)
AInt[3*(r-1)+1,2*(l-1)+1,s]=Intup(0.5,0.5,c)

AInt[3*(r-1)+2,2*(l-1),s]=Intlow(1,1,c)
AInt[3*(r-1)+2,2*(l-1)+1,s]=Intup(1,1,c)

AInt[3*(r-1)+3,2*(l-1),s]=Intlow(2,2,c)
AInt[3*(r-1)+3,2*(l-1)+1,s]=Intup(2,2,c)

}
}
}
round(AInt,4)
BInt=array (,c(3*(n-2),2*(n-1),n-1))
for (r in 1:(n-2)){
for (s in (r+2):n){
for (l in (r+1):(s-1)){     

BInt[3*(r-1)+1,2*(s-1)-1,l]=AInt[3*(r-1)+1,2*(l-1),s]
BInt[3*(r-1)+2,2*(s-1)-1,l]=AInt[3*(r-1)+2,2*(l-1),s]
BInt[3*(r-1)+3,2*(s-1)-1,l]=AInt[3*(r-1)+3,2*(l-1),s]
BInt[3*(r-1)+1,2*(s-1),l]=AInt[3*(r-1)+1,2*(l-1)+1,s]
BInt[3*(r-1)+2,2*(s-1),l]=AInt[3*(r-1)+2,2*(l-1)+1,s]
BInt[3*(r-1)+3,2*(s-1),l]=AInt[3*(r-1)+3,2*(l-1)+1,s]

}
}
}
round(BInt,4)

CInt=array (,c(3*(n-2),n-1,n))
for (r in 1:(n-2)){
for (s in (r+2):n){
for (l in (r+1):(s-1)){     

CInt[3*(r-1)+1,l,s]=BInt[3*(r-1)+1,2*(s-1),l]-BInt[3*(r-1)+1,2*(s-1)-1,l]
CInt[3*(r-1)+2,l,s]=BInt[3*(r-1)+2,2*(s-1),l]-BInt[3*(r-1)+2,2*(s-1)-1,l]
CInt[3*(r-1)+3,l,s]=BInt[3*(r-1)+3,2*(s-1),l]-BInt[3*(r-1)+3,2*(s-1)-1,l]

}
}
}
round(CInt,4)

******************************

Example 2: GSCS-2010-0395

******************************

n<-19
y<-array(n)
y<-c(162,20,271,320,393,508,539,629,709,777,884,1008,1101,1182,1463,1603,1984,2355,2880)
f11=function(r){
xx=0
i=0
while (i<=(r-1)) {xx=c(xx,y[i]);i=i+1}
sum(xx)
}
f12=function(s){
xx=0
i=(s+1)
while (i<=n) {xx=c(xx,y[i]);i=i+1}
sum(xx)
}
***********************************************
Exapme 2: Point reconstruction
=======================
g11=function(r,l,s,k,m,b,g,h) {(y[s]*(h-n*b+f11(r)+f12(s)+(l-r-k)*y[r]+(s-l+k+1)*y[s])^(1-(n-s+r+g+2))-
y[r]*(h-n*b+f11(r)+f12(s)+(l-r+m+1)*y[r]+(s-l-m)*y[s])^(1-(n-s+r+g+2)))}
g12=function(r,l,s,k,m,b,g,h) {((h-n*b+f11(r)+f12(s)+(l-r-k)*y[r]+(s-l+k+1)*y[s])^(2-(n-s+r+g+2))-
(h-n*b+f11(r)+f12(s)+(l-r+m+1)*y[r]+(s-l-m)*y[s])^(2-(n-s+r+g+2)))/(k+m+1)/(2-(n-s+r+g+2))}
g1=function(r,l,s,k,b,g,h){
xx=0
m=0
while (m<=(s-l-1)) {xx=c(xx,choose(s-l-1,m)*(-1)^(s-l-1-m)*
(g11(r,l,s,k,m,b,g,h)-g12(r,l,s,k,m,b,g,h))/(k+m+1)/(1-(n-s+r+g+2)));m=m+1}
sum(xx)
}
gg=function(r,l,s,b,g,h){
xx=0
k=0
while (k<=(l-r-1)) {xx=c(xx,choose(l-r-1,k)*(-1)^k*g1(r,l,s,k,b,g,h));k=k+1}
sum(xx)
}
h1=function(r,l,s,k,b,g,h){
xx=0
m=0
while (m<=(s-l-1)) {xx=c(xx,choose(s-l-1,m)*(-1)^(s-l-1-m)*
((h-n*b+f11(r)+f12(s)+(l-r-k)*y[r]+(s-l+k+1)*y[s])^(1-(n-s+r+g+2))-
(h-n*b+f11(r)+f12(s)+(l-r+m+1)*y[r]+(s-l-m)*y[s])^(1-(n-s+r+g+2)))/(k+m+1)/(1-(n-s+r+g+2)));m=m+1}
sum(xx)
}
hh=function(r,l,s,b,g,h){
xx=0
k=0
while (k<=(l-r-1)) {xx=c(xx,choose(l-r-1,k)*(-1)^k*h1(r,l,s,k,b,g,h));k=k+1}
sum(xx)
}
for (r in 1:(n-2)){
for (s in (r+2):n){
for (l in (r+1):(s-1)){     

yhat=function(b,g,h)   gg(r,l,s,b,g,h)/hh(r,l,s,b,g,h)

}
}
}
A=array (,c(3*(n-2),(n-1),n))
for (r in 1:(n-2)){
for (s in (r+2):n){
for (l in (r+1):(s-1)){     

A[3*(r-1)+1,l,s]=yhat(0.10263,0.5,1)
A[3*(r-1)+2,l,s]=yhat(0.2873,0.625,1.25)
A[3*(r-1)+3,l,s]=yhat(0.16697,0.83333,1.6667)

}
}
}
round(A,4)
***********************************************
Example 2: Interval reconstruction
=======================
Intg1=function(r,l,s,x,k,b,g,h){
xx=0
m=0
while (m<=(s-l-1)) {xx=c(xx,choose(s-l-1,m)*(-1)^(s-l-1-m)*
((h-n*b+f11(r)+f12(s)+(l-r-k)*y[r]+(s-l+k+1)*y[s])^(1-(n-s+r+g+2))-
(h-n*b+f11(r)+f12(s)+(l-r-k)*y[r]+(k+m+1)*x+(s-l-m)*y[s])^(1-(n-s+r+g+2)))/(k+m+1)/(1-(n-s+r+g+2)));m=m+1}
sum(xx)
}
Intg=function(r,l,s,x,b,g,h){
xx=0
k=0
while (k<=(l-r-1)) {xx=c(xx,choose(l-r-1,k)*(-1)^k*Intg1(r,l,s,x,k,b,g,h));k=k+1}
sum(xx)
}
Inth1=function(r,l,s,k,b,g,h){
xx=0
m=0
while (m<=(s-l-1)) {xx=c(xx,choose(s-l-1,m)*(-1)^(s-l-1-m)*
((h-n*b+f11(r)+f12(s)+(l-r-k)*y[r]+(s-l+k+1)*y[s])^(1-(n-s+r+g+2))-
(h-n*b+f11(r)+f12(s)+(l-r+m+1)*y[r]+(s-l-m)*y[s])^(1-(n-s+r+g+2)))/(k+m+1)/(1-(n-s+r+g+2)));m=m+1}
sum(xx)
}
Inth=function(r,l,s,b,g,h){
xx=0
k=0
while (k<=(l-r-1)) {xx=c(xx,choose(l-r-1,k)*(-1)^k*Inth1(r,l,s,k,b,g,h));k=k+1}
sum(xx)
}
for (r in 1:(n-2)){
for (s in (r+2):n){
for (l in (r+1):(s-1)){     


Intlow=function(b,g,h,c) {

Intf1=function(x) { Intg(r,l,s,x,b,g,h)/Inth(r,l,s,b,g,h)-(1-c/2)}

uniroot(Intf1, lower=y[r] , upper=y[s])$root
}
Intup=function(b,g,h,c) {
Intf2=function(x) { Intg(r,l,s,x,b,g,h)/Inth(r,l,s,b,g,h)-(c/2)}
uniroot(Intf2, lower=y[r] , upper=y[s])$root
}
}
}
}
c=0.05
Intlow(0.10263,0.5,1,c)
AInt=array (,c(3*(n-2),2*(n-1),n))
for (r in 1:(n-2)){
for (s in (r+2):n){
for (l in (r+1):(s-1)){     

AInt[3*(r-1)+1,2*(l-1),s]=Intlow(0.10263,0.5,1,c)
AInt[3*(r-1)+1,2*(l-1)+1,s]=Intup(0.10263,0.5,1,c)
AInt[3*(r-1)+2,2*(l-1),s]=Intlow(0.2873,0.625,1.25,c)
AInt[3*(r-1)+2,2*(l-1)+1,s]=Intup(0.2873,0.625,1.25,c)
AInt[3*(r-1)+3,2*(l-1),s]=Intlow(0.16697,0.83333,1.6667,c)
AInt[3*(r-1)+3,2*(l-1)+1,s]=Intup(0.16697,0.83333,1.6667,c)

}
}
}
round(AInt,4)
BInt=array (,c(3*(n-2),2*(n-1),n-1))
for (r in 1:(n-2)){
for (s in (r+2):n){
for (l in (r+1):(s-1)){     

BInt[3*(r-1)+1,2*(s-1)-1,l]=AInt[3*(r-1)+1,2*(l-1),s]
Int[3*(r-1)+2,2*(s-1)-1,l]=AInt[3*(r-1)+2,2*(l-1),s]
BInt[3*(r-1)+3,2*(s-1)-1,l]=AInt[3*(r-1)+3,2*(l-1),s]
BInt[3*(r-1)+1,2*(s-1),l]=AInt[3*(r-1)+1,2*(l-1)+1,s]
BInt[3*(r-1)+2,2*(s-1),l]=AInt[3*(r-1)+2,2*(l-1)+1,s]
BInt[3*(r-1)+3,2*(s-1),l]=AInt[3*(r-1)+3,2*(l-1)+1,s]

}
}
}
round(BInt,4)
CInt=array (,c(3*(n-2),n-1,n))
for (r in 1:(n-2)){
for (s in (r+2):n){
for (l in (r+1):(s-1)){     

CInt[3*(r-1)+1,l,s]=BInt[3*(r-1)+1,2*(s-1),l]-BInt[3*(r-1)+1,2*(s-1)-1,l]
CInt[3*(r-1)+2,l,s]=BInt[3*(r-1)+2,2*(s-1),l]-BInt[3*(r-1)+2,2*(s-1)-1,l]
CInt[3*(r-1)+3,l,s]=BInt[3*(r-1)+3,2*(s-1),l]-BInt[3*(r-1)+3,2*(s-1)-1,l]

}
}
}
round(CInt,4)

*************************

Simulation results: GSCS-2010-0395

****************************

w=10000
n=20
a=0.05
g=0.5
h=1
b=0.10263
c=1.05410
ss=rgamma(1,g,1/h)
mm=rexp(1,c*ss)+b
ss
mm
y<-array(n)
f11=function(r){
xx=0
i=0
while (i<=(r-1)) {xx=c(xx,y[i]);i=i+1}
sum(xx)
}
f12=function(s){
xx=0
i=(s+1)
while (i<=n) {xx=c(xx,y[i]);i=i+1}
sum(xx)
}
g11=function(r,l,s,k,m) {(y[s]*(h-n*b+f11(r)+f12(s)+(l-r-k)*y[r]+(s-l+k+1)*y[s])^(1-(n-s+r+g+2))-
y[r]*(h-n*b+f11(r)+f12(s)+(l-r+m+1)*y[r]+(s-l-m)*y[s])^(1-(n-s+r+g+2)))}

g12=function(r,l,s,k,m) {((h-n*b+f11(r)+f12(s)+(l-r-k)*y[r]+(s-l+k+1)*y[s])^(2-(n-s+r+g+2))-
(h-n*b+f11(r)+f12(s)+(l-r+m+1)*y[r]+(s-l-m)*y[s])^(2-(n-s+r+g+2)))/(k+m+1)/(2-(n-s+r+g+2))}
g1=function(r,l,s,k){
xx=0
m=0
while (m<=(s-l-1)) {xx=c(xx,choose(s-l-1,m)*(-1)^(s-l-1-m)*
(g11(r,l,s,k,m)-g12(r,l,s,k,m))/(k+m+1)/(1-(n-s+r+g+2)));m=m+1}
sum(xx)
}
gg=function(r,l,s){
xx=0
k=0
while (k<=(l-r-1)) {xx=c(xx,choose(l-r-1,k)*(-1)^k*g1(r,l,s,k));k=k+1}
sum(xx)
}
h1=function(r,l,s,k){
xx=0
m=0
while (m<=(s-l-1)) {xx=c(xx,choose(s-l-1,m)*(-1)^(s-l-1-m)*
(h-n*b+f11(r)+f12(s)+(l-r-k)*y[r]+(s-l+k+1)*y[s])^(1-(n-s+r+g+2))-
(h-n*b+f11(r)+f12(s)+(l-r+m+1)*y[r]+(s-l-m)*y[s])^(1-(n-s+r+g+2)))/(k+m+1)/(1-(n-s+r+g+2)));m=m+1}
sum(xx)
}
hh=function(r,l,s){
xx=0
k=0
while (k<=(l-r-1)) {xx=c(xx,choose(l-r-1,k)*(-1)^k*h1(r,l,s,k));k=k+1}
sum(xx)
}
gg1=function(r,l,s,x,k){
xx=0
m=0
while (m<=(s-l-1)) {xx=c(xx,choose(s-l-1,m)*(-1)^(s-l-1-m)*
((h-n*b+f11(r)+f12(s)+(l-r-k)*y[r]+(s-l+k+1)*y[s])^(1-(n-s+r+g+2))-
(h-n*b+f11(r)+f12(s)+(l-r-k)*y[r]+(k+m+1)*x+(s-l-m)*y[s])^(1-(n-s+r+g+2)))/(k+m+1)/(1-(n-s+r+g+2)));m=m+1}
sum(xx)
}
ggg=function(r,l,s,x){
xx=0
k=0
while (k<=(l-r-1)) {xx=c(xx,choose(l-r-1,k)*(-1)^k*gg1(r,l,s,x,k));k=k+1}
sum(xx)
}
hh1=function(r,l,s,k){
xx=0
m=0
while (m<=(s-l-1)) {xx=c(xx,choose(s-l-1,m)*(-1)^(s-l-1-m)*
((h-n*b+f11(r)+f12(s)+(l-r-k)*y[r]+(s-l+k+1)*y[s])^(1-(n-s+r+g+2))-
(h-n*b+f11(r)+f12(s)+(l-r+m+1)*y[r]+(s-l-m)*y[s])^(1-(n-s+r+g+2)))/(k+m+1)/(1-(n-s+r+g+2)));m=m+1}
sum(xx)
}
hhh=function(r,l,s){
xx=0
k=0
while (k<=(l-r-1)) {xx=c(xx,choose(l-r-1,k)*(-1)^k*hh1(r,l,s,k));k=k+1}
sum(xx)
}
A=array (0,c(n-2,n-1,n))
B=array (0,c(n-2,n-1,n))
C=array (0,c(n-2,2*(n-1),n))
for (p in 1:w){
x<-rexp(n,ss)+mm
y<-sort(x)
for (r in c(1,2,3,5,8,9)){
for (s in c(13,14,15,17,19,20)){
for (l in 10:12){     
yhat=0
yhat=gg(r,l,s)/hh(r,l,s)
A[r,l,s]=A[r,l,s]+yhat
B[r,l,s]=B[r,l,s]+(yhat-y[l])^2
fff1=function(x)  ggg(r,l,s,x)/hhh(r,l,s)-(1-a/2)
fff2=function(x)  ggg(r,l,s,x)/hhh(r,l,s)-(a/2)
Intlow=uniroot(fff1, lower=y[r] , upper=y[s])$root
Intup=uniroot(fff2, lower=y[r] , upper=y[s])$root
C[r,2*(l-1),s]=C[r,2*(l-1),s]+Intlow
C[r,2*(l-1)+1,s]=C[r,2*(l-1)+1,s]+Intup
}
}
}
}
for (r in c(1,2,3,5,8,9)){
for (s in c(13,14,15,17,19,20)){
for (l in 10:12){ 

A[r,l,s]=A[r,l,s]/w
B[r,l,s]=B[r,l,s]/w
C[r,2*(l-1),s]=C[r,2*(l-1),s]/w
C[r,2*(l-1)+1,s]=C[r,2*(l-1)+1,s]/w
}
}

}
round(A,4)
round(B,4)
round(C,4)
D=array (,c((n-2),n,(n-1)))
E=array (,c((n-2),n,(n-1)))
F=array (,c((n-2),2*(n-1),n-1))
for (r in c(1,2,3,5,8,9)){
for (s in c(13,14,15,17,19,20)){
for (l in 10:12){ 
D[r,s,l]=A[r,l,s]  
E[r,s,l]=B[r,l,s]
F[r,2*(s-1)-1,l]=C[r,2*(l-1),s]
F[r,2*(s-1),l]=C[r,2*(l-1)+1,s]
}
}
}
round(D,4)
round(E,4)
round(F,4)

 

Last Updated on Friday, 16 May 2014 21:59