# simul params N=100 nsim=10000 #input position, set & mass # 1: empty 0.05 # 2: 1 0.1 # 3: 2 0.05 # 4: 1,3 0.1 # 5: 1,4 0.1 # 6: 2,4 0.15 # 7: 1,2,4 0.1 # 8: 1,4,5 0.05 # 9: 2,4,5 0.05 # 10:1,2,3,4 0.15 # 11:2,3,4,5 0.05 # 12:1,2,3,4,5 0.05 # output position & sets # 1:empty # 2:1 # 3:2 # 4:3 # 5:4 # 6:5 # 7:1,2 # 8:1,3 # 9:1,4 # 10:1,5 # 11:2,3 # 12:2,4 # 13:2,5 # 14:3,4 # 15:3,5 # 16:4,5 # 17:1,2,3 # 18:1,2,4 # 19:1,2,5 # 20:1,3,4 # 21:1,3,5 # 22:1,4,5 # 23:2,3,4 # 24:2,3,5 # 25:2,4,5 # 26:3,4,5 # 27:1,2,3,4 # 28:1,2,3,5 # 29:1,2,4,5 # 30:1,3,4,5 # 31:2,3,4,5 # 32:1,2,3,4,5 # mass function p=rep(.05, 12) p[2]=.1 p[4]=.1 p[5]=.1 p[6]=.15 p[7]=.1 p[10]=.15 K=2.5 w=matrix(0, nrow=32, ncol=12) w[2,]=c(0,1,0,1,1,0,1,1,0,1,0,1) w[3,]=c(0,0,1,0,0,1,1,0,1,1,1,1) w[4,]=c(0,0,0,1,0,0,0,0,0,1,1,1) w[5,]=c(0,0,0,0,1,1,1,1,1,1,1,1) w[6,]=c(0,0,0,0,0,0,0,1,1,0,1,1) w[7,]=c(0,1,1,1,1,1,2,0,1,2,1,2) w[8,]=c(0,1,0,2,1,0,1,0,0,2,1,2) w[9,]=c(0,1,0,1,2,1,2,2,1,2,1,2) w[10,]=c(0,1,0,1,1,0,1,2,1,1,1,2) w[11,]=c(0,0,1,1,0,1,1,0,1,2,2,2) w[12,]=c(0,0,1,0,1,2,2,1,2,2,2,2) w[13,]=c(0,0,1,0,0,1,1,1,2,1,2,2) w[14,]=c(0,0,0,1,1,1,1,1,1,2,2,2) w[15,]=c(0,0,0,1,0,0,0,1,1,1,2,2) w[16,]=c(0,0,0,0,1,1,1,2,2,1,2,2) w[17,]=c(0,1,1,2,1,1,2,1,1,3,2,3) w[18,]=c(0,1,1,1,2,2,3,2,2,3,2,3) w[19,]=c(0,1,1,1,1,1,2,1,1,2,2,3) w[20,]=c(0,1,0,2,2,1,1,2,1,3,2,3) w[21,]=c(0,1,0,2,1,0,1,2,0,2,2,3) w[22,]=c(0,1,0,0,2,1,2,3,2,2,2,3) w[23,]=c(0,0,1,1,1,2,2,1,2,3,3,3) w[24,]=c(0,0,1,1,0,1,1,1,2,2,3,3) w[25,]=c(0,0,1,0,1,2,2,2,3,2,3,3) w[26,]=c(0,0,0,1,1,1,1,2,2,2,3,3) w[27,]=c(0,1,1,2,2,2,3,2,2,4,3,4) w[28,]=c(0,1,1,2,2,1,2,2,2,3,3,4) w[29,]=c(0,1,1,1,2,2,3,3,3,3,3,4) w[30,]=c(0,1,0,2,2,1,2,3,2,3,3,4) w[31,]=c(0,0,1,1,1,2,2,2,3,3,4,4) w[32,]=c(0,1,1,2,2,2,3,3,3,4,4,5) erg_cp=matrix(0, nrow=nsim+1, ncol=32) erg_pl=matrix(0, nrow=nsim+1, ncol=32) erg_pl2=matrix(0, nrow=nsim+1, ncol=32) erg_plmin=matrix(0, nrow=nsim+1, ncol=32) erg_bel=matrix(0, nrow=nsim+1, ncol=32) erg_belmin=matrix(0, nrow=nsim+1, ncol=32) mean_cp=1:32 mean_pl=1:32 mean_pl2=1:32 mean_plmin=1:32 mean_bel=1:32 mean_belmin=1:32 cp <- function(p, w, K) { r=rep(0, 32) size = c(0,1,1,2,2,2,3,3,3,4,4,5) for (i in 2:32) { for (j in 1:12) r[i] = r[i] + p[j]*w[i,j] } K=0 for (j in 1:12) K= K + size[j]*p[j] r= r / K return (r) } pl <- function(p) { r=rep(1-p[1], 32) r[1]=0 r[2]=p[2]+p[4]+p[5]+p[7]+p[8]+p[10]+p[12] r[3]=p[3]+p[6]+p[7]+p[9]+p[10]+p[11]+p[12] r[4]=p[4]+p[10]+p[11]+p[12] r[5]=p[5]+p[6]+p[8]+p[9]+p[10]+p[11]+p[12] r[6]=p[8]+p[9]+p[11]+p[12] r[7]=r[7]-p[9] r[8]=r[8]-p[3]-p[6]-p[9] r[9]=r[7] r[10]=r[10]-p[3]-p[6] r[11]=r[11]-p[2]-p[5]-p[7]-p[8] r[12]=r[12]-p[2]-p[4]-p[5] r[13]=r[12] r[14]=p[4]+p[6]+p[8]+p[10]+p[11]+p[12] r[15]=p[4]+p[8]+p[9]+p[10]+p[11]+p[12] r[16]=r[16]-p[2]-p[3]-p[4] r[20]=r[20]-p[3] r[21]=r[10] r[22]=r[21] r[23]=r[23]-p[2] r[24]=r[24]-p[2]-p[5] r[25]=r[24] r[26]=r[26]-p[2]-p[3] r[30]=r[20] r[31]=r[23] return (r) } pl2 <- function(p) { r=rep(1-p[1]-p[2]-p[3], 32) r[1]=0 r[2]=p[2]+p[4]+p[5]+p[7]+p[8]+p[10]+p[12] r[3]=p[3]+p[6]+p[7]+p[9]+p[10]+p[11]+p[12] r[4]=p[4]+p[10]+p[11]+p[12] r[5]=p[5]+p[6]+p[8]+p[9]+p[10]+p[11]+p[12] r[6]=p[8]+p[9]+p[11]+p[12] r[7]=p[7]+p[10]+p[12] r[8]=p[4]+p[10]+p[12] r[9]=p[4]+p[7]+p[8]+p[10]+p[12] r[10]=p[8]+p[12] r[11]=p[10]+p[11]+p[12] r[12]=p[6]+p[7]+p[9]+p[10]+p[11]+p[12] r[13]=p[9]+p[11]+p[12] r[14]=r[11] r[15]=p[11]+p[12] r[16]=r[6] r[17]=p[4]+p[7]+p[10]+p[11]+p[12] r[18]=p[5]+p[6]+p[7]+p[8]+p[9]+p[10]+p[11]+p[12] r[19]=p[7]+p[8]+p[9]+p[10]+p[11]+p[12] r[20]=p[4]+p[5]+p[7]+p[8]+p[10]+p[11]+p[12] r[21]=p[4]+p[8]+p[10]+p[11]+p[12] r[22]=p[5]+p[7]+p[8]+p[9]+p[10]+p[11]+p[12] r[23]=p[6]+p[7]+p[9]+p[10]+p[11]+p[12] r[24]=p[9]+p[10]+p[11]+p[12] r[25]=p[6]+p[7]+p[8]+p[9]+p[10]+p[11]+p[12] r[26]=p[8]+p[9]+p[10]+p[11]+p[12] return (r) } plmin <- function(p) { r=rep(p[12], 32) r[1]=0 r[2]=r[2]+p[2]+p[4]+p[5]+p[7]+p[8]+p[10] r[3]=r[3]+p[3]+p[6]+p[7]+p[9]+p[10]+p[11] r[4]=r[4]+p[4]+p[10]+p[11] r[5]=r[5]+p[5]+p[6]+p[8]+p[9]+p[10]+p[11] r[6]=r[6]+p[8]+p[9]+p[11] r[7]=r[7]+p[7]+p[10] r[8]=r[8]+p[4]+p[10] r[9]=r[9]+p[5]+p[7]+p[8]+p[10] r[10]=r[10]+p[8] r[11]=r[11]+p[10]+p[11] r[12]=r[12]+p[6]+p[7]+p[9]+p[10]+p[11] r[13]=r[13]+p[9]+p[11] r[14]=r[14]+p[10]+p[11] r[15]=r[15]+p[11] r[16]=r[16]+p[8]+p[9]+p[11] r[17]=r[17]+p[10] r[18]=r[18]+p[10]+p[7] r[20]=r[17] r[22]=r[22]+p[8] r[23]=r[23]+p[10]+p[11] r[24]=r[15] r[25]=r[25]+p[9]+p[11] r[26]=r[15] r[27]=r[17] r[31]=r[15] return (r) } bel <- function(p) { r=rep(0, 32) r[2]=p[2] r[3]=p[3] r[7]=p[2]+p[3] r[8]=p[2]+p[4] r[9]=p[2]+p[5] r[10]=r[2] r[11]=r[3] r[12]=p[3]+p[6] r[13]=r[3] r[17]=p[2]+p[3]+p[4] r[18]=p[2]+p[3]+p[5]+p[6]+p[7] r[19]=r[7] r[20]=p[2]+p[4]+p[5] r[21]=r[8] r[22]=p[2]+p[5]+p[8] r[23]=r[12] r[24]=r[3] r[25]=p[3]+p[6]+p[9] r[27]=p[2]+p[3]+p[4]+p[5]+p[6]+p[7]+p[10] r[28]=r[17] r[29]=p[2]+p[3]+p[6]+p[7]+p[9] r[30]=p[2]+p[4]+p[8] r[31]=p[3]+p[6]+p[9]+p[11] r[32]=1-p[1] return (r) } belmin <- function(p) { r=rep(0, 32) r[1]=p[1] r[2]=p[2] r[3]=p[3] r[8]=p[4] r[9]=p[5] r[12]=p[6] r[18]=p[7] r[22]=p[8] r[25]=p[9] r[27]=p[10] r[31]=p[11] r[32]=p[12] return (r) } erg_cp[1,]=cp(p, w, K) erg_pl[1,]=pl(p) erg_pl2[1,]=pl2(p) erg_plmin[1,]=plmin(p) erg_bel[1,]=bel(p) erg_belmin[1,]=belmin(p) zuf=rmultinom(nsim, N, p)/N for (i in 1:nsim) { erg_cp[i+1,]=cp(zuf[,i], w, K) erg_pl[i+1,]=pl(zuf[,i]) erg_pl2[i+1,]=pl2(zuf[,i]) erg_plmin[i+1,]=plmin(zuf[,i]) erg_bel[i+1,]=bel(zuf[,i]) erg_belmin[i+1,]=belmin(zuf[,i]) } for (j in 1:32) { mean_cp[j]=mean(erg_cp[,j]) mean_pl[j]=mean(erg_pl[,j]) mean_pl2[j]=mean(erg_pl2[,j]) mean_plmin[j]=mean(erg_plmin[,j]) mean_bel[j]=mean(erg_bel[,j]) mean_belmin[j]=mean(erg_belmin[,j]) } plot(2:32,mean_pl[2:32], type="l", ylim=c(0, 1.2), xlab="31 sets of 2^Q (sets with one element left; Q right)", ylab="bel / pl") lines(2:32, mean_pl2[2:32], col="grey66") lines(2:32, mean_plmin[2:32], col="blue") lines(2:32, mean_bel[2:32], col="red") lines(2:32, mean_belmin[2:32], col="green") lines(2:32, mean_cp[2:32], lty=2) legend(0, 1.24, c("pl", "pl-min", "bel", "bel-min", "pl(k=2)", "cp"), col = c("black", "blue", "red", "green", "grey66", "black"), lty = c(1,1,1,1,1,2) ) mean_x=data.frame(mean_pl, mean_plmin, mean_bel, mean_belmin, mean_pl2, mean_cp) write.table(mean_x,file="X:/meanvalN100.txt",quote=FALSE,sep=",",col.names=NA) # table structure of meanval..-type tables # rows 1 - 32 the output sets (see top) # cols as descibed in the table q_cp=matrix(0, 7, 32) q_pl=matrix(0, 7, 32) q_pl2=matrix(0, 7, 32) q_plmin=matrix(0, 7, 32) q_bel=matrix(0, 7, 32) q_belmin=matrix(0, 7, 32) for (j in 1:32) { q_cp[,j]=quantile(erg_cp[,j], probs = c(.025, 0.05, .25, 0.5, 0.75, .95, .975)) q_pl[,j]=quantile(erg_pl[,j], probs = c(.025, 0.05, .25, 0.5, 0.75, .95, .975)) q_pl2[,j]=quantile(erg_pl[,j], probs = c(.025, 0.05, .25, 0.5, 0.75, .95, .975)) q_plmin[,j]=quantile(erg_plmin[,j], probs = c(.025, 0.05, .25, 0.5, 0.75, .95, .975)) q_bel[,j]=quantile(erg_bel[,j], probs = c(.025, 0.05, .25, 0.5, 0.75, .95, .975)) q_belmin[,j]=quantile(erg_belmin[,j], probs = c(.025, 0.05, .25, 0.5, 0.75, .95, .975)) } write.table(q_cp,file="X:/q_cpN100.txt",quote=FALSE,sep=",",col.names=NA) write.table(q_pl,file="X:/q_plN100.txt",quote=FALSE,sep=",",col.names=NA) write.table(q_pl2,file="X:/q_pl2N100.txt",quote=FALSE,sep=",",col.names=NA) write.table(q_plmin,file="X:/q_plminN100.txt",quote=FALSE,sep=",",col.names=NA) write.table(q_bel,file="X:/q_belN100.txt",quote=FALSE,sep=",",col.names=NA) write.table(q_belmin,file="X:/q_belminN100.txt",quote=FALSE,sep=",",col.names=NA) # table structure of q_..-type tables # cols 1 - 32 the output sets (see top) # row 1: .025-Quantile # row 2: .05-Quantile # row 3: .25-Quantile # row 4: .5-Quantile # row 5: .75-Quantile # row 6: .95-Quantile # row 7: .975-Quantile