cut <- function(poss,s){ # s-cut of a triangular possibility distribution u<- poss$a+s*(poss$b - poss$a) v<- poss$b + (1-s)*(poss$c - poss$b) return(list(u=u,v=v)) } inter <- function(I1,I2){ # intersection of two intervals u<-max(I1$u,I2$u) v<-min(I1$v,I2$v) if(u>v){ I<-NULL # returns NULL if intersection is empty } else { I<-list(u=u,v=v) } return(I) } dempster<-function(poss1,poss2,N){ # Dempster's combination of 2 triangular possibility distributions F<-matrix(0,N,2) i<-0 while(i < N){ s1<-runif(1) s2<-runif(2) I1<-cut(poss1,s1) I2<-cut(poss2,s2) I12<-inter(I1,I2) if(!is.null(I12)){ i<-i+1 F[i,]<-c(I12$u,I12$v) } } return(F) } # Example poss1<-list(a=0,b=2,c=3) poss2<-list(a=1,b=1.5,c=3.5) B<-dempster(poss1,poss2,100000) # Lower and upper cdfs Finf<-ecdf(B[,2]) Fsup<-ecdf(B[,1]) plot(Finf,xlim=c(1,3)) lines(Fsup) # Contour function x<-seq(1,3,0.01) n<-length(x) pl<-rep(0,n) for(i in 1:n){ pl[i]<-mean((B[,1]<= x[i]) & (B[,2]>=x[i])) } plot(x,pl,type="l")