###################################################################### # # # Membership Function Sensitivity of Descriptive # # Statistics in Fuzzy-Set Relations # # # # Author: Alrik Thiem # # # ###################################################################### set.seed(10) X <- rnorm(30) # 30 draws from normal round(X <- X + abs(min(X)), 3) # positivize data (tau.e <- as.vector(quantile(X, 0.05))) # exclusion anchor (tau.i <- as.vector(quantile(X, 0.95))) # inclusion anchor (tau.l <- tau.e + 1/4 * (tau.i - tau.e)) # low crossover anchor (tau.m <- tau.e + 1/2 * (tau.i - tau.e)) # medium crossover anchor (tau.h <- tau.e + 3/4 * (tau.i - tau.e)) # high crossover anchor ######################## # MEMBERSHIP FUNCTIONS # ######################## mu.g <- function(X, p, tau.e, tau.c, tau.i){ ifelse (X <= tau.e, 0, ifelse((tau.e < X) & (X <= tau.c), 1/2 * ((tau.e - X)/(tau.e - tau.c))^p, ifelse((tau.c < X) & (X <= tau.i), 1 - 1/2 * ((tau.i - X)/(tau.i - tau.c))^p, 1))) } mu.l <- function(X, tau.e, tau.c, tau.i){ ifelse(X < tau.c, 1 / (1 + exp(-((X - tau.c) * ((-log(19))/(tau.e - tau.c))))), 1 / (1 + exp(-((X - tau.c) * ((log(19))/(tau.i - tau.c)))))) } round(matrix(c(fsvals.ln.l <- mu.g(X, 1, tau.e, tau.l, tau.i), fsvals.ln.m <- mu.g(X, 1, tau.e, tau.m, tau.i), fsvals.ln.h <- mu.g(X, 1, tau.e, tau.h, tau.i)), nrow = 30, ncol = 3, dimnames = list(c(), c("tau.l", "tau.m", "tau.h"))), 2) round(matrix(c(fsvals.qu.l <- mu.g(X, 2, tau.e, tau.l, tau.i), fsvals.qu.m <- mu.g(X, 2, tau.e, tau.m, tau.i), fsvals.qu.h <- mu.g(X, 2, tau.e, tau.h, tau.i)), nrow = 30, ncol = 3, dimnames = list(c(), c("tau.l", "tau.m", "tau.h"))), 2) round(matrix(c(fsvals.ro.l <- mu.g(X, 0.5, tau.e, tau.l, tau.i), fsvals.ro.m <- mu.g(X, 0.5, tau.e, tau.m, tau.i), fsvals.ro.h <- mu.g(X, 0.5, tau.e, tau.h, tau.i)), nrow = 30, ncol = 3, dimnames = list(c(), c("tau.l", "tau.m", "tau.h"))), 2) round(matrix(c(fsvals.lg.l <- mu.l(X, tau.e, tau.l, tau.i), fsvals.lg.m <- mu.l(X, tau.e, tau.m, tau.i), fsvals.lg.h <- mu.l(X, tau.e, tau.h, tau.i)), nrow = 30, ncol = 3, dimnames = list(c(), c("tau.l", "tau.m", "tau.h"))), 2) ################### # BIVARIATE PLOTS # ################### fsvals.mat.ln <- cbind(fsvals.ln.l, fsvals.ln.m, fsvals.ln.h) fsvals.mat.qu <- cbind(fsvals.qu.l, fsvals.qu.m, fsvals.qu.h) fsvals.mat.ro <- cbind(fsvals.ro.l, fsvals.ro.m, fsvals.ro.h) fsvals.mat.lg <- cbind(fsvals.lg.l, fsvals.lg.m, fsvals.lg.h) par(mfcol = c(2, 2), mgp = c(2.1, 0.4, 0), tcl = -0.3, mar = c(3.5, 3.5, 2, 1.5), las = 1) matplot(X, fsvals.mat.ln, type = c("p"), pch = 20, col = gray(c(0.85, 0.5, 0.15)), cex = 1.25, main = expression(paste("2a) ", bolditalic(mu)[lin])), xlab = "Base Variable X", ylab = expression(paste("Membership in ", bold(A)))) segments(tau.e, -1, tau.e, mu.g(tau.e, 1, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(tau.l, -1, tau.l, mu.g(tau.l, 1, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(tau.m, -1, tau.m, mu.g(tau.m, 1, tau.e, tau.m, tau.i), lty = 2, col = gray(0.5)) segments(tau.h, -1, tau.h, mu.g(tau.h, 1, tau.e, tau.h, tau.i), lty = 2, col = gray(0.5)) segments(tau.i, -1, tau.i, mu.g(tau.i, 1, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(min(X) - 1, mu.g(tau.e, 1, tau.e, tau.l, tau.i), tau.e, mu.g(tau.e, 1, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(min(X) - 1, mu.g(tau.h, 1, tau.e, tau.h, tau.i), tau.h, mu.g(tau.h, 1, tau.e, tau.h, tau.i), lty = 2, col = gray(0.5)) segments(min(X) - 1, mu.g(tau.i, 1, tau.e, tau.l, tau.i), tau.i, mu.g(tau.i, 1, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) text(min(X) - 0.2, mu.g(tau.e, 1, tau.e, tau.l, tau.i) + 0.075, labels = expression(alpha[0](tau[5])), pos = 4) text(min(X) - 0.2, mu.g(tau.l, 1, tau.e, tau.l, tau.i) + 0.05, labels = expression(alpha[0.5](tau[c])), pos = 4) text(tau.l, mu.g(tau.l, 1, tau.e, tau.l, tau.i) + 0.05, labels = "l", pos = 2, cex = 0.8) text(tau.m, mu.g(tau.l, 1, tau.e, tau.l, tau.i) + 0.05, labels = "m", pos = 2, cex = 0.8) text(tau.h + 0.5, mu.g(tau.l, 1, tau.e, tau.l, tau.i) + 0.05, labels = "h", pos = 2, cex = 0.8) text(min(X) - 0.2, mu.g(tau.i, 1, tau.e, tau.l, tau.i) - 0.1, labels = expression(alpha[1](tau[95])), pos = 4) matplot(X, fsvals.mat.ro, type = c("p"), pch = 20, col = gray(c(0.85,0.5,0.15)), cex = 1.25, main = expression(paste("2c) ", bolditalic(mu)[root])), xlab = "Base Variable X", ylab = expression(paste("Membership in ", bold(A)))) segments(tau.e, -1, tau.e, mu.g(tau.e, 0.5, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(tau.l, -1, tau.l, mu.g(tau.l, 0.5, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(tau.m, -1, tau.m, mu.g(tau.m, 0.5, tau.e, tau.m, tau.i), lty = 2, col = gray(0.5)) segments(tau.h, -1, tau.h, mu.g(tau.h, 0.5, tau.e, tau.h, tau.i), lty = 2, col = gray(0.5)) segments(tau.i, -1, tau.i, mu.g(tau.i, 0.5, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(min(X) - 1, mu.g(tau.e, 0.5, tau.e, tau.l, tau.i), tau.e, mu.g(tau.e, 0.5, tau.l, tau.e, tau.i), lty = 2, col = gray(0.5)) segments(min(X) - 1, mu.g(tau.h, 0.5, tau.e, tau.h, tau.i), tau.h, mu.g(tau.h, 0.5, tau.e, tau.h, tau.i), lty = 2, col = gray(0.5)) segments(min(X) - 1, mu.g(tau.i, 0.5, tau.e, tau.l, tau.i), tau.i, mu.g(tau.i, 0.5, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) text(min(X) - 0.2, mu.g(tau.e, 0.5, tau.e, tau.l, tau.i) + 0.075, labels = expression(alpha[0](tau[5])), pos = 4) text(min(X) - 0.2, mu.g(tau.l, 0.5, tau.e, tau.l, tau.i) + 0.05, labels = expression(alpha[0.5](tau[c])), pos = 4) text(tau.l, mu.g(tau.l, 0.5, tau.e, tau.l, tau.i) + 0.05, labels = "l", pos = 2, cex = 0.8) text(tau.m + 0.1, mu.g(tau.l, 0.5, tau.e, tau.l, tau.i) + 0.05, labels = "m", pos = 2, cex = 0.8) text(tau.h + 0.5, mu.g(tau.l, 0.5, tau.e, tau.l, tau.i) + 0.05, labels = "h", pos = 2, cex = 0.8) text(min(X) - 0.2, mu.g(tau.i, 0.5, tau.e, tau.l, tau.i) - 0.1, labels = expression(alpha[1](tau[95])), pos = 4) matplot(X, fsvals.mat.qu, type = c("p"), pch = 20, col = gray(c(0.85,0.5,0.15)), cex = 1.25, main = expression(paste("2b) ", bolditalic(mu)[quad])), xlab = "Base Variable X", ylab = expression(paste("Membership in ", bold(A)))) segments(tau.e, -1, tau.e, mu.g(tau.e, 2, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(tau.l, -1, tau.l, mu.g(tau.l, 2, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(tau.m, -1, tau.m, mu.g(tau.m, 2, tau.e, tau.m, tau.i), lty = 2, col = gray(0.5)) segments(tau.h, -1, tau.h, mu.g(tau.h, 2, tau.e, tau.h, tau.i), lty = 2, col = gray(0.5)) segments(tau.i, -1, tau.i, mu.g(tau.i, 2, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(min(X) - 1, mu.g(tau.e, 2, tau.e, tau.l, tau.i), tau.e, mu.g(tau.e, 2, tau.l, tau.e, tau.i), lty = 2, col = gray(0.5)) segments(min(X) - 1, mu.g(tau.h, 2, tau.e, tau.h, tau.i), tau.h, mu.g(tau.h, 2, tau.e, tau.h, tau.i), lty = 2, col = gray(0.5)) segments(min(X) - 1, mu.g(tau.i, 2, tau.e, tau.l, tau.i), tau.i, mu.g(tau.i, 2, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) text(min(X) - 0.2, mu.g(tau.e, 2, tau.e, tau.l, tau.i) + 0.075, labels = expression(alpha[0](tau[5])), pos = 4) text(min(X) - 0.2, mu.g(tau.l, 2, tau.e, tau.l, tau.i) + 0.05, labels = expression(alpha[0.5](tau[c])), pos = 4) text(tau.l, mu.g(tau.l, 2, tau.e, tau.l, tau.i) + 0.05, labels = "l", pos = 2, cex = 0.8) text(tau.m, mu.g(tau.l, 2, tau.e, tau.l, tau.i) + 0.05, labels = "m", pos = 2, cex = 0.8) text(tau.h + 0.5, mu.g(tau.l, 2, tau.e, tau.l, tau.i) + 0.05, labels = "h", pos = 2, cex = 0.8) text(min(X) - 0.2, mu.g(tau.i, 2, tau.e, tau.l, tau.i) - 0.1, labels = expression(alpha[1](tau[95])), pos = 4) matplot(X, fsvals.mat.lg, type = c("p"), pch = 20, col = gray(c(0.85,0.5,0.15)), cex = 1.25, main = expression(paste("2d) ", bolditalic(mu)[log])), xlab = "Base Variable X", ylab = expression(paste("Membership in ", bold(A)))) segments(tau.e, -1, tau.e, mu.l(tau.e, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(tau.i, -1, tau.i, mu.l(tau.i, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(tau.l, -1, tau.l, mu.l(tau.l, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(tau.m, -1, tau.m, mu.l(tau.l, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(tau.h, -1, tau.h, mu.l(tau.l, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(min(X) - 1, mu.l(tau.e, tau.e, tau.l, tau.i), tau.e, mu.l(tau.e, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(min(X) - 1, mu.l(tau.h, tau.e, tau.h, tau.i), tau.h, mu.l(tau.h, tau.e, tau.h, tau.i), lty = 2, col = gray(0.5)) segments(min(X) - 1, mu.l(tau.i, tau.e, tau.l, tau.i), tau.i, mu.l(tau.i, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) text(min(X) - 0.2, mu.l(tau.e, tau.e, tau.l, tau.i) + 0.075, labels = expression(alpha[0](tau[5])), pos = 4) text(min(X) - 0.2, mu.l(tau.m, tau.e, tau.m, tau.i) + 0.05, labels = expression(alpha[0.5](tau[c])), pos = 4) text(tau.l, mu.l(tau.m, tau.e, tau.m, tau.i) + 0.05, labels = "l", pos = 2, cex = 0.8) text(tau.m, mu.l(tau.m, tau.e, tau.m, tau.i) + 0.05, labels = "m", pos = 2, cex = 0.8) text(tau.h + 0.5, mu.l(tau.m, tau.e, tau.m, tau.i) + 0.05, labels = "h", pos = 2, cex = 0.8) text(min(X) - 0.2, mu.l(tau.i, tau.e, tau.h, tau.i) - 0.1, labels = expression(alpha[1](tau[95])), pos = 4) ############################# # GENERATE DATA FOR OUTCOME # ############################# set.seed(1234) r <- runif(30, 0, 1) max.vals.l <- pmax(fsvals.ln.l, fsvals.qu.l, fsvals.ro.l, fsvals.lg.l) max.vals.m <- pmax(fsvals.ln.m, fsvals.qu.m, fsvals.ro.m, fsvals.lg.m) max.vals.h <- pmax(fsvals.ln.h, fsvals.qu.h, fsvals.ro.h, fsvals.lg.h) outvals <- function(x, r){ ifelse((0 < x) & (x < 1), x + r * (1 - x), 1) } fsvals.ln.out.l <- fsvals.qu.out.l <- fsvals.ro.out.l <- fsvals.lg.out.l <- outvals(max.vals.l, r) fsvals.ln.out.m <- fsvals.qu.out.m <- fsvals.ro.out.m <- fsvals.lg.out.m <- outvals(max.vals.m, r) fsvals.ln.out.h <- fsvals.qu.out.h <- fsvals.ro.out.h <- fsvals.lg.out.h <- outvals(max.vals.h, r) ################### # COVERAGE SCORES # ################### coverage <- function(C, O){ sum(pmin(C, O)) / sum(O) } round(matrix(c( coverage.ln.l <- coverage(fsvals.ln.l, fsvals.ln.out.l), coverage.qu.l <- coverage(fsvals.qu.l, fsvals.qu.out.l), coverage.ro.l <- coverage(fsvals.ro.l, fsvals.ro.out.l), coverage.lg.l <- coverage(fsvals.lg.l, fsvals.lg.out.l), coverage.ln.m <- coverage(fsvals.ln.m, fsvals.ln.out.m), coverage.qu.m <- coverage(fsvals.qu.m, fsvals.qu.out.m), coverage.ro.m <- coverage(fsvals.ro.m, fsvals.ro.out.m), coverage.lg.m <- coverage(fsvals.lg.m, fsvals.lg.out.m), coverage.ln.h <- coverage(fsvals.ln.h, fsvals.ln.out.h), coverage.qu.h <- coverage(fsvals.qu.h, fsvals.qu.out.h), coverage.ro.h <- coverage(fsvals.ro.h, fsvals.ro.out.h), coverage.lg.h <- coverage(fsvals.lg.h, fsvals.lg.out.h)), nrow = 4, ncol = 3, dimnames = list(c("ln", "qu", "ro", "lg"), c("tau.l", "tau.m", "tau.h"))), 3) ########################### # CONDITION-OUTCOME PLOTS # ########################### X11() par(mfcol = c(1, 3), tcl = -0.3, las = 1, cex.lab = 1.25, cex.main = 2) plot(fsvals.ln.l, fsvals.ln.out.l, type = "n", xlim = c(-0.006, 1.006), ylim = c(-0.006, 1.006), main = expression(paste("3a) ", tau[l])), xlab = expression(paste("Condition ", bold(C))), ylab = expression(paste("Outcome ", bold(O)))) points(fsvals.ln.l, fsvals.ln.out.l, pch = 19, cex = 1.75) points(fsvals.qu.l, fsvals.qu.out.l, pch = 21, col = grey(0.5), bg = grey(0.5), cex = 1.75) points(fsvals.ro.l, fsvals.ro.out.l, pch = 21, col = grey(0.75), bg = grey(0.75), cex = 1.75) points(fsvals.lg.l, fsvals.lg.out.l, pch = 21, cex = 1.75) abline(0, 1) text(1, 0.35, labels = expression(paste(omega(mu[lin])%~~%0.792)), cex = 1.5, pos = 2) text(1, 0.26, labels = expression(paste(omega(mu[quad])%~~%0.864)), cex = 1.5, pos = 2) text(1, 0.17, labels = expression(paste(omega(mu[root])%~~%0.725)), cex = 1.5, pos = 2) text(1, 0.08, labels = expression(paste(omega(mu[log])%~~%0.82,"0")), cex = 1.5, pos = 2) rug(c(fsvals.ln.l, fsvals.qu.l, fsvals.ro.l, fsvals.lg.l), ticksize = 0.05, side = 1, lwd = 0.1) rug(fsvals.ln.out.l, ticksize = -0.05, side = 4, lwd = 0.1) points(0.55, 0.351, pch = 19, cex = 1.75) points(0.55, 0.261, pch = 21, col = grey(0.5), bg = grey(0.5), cex = 1.75) points(0.55, 0.171, pch = 21, col = grey(0.75), bg = grey(0.75), cex = 1.75) points(0.55, 0.081, pch = 21, cex = 1.75) plot(fsvals.ln.m, fsvals.ln.out.m, xlim = c(-0.006, 1.006), ylim = c(-0.006, 1.006), pch = 19, cex = 1.75, main = expression(paste("3b) ", tau[m])), xlab = expression(paste("Condition ", bold(C))), ylab = expression(paste("Outcome ", bold(O)))) points(fsvals.qu.m, fsvals.qu.out.m, pch = 21, col = grey(0.5), bg = grey(0.5), cex = 1.75) points(fsvals.ro.m, fsvals.ro.out.m, pch = 21, col = grey(0.75), bg = grey(0.75), cex = 1.75) points(fsvals.lg.m, fsvals.lg.out.m, pch = 21, cex = 1.75) abline(0, 1) text(1, 0.35, labels = expression(paste(omega(mu[lin])%~~%0.713)), cex = 1.5, pos = 2) text(1, 0.26, labels = expression(paste(omega(mu[quad])%~~%0.729)), cex = 1.5, pos = 2) text(1, 0.17, labels = expression(paste(omega(mu[root])%~~%0.696)), cex = 1.5, pos = 2) text(1, 0.08, labels = expression(paste(omega(mu[log])%~~%0.716)), cex = 1.5, pos = 2) rug(c(fsvals.ln.m, fsvals.qu.m, fsvals.ro.m, fsvals.lg.m), ticksize = 0.05, side = 1, lwd = 0.1) rug(fsvals.ln.out.m, ticksize = -0.05, side = 4, lwd = 0.1) points(0.55, 0.351, pch = 19, cex = 1.75) points(0.55, 0.261, pch = 21, col = grey(0.5), bg = grey(0.5), cex = 1.75) points(0.55, 0.171, pch = 21, col = grey(0.75), bg = grey(0.75), cex = 1.75) points(0.55, 0.081, pch = 21, cex = 1.75) plot(fsvals.ln.h, fsvals.ln.out.h, xlim = c(-0.006, 1.006), ylim = c(-0.006, 1.006), pch = 19, cex = 1.75, main = expression(paste("3c) ", tau[h])), xlab = expression(paste("Condition ", bold(C))), ylab = expression(paste("Outcome ", bold(O)))) points(fsvals.qu.h, fsvals.qu.out.h, pch = 21, col = grey(0.5), bg = grey(0.5), cex = 1.75) points(fsvals.ro.h, fsvals.ro.out.h, pch = 21, col = grey(0.75), bg = grey(0.75), cex = 1.75) points(fsvals.lg.h, fsvals.lg.out.h, pch = 21, cex = 1.75) abline(0, 1) text(1, 0.35, labels = expression(paste(omega(mu[lin])%~~%0.596)), cex = 1.5, pos = 2) text(1, 0.26, labels = expression(paste(omega(mu[quad])%~~%0.529)), cex = 1.5, pos = 2) text(1, 0.17, labels = expression(paste(omega(mu[root])%~~%0.647)), cex = 1.5, pos = 2) text(1, 0.08, labels = expression(paste(omega(mu[log])%~~%0.561)), cex = 1.5, pos = 2) rug(c(fsvals.ln.h, fsvals.qu.h, fsvals.ro.h, fsvals.lg.h), ticksize = 0.05, side = 1, lwd = 0.1) rug(fsvals.ln.out.h, ticksize = -0.05, side = 4, lwd = 0.1) points(0.55, 0.351, pch = 19, cex = 1.75) points(0.55, 0.261, pch = 21, col = grey(0.5), bg = grey(0.5), cex = 1.75) points(0.55, 0.171, pch = 21, col = grey(0.75), bg = grey(0.75), cex = 1.75) points(0.55, 0.081, pch = 21, cex = 1.75) ################################## # SIMULATION OF CONDITION VALUES # ################################## # draw 1000 random samples of size 30 from normal distribution set.seed(10) X.sim <- matrix(rnorm(30000), 30, 1000) # make all samples values positive X.sim <- apply(X.sim, 2, function(X.sim) X.sim + abs(min(X.sim))) # determine exclusion and inclusion thresholds for each sample tau.e.sim <- apply(X.sim, 2, function(X.sim){quantile(X.sim, 0.05)}) tau.i.sim <- apply(X.sim, 2, function(X.sim){quantile(X.sim, 0.95)}) # functions for low, medium and high crossover thresholds cross.l <- function(tau.e, tau.i){ tau.e + 1/4 * (tau.i - tau.e) } cross.m <- function(tau.e, tau.i){ tau.e + 1/2 * (tau.i - tau.e) } cross.h <- function(tau.e, tau.i){ tau.e + 3/4 * (tau.i - tau.e) } # sample crossover thresholds tau.l.sim <- cross.l(tau.e.sim, tau.i.sim) tau.m.sim <- cross.m(tau.e.sim, tau.i.sim) tau.h.sim <- cross.h(tau.e.sim, tau.i.sim) # initialize 4 (functions) x 3 (thresholds) = 12 matrices for fuzzy-set values fsvals.ln.sim.l <- fsvals.qu.sim.l <- fsvals.ro.sim.l <- fsvals.lg.sim.l <- fsvals.ln.sim.m <- fsvals.qu.sim.m <- fsvals.ro.sim.m <- fsvals.lg.sim.m <- fsvals.ln.sim.h <- fsvals.qu.sim.h <- fsvals.ro.sim.h <- fsvals.lg.sim.h <- matrix(numeric(30000), 30, 1000) # apply the respective membership function over the columns of each of the 12 matrices # to produce condition set membership scores for(i in 1:1000){ fsvals.ln.sim.l[, i] <- mu.g(X.sim[, i], 1, tau.e.sim[i], tau.l.sim[i], tau.i.sim[i]) fsvals.qu.sim.l[, i] <- mu.g(X.sim[, i], 2, tau.e.sim[i], tau.l.sim[i], tau.i.sim[i]) fsvals.ro.sim.l[, i] <- mu.g(X.sim[, i], 0.5, tau.e.sim[i], tau.l.sim[i], tau.i.sim[i]) fsvals.lg.sim.l[, i] <- mu.l(X.sim[, i], tau.e.sim[i], tau.l.sim[i], tau.i.sim[i]) fsvals.ln.sim.m[, i] <- mu.g(X.sim[, i], 1, tau.e.sim[i], tau.m.sim[i], tau.i.sim[i]) fsvals.qu.sim.m[, i] <- mu.g(X.sim[, i], 2, tau.e.sim[i], tau.m.sim[i], tau.i.sim[i]) fsvals.ro.sim.m[, i] <- mu.g(X.sim[, i], 0.5, tau.e.sim[i], tau.m.sim[i], tau.i.sim[i]) fsvals.lg.sim.m[, i] <- mu.l(X.sim[, i], tau.e.sim[i], tau.m.sim[i], tau.i.sim[i]) fsvals.ln.sim.h[, i] <- mu.g(X.sim[, i], 1, tau.e.sim[i], tau.h.sim[i], tau.i.sim[i]) fsvals.qu.sim.h[, i] <- mu.g(X.sim[, i], 2, tau.e.sim[i], tau.h.sim[i], tau.i.sim[i]) fsvals.ro.sim.h[, i] <- mu.g(X.sim[, i], 0.5, tau.e.sim[i], tau.h.sim[i], tau.i.sim[i]) fsvals.lg.sim.h[, i] <- mu.l(X.sim[, i], tau.e.sim[i], tau.h.sim[i], tau.i.sim[i]) } ################################ # SIMULATION OF OUTCOME VALUES # ################################ # find the maximum values for each crossover threshold max.vals.sim.l <- pmax(fsvals.ln.sim.l, fsvals.qu.sim.l, fsvals.ro.sim.l, fsvals.lg.sim.l) max.vals.sim.m <- pmax(fsvals.ln.sim.m, fsvals.qu.sim.m, fsvals.ro.sim.m, fsvals.lg.sim.m) max.vals.sim.h <- pmax(fsvals.ln.sim.h, fsvals.qu.sim.h, fsvals.ro.sim.h, fsvals.lg.sim.h) # generate an outcome value matrix for each of the condition value matrices fsvals.ln.out.sim.l <- fsvals.qu.out.sim.l <- fsvals.ro.out.sim.l <- fsvals.lg.out.sim.l <- matrix(outvals(max.vals.sim.l, r), nrow = 30, ncol = 1000) fsvals.ln.out.sim.m <- fsvals.qu.out.sim.m <- fsvals.ro.out.sim.m <- fsvals.lg.out.sim.m <- matrix(outvals(max.vals.sim.m, r), nrow = 30, ncol = 1000) fsvals.ln.out.sim.h <- fsvals.qu.out.sim.h <- fsvals.ro.out.sim.h <- fsvals.lg.out.sim.h <- matrix(outvals(max.vals.sim.h, r), nrow = 30, ncol = 1000) # compute coverage vectors (length 1000) for each of the 12 matrices coverage.ln.sim.l <- apply(fsvals.ln.sim.l, 2, sum) / apply(fsvals.ln.out.sim.l, 2, sum) coverage.qu.sim.l <- apply(fsvals.qu.sim.l, 2, sum) / apply(fsvals.qu.out.sim.l, 2, sum) coverage.ro.sim.l <- apply(fsvals.ro.sim.l, 2, sum) / apply(fsvals.ro.out.sim.l, 2, sum) coverage.lg.sim.l <- apply(fsvals.lg.sim.l, 2, sum) / apply(fsvals.lg.out.sim.l, 2, sum) coverage.ln.sim.m <- apply(fsvals.ln.sim.m, 2, sum) / apply(fsvals.ln.out.sim.m, 2, sum) coverage.qu.sim.m <- apply(fsvals.qu.sim.m, 2, sum) / apply(fsvals.qu.out.sim.m, 2, sum) coverage.ro.sim.m <- apply(fsvals.ro.sim.m, 2, sum) / apply(fsvals.ro.out.sim.m, 2, sum) coverage.lg.sim.m <- apply(fsvals.lg.sim.m, 2, sum) / apply(fsvals.lg.out.sim.m, 2, sum) coverage.ln.sim.h <- apply(fsvals.ln.sim.h, 2, sum) / apply(fsvals.ln.out.sim.h, 2, sum) coverage.qu.sim.h <- apply(fsvals.qu.sim.h, 2, sum) / apply(fsvals.qu.out.sim.h, 2, sum) coverage.ro.sim.h <- apply(fsvals.ro.sim.h, 2, sum) / apply(fsvals.ro.out.sim.h, 2, sum) coverage.lg.sim.h <- apply(fsvals.lg.sim.h, 2, sum) / apply(fsvals.lg.out.sim.h, 2, sum) ######################################### # PLOT: COVERAGE SCORES FROM SIMULATION # ######################################### # bind together coverage vectors for each crossover threshold coverage.mat.l <- cbind(coverage.ln.sim.l, coverage.qu.sim.l, coverage.ro.sim.l, coverage.lg.sim.l) coverage.mat.m <- cbind(coverage.ln.sim.m, coverage.qu.sim.m, coverage.ro.sim.m, coverage.lg.sim.m) coverage.mat.h <- cbind(coverage.ln.sim.h, coverage.qu.sim.h, coverage.ro.sim.h, coverage.lg.sim.h) # produce kernel density plots of coverage scores X11() par(mfcol = c(1, 3), tcl = -0.3, las = 1, cex.lab = 1.25, cex.main = 2) plot(density(coverage.ln.sim.l), type = "n", main = expression(paste("4a) ", tau[l])), xlab = expression(bold(omega)), ylim = c(0, 20), xlim = c(0.55, 0.95)) segments(mean(coverage.ln.sim.l), 0, mean(coverage.ln.sim.l), 12, col = gray(0.5), lty = 2) segments(mean(coverage.qu.sim.l), 0, mean(coverage.qu.sim.l), 8, col = gray(0.5), lty = 2) segments(mean(coverage.ro.sim.l), 0, mean(coverage.ro.sim.l), 19, col = gray(0.5), lty = 2) segments(mean(coverage.lg.sim.l), 0, mean(coverage.lg.sim.l), 9.8, col = gray(0.5), lty = 2) lines(density(coverage.ln.sim.l), lwd = 2, col = gray(0.5)) lines(density(coverage.qu.sim.l), lwd = 2, col = gray(0)) lines(density(coverage.ro.sim.l), lwd = 2, col = gray(0.75)) lines(density(coverage.lg.sim.l), lwd = 2, col = gray(0.25)) text(0.63, 18, labels = expression(paste(bar(omega)%~~%0.69))) text(0.74, 13, labels = 0.74) text(0.78, 12, labels = 0.77) text(0.82, 10, labels = 0.80) legend(max(coverage.mat.l) - 0.08, 20, bg = gray(0.975), y.intersp = 1.25, x.intersp = 0.75, lty = rep(1, 4), lwd = 4, c(expression(bold(mu)[quad]), expression(bold(mu)[log]), expression(bold(mu)[lin]), expression(bold(mu)[root])), col = c(gray(0), gray(0.25), gray(0.5), gray(0.75)), cex = 1.25) plot(density(coverage.ln.sim.m), type = "n", main = expression(paste("4b) ", tau[m])), xlab = expression(bold(omega)), ylim = c(0, 20), xlim = c(0.35, 0.9)) segments(mean(coverage.ro.sim.m), 0, mean(coverage.ro.sim.m), 15, col = gray(0.5), lty = 2) lines(density(coverage.ln.sim.m), lwd = 2, col = gray(0.5)) lines(density(coverage.qu.sim.m), lwd = 2, col = gray(0)) lines(density(coverage.ro.sim.m), lwd = 2, col = gray(0.75)) lines(density(coverage.lg.sim.m), lwd = 2, col = gray(0.25)) text(0.65, 16.5, labels = expression(paste(bar(omega)%~~%0.65))) legend(max(coverage.mat.m) - 0.1, 20, bg = gray(0.975), y.intersp = 1.25, x.intersp = 0.75, lty = rep(1, 4), lwd = 4, c(expression(bold(mu)[quad]), expression(bold(mu)[log]), expression(bold(mu)[lin]), expression(bold(mu)[root])), col = c(gray(0), gray(0.25), gray(0.5), gray(0.75)), cex = 1.25) plot(density(coverage.ln.sim.h), type = "n", main = expression(paste("4c) ", tau[h])), xlab = expression(bold(omega)), ylim = c(0, 20), xlim = c(0.2, 0.75)) segments(mean(coverage.ln.sim.h), 0, mean(coverage.ln.sim.h), 8.4, col = gray(0.5), lty = 2) segments(mean(coverage.qu.sim.h), 0, mean(coverage.qu.sim.h), 5.6, col = gray(0.5), lty = 2) segments(mean(coverage.ro.sim.h), 0, mean(coverage.ro.sim.h), 13.8, col = gray(0.5), lty = 2) segments(mean(coverage.lg.sim.h), 0, mean(coverage.lg.sim.h), 7.3, col = gray(0.5), lty = 2) lines(density(coverage.ln.sim.h), lwd = 2, col = gray(0.5)) lines(density(coverage.qu.sim.h), lwd = 2, col = gray(0)) lines(density(coverage.ro.sim.h), lwd = 2, col = gray(0.75)) lines(density(coverage.lg.sim.h), lwd = 2, col = gray(0.25)) text(0.42, 6.5, labels = expression(paste(bar(omega)%~~%0.48))) text(0.48, 8, labels = 0.52) text(0.53, 9.5, labels = 0.54) text(0.61, 15.25, labels = 0.61) legend(max(coverage.mat.h) - 0.47, 20, bg = gray(0.975), y.intersp = 1.25, x.intersp = 0.75, lty = rep(1, 4), lwd = 4, c(expression(bold(mu)[quad]), expression(bold(mu)[log]), expression(bold(mu)[lin]), expression(bold(mu)[root])), col = c(gray(0), gray(0.25), gray(0.5), gray(0.75)), cex = 1.25) ######################################## # COVERAGE SCORES OVER CROSSOVER RANGE # ######################################## # function for range of crossover thresholds cross.range <- function(tau.e, tau.i){ tau.e + seq(1/10, 9/10, by = 1/20) * (tau.i - tau.e) } # initialize matrix of 17 (crossover thresholds) x 1000 (samples) = 17000 elements tau.c.range.sim <- matrix(numeric(17000), 17, 1000) for(i in 1:1000){ tau.c.range.sim[, i] <- cross.range(tau.e.sim[i], tau.i.sim[i]) } # initialize four matrices of dimension: # rows : 17 crossover thresholds x 30 condition values = 510 # columns: 1000 samples fsvals.ln.range.sim <- fsvals.qu.range.sim <- fsvals.ro.range.sim <- fsvals.lg.range.sim <- matrix(numeric(510000), 510, 1000) # apply the respective membership function 17 times over each column, every time # for the next 30 values to produce condition set membership scores; # may take some seconds to process for(p in 1:1000) { for(i in 1:17) { fsvals.ln.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], 1, tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) fsvals.qu.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], 2, tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) fsvals.ro.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], 0.5, tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) fsvals.lg.range.sim[(-29+30*i):(i*30), p] <- mu.l(X.sim[, p], tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) } } # find the maximum condition values max.vals.range.sim <- pmax(fsvals.ln.range.sim, fsvals.qu.range.sim, fsvals.ro.range.sim, fsvals.lg.range.sim) # generate an outcome value matrix for each of the four condition value matrices fsvals.ln.out.range.sim <- fsvals.qu.out.range.sim <- fsvals.ro.out.range.sim <- fsvals.lg.out.range.sim <- matrix(numeric(510000), 510, 1000) for(p in 1:1000){ for(i in 1:17){ fsvals.ln.out.range.sim[(-29+30*i):(i*30), p] <- fsvals.qu.out.range.sim[(-29+30*i):(i*30), p] <- fsvals.ro.out.range.sim[(-29+30*i):(i*30), p] <- fsvals.lg.out.range.sim[(-29+30*i):(i*30), p] <- outvals(max.vals.range.sim[(-29+30*i):(i*30), p], r) } } # initialize a coverage matrix of dimension 17 x 1000 for each membership function coverage.ln.range.sim <- coverage.qu.range.sim <- coverage.ro.range.sim <- coverage.lg.range.sim <- matrix(numeric(17000), 17, 1000) for(p in 1:1000){ for(i in 1:17){ coverage.ln.range.sim[i, p] <- sum(fsvals.ln.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.ln.out.range.sim[(-29+30*i):(i*30), p]) coverage.qu.range.sim[i, p] <- sum(fsvals.qu.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.qu.out.range.sim[(-29+30*i):(i*30), p]) coverage.ro.range.sim[i, p] <- sum(fsvals.ro.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.ro.out.range.sim[(-29+30*i):(i*30), p]) coverage.lg.range.sim[i, p] <- sum(fsvals.lg.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.lg.out.range.sim[(-29+30*i):(i*30), p]) } } # average coverage scores from all 1000 samples for each crossover threshold cov.ln.range.sim.mean <- apply(coverage.ln.range.sim, 1, mean) cov.qu.range.sim.mean <- apply(coverage.qu.range.sim, 1, mean) cov.ro.range.sim.mean <- apply(coverage.ro.range.sim, 1, mean) cov.lg.range.sim.mean <- apply(coverage.lg.range.sim, 1, mean) k <- seq(1/10, 9/10, by = 1/20) X11() par(mar = c(4.5, 4, 0.5, 0.5), las = 1) plot(k, cov.qu.range.sim.mean, type="n", ylim = c(0.36, 0.9), ylab = expression(bar(bolditalic(omega))), xlab = expression(italic(k)), xaxt = "n", las = 1) abline(h = seq(0.4, 0.8, by = 0.1), col = gray(0.95), lwd = 0.5) segments(k[4], 0.32, k[4], cov.qu.range.sim.mean[4], col = gray(0.5), lty = 2) segments(k[9], 0.32, k[9], cov.ln.range.sim.mean[9], col = gray(0.5), lty = 2) segments(k[14], 0.32, k[14], cov.ro.range.sim.mean[14], col = gray(0.5), lty = 2) segments(k[4], cov.qu.range.sim.mean[4], 0, cov.qu.range.sim.mean[4], col = gray(0.5), lty = 2) segments(k[9], cov.ln.range.sim.mean[9], 0, cov.ln.range.sim.mean[9], col = gray(0.5), lty = 2) segments(k[14], cov.ro.range.sim.mean[14], 0, cov.ro.range.sim.mean[14], col = gray(0.5), lty = 2) text(k[4] + 0.02, 0.36, labels = expression(tau[l])) text(k[9] + 0.02, 0.36, labels = expression(tau[m])) text(k[14] + 0.02, 0.36, labels = expression(tau[h])) axis(1, at = k[seq(1, 17, by = 2)], lab = seq(0.1, 0.9, by = 0.1)) lines(k, cov.qu.range.sim.mean, col = gray(0), lwd = 5) lines(k, cov.lg.range.sim.mean, col = gray(0.25), lwd = 5) lines(k, cov.ln.range.sim.mean, col = gray(0.5), lwd = 5) lines(k, cov.ro.range.sim.mean, col = gray(0.75), lwd = 5) arrows(k[1], 0.9, k[9], angle = 20, code = 3) arrows(k[9], 0.9, k[17], angle = 20, code = 3) text(k[5], 0.87, labels = expression(paste(bar(omega)[quad], " > ", bar(omega)[log], " > ", bar(omega)[lin], " > ", bar(omega)[root]))) text(k[13], 0.87, labels = expression(paste(bar(omega)[root], " > ", bar(omega)[lin], " > ", bar(omega)[log], " > ", bar(omega)[quad]))) legend(0.76, 0.83, bg = gray(0.975), y.intersp = 1.25, x.intersp = 1, lty = rep(1, 4), lwd = 5, c(expression(bold(mu)[quad]), expression(bold(mu)[log]), expression(bold(mu)[lin]), expression(bold(mu)[root])), col = c(gray(0), gray(0.25), gray(0.5), gray(0.75))) ##################################################### # COVERAGE SCORES OVER CROSSOVER AND FUNCTION RANGE # ##################################################### fsvals.1.range.sim <- fsvals.2.range.sim <- fsvals.3.range.sim <- fsvals.4.range.sim <- fsvals.5.range.sim <- fsvals.6.range.sim <- fsvals.7.range.sim <- fsvals.8.range.sim <- fsvals.9.range.sim <- fsvals.10.range.sim <- fsvals.11.range.sim <- matrix(numeric(510000), 510, 1000) exp.range <- c(1/seq(2, 1, -0.2), seq(1.2, 2, 0.2)) for(p in 1:1000) { for(i in 1:17) { fsvals.1.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], exp.range[1], tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) fsvals.2.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], exp.range[2], tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) fsvals.3.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], exp.range[3], tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) fsvals.4.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], exp.range[4], tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) fsvals.5.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], exp.range[5], tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) fsvals.6.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], exp.range[6], tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) fsvals.7.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], exp.range[7], tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) fsvals.8.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], exp.range[8], tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) fsvals.9.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], exp.range[9], tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) fsvals.10.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], exp.range[10], tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) fsvals.11.range.sim[(-29+30*i):(i*30), p] <- mu.g(X.sim[, p], exp.range[11], tau.e.sim[p], tau.c.range.sim[i, p], tau.i.sim[p]) } } max.vals.p.range.sim <- pmax(fsvals.1.range.sim, fsvals.2.range.sim, fsvals.3.range.sim, fsvals.4.range.sim, fsvals.5.range.sim, fsvals.6.range.sim, fsvals.7.range.sim, fsvals.8.range.sim, fsvals.9.range.sim, fsvals.10.range.sim, fsvals.11.range.sim) fsvals.1.out.range.sim <- fsvals.2.out.range.sim <- fsvals.3.out.range.sim <- fsvals.4.out.range.sim <- fsvals.5.out.range.sim <- fsvals.6.out.range.sim <- fsvals.7.out.range.sim <- fsvals.8.out.range.sim <- fsvals.9.out.range.sim <- fsvals.10.out.range.sim <- fsvals.11.out.range.sim <- matrix(numeric(510000), 510, 1000) for(p in 1:1000){ for(i in 1:17){ fsvals.1.out.range.sim[(-29+30*i):(i*30), p] <- fsvals.2.out.range.sim[(-29+30*i):(i*30), p] <- fsvals.3.out.range.sim[(-29+30*i):(i*30), p] <- fsvals.4.out.range.sim[(-29+30*i):(i*30), p] <- fsvals.5.out.range.sim[(-29+30*i):(i*30), p] <- fsvals.6.out.range.sim[(-29+30*i):(i*30), p] <- fsvals.7.out.range.sim[(-29+30*i):(i*30), p] <- fsvals.8.out.range.sim[(-29+30*i):(i*30), p] <- fsvals.9.out.range.sim[(-29+30*i):(i*30), p] <- fsvals.10.out.range.sim[(-29+30*i):(i*30), p] <- fsvals.11.out.range.sim[(-29+30*i):(i*30), p] <- outvals(max.vals.range.sim[(-29+30*i):(i*30), p], r) } } coverage.1.range.sim <- coverage.2.range.sim <- coverage.3.range.sim <- coverage.4.range.sim <- coverage.5.range.sim <- coverage.6.range.sim <- coverage.7.range.sim <- coverage.8.range.sim <- coverage.9.range.sim <- coverage.10.range.sim <- coverage.11.range.sim <- matrix(numeric(17000), 17, 1000) for(p in 1:1000){ for(i in 1:17){ coverage.1.range.sim[i, p] <- sum(fsvals.1.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.1.out.range.sim[(-29+30*i):(i*30), p]) coverage.2.range.sim[i, p] <- sum(fsvals.2.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.2.out.range.sim[(-29+30*i):(i*30), p]) coverage.3.range.sim[i, p] <- sum(fsvals.3.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.3.out.range.sim[(-29+30*i):(i*30), p]) coverage.4.range.sim[i, p] <- sum(fsvals.4.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.4.out.range.sim[(-29+30*i):(i*30), p]) coverage.5.range.sim[i, p] <- sum(fsvals.5.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.5.out.range.sim[(-29+30*i):(i*30), p]) coverage.6.range.sim[i, p] <- sum(fsvals.6.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.6.out.range.sim[(-29+30*i):(i*30), p]) coverage.7.range.sim[i, p] <- sum(fsvals.7.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.7.out.range.sim[(-29+30*i):(i*30), p]) coverage.8.range.sim[i, p] <- sum(fsvals.8.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.8.out.range.sim[(-29+30*i):(i*30), p]) coverage.9.range.sim[i, p] <- sum(fsvals.9.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.9.out.range.sim[(-29+30*i):(i*30), p]) coverage.10.range.sim[i, p] <- sum(fsvals.10.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.10.out.range.sim[(-29+30*i):(i*30), p]) coverage.11.range.sim[i, p] <- sum(fsvals.11.range.sim[(-29+30*i):(i*30), p]) / sum(fsvals.11.out.range.sim[(-29+30*i):(i*30), p]) } } dat <- matrix(numeric(561), 187, 3, dimnames = list(1:187, c("coverage", "p", "k"))) dat[, 1] <- rbind( cov.1.range.sim.mean <- apply(coverage.1.range.sim, 1, mean), cov.2.range.sim.mean <- apply(coverage.2.range.sim, 1, mean), cov.3.range.sim.mean <- apply(coverage.3.range.sim, 1, mean), cov.4.range.sim.mean <- apply(coverage.4.range.sim, 1, mean), cov.5.range.sim.mean <- apply(coverage.5.range.sim, 1, mean), cov.6.range.sim.mean <- apply(coverage.6.range.sim, 1, mean), cov.7.range.sim.mean <- apply(coverage.7.range.sim, 1, mean), cov.8.range.sim.mean <- apply(coverage.8.range.sim, 1, mean), cov.9.range.sim.mean <- apply(coverage.9.range.sim, 1, mean), cov.10.range.sim.mean <- apply(coverage.10.range.sim, 1, mean), cov.11.range.sim.mean <- apply(coverage.11.range.sim, 1, mean) ) dat[, 2] <- rep(exp.range, 17) dat[, 3] <- rep(k, 1, each = 11) dat <- as.data.frame(dat) col.regions = rev(grey.colors(100, start = 0, end = 1)) library(lattice) levelplot(coverage ~ k * p, data = dat, col.regions = col.regions, xlab = expression(italic(k)), ylab = expression(italic(p)), at = seq(0.35, 0.9, 0.025), colorkey = list(tick.number = 9), contour = T, labels = F, scales = list(x = list(at = seq(0.1, 0.9, 0.1)), y = list(at = exp.range, labels = c(expression(2.0^(-1)), expression(1.8^(-1)), expression(1.6^(-1)), expression(1.4^(-1)), expression(1.2^(-1)), 1.0,1.2,1.4,1.6,1.8,2.0)))) ################################### # MAXIMUM DIFFERENCES IN COVERAGE # ################################### x.pos <- function(tau.e, tau.c){ tau.e + ((-tau.e^3 + 3*tau.e^2*tau.c - 3*tau.e*tau.c^2 + tau.c^3)^(1/3))/(2*2^(1/3)) } (x.pos.tau.l <- x.pos(tau.e, tau.l)) (x.pos.tau.m <- x.pos(tau.e, tau.m)) (x.pos.tau.h <- x.pos(tau.e, tau.h)) x.neg <- function(tau.c, tau.i){ tau.i + (-(-(tau.c^3 - 3*tau.c^2*tau.i + 3*tau.c*tau.i^2 - tau.i^3))^(1/3))/(2*2^(1/3)) } (x.neg.tau.l <- x.neg(tau.l, tau.i)) (x.neg.tau.m <- x.neg(tau.m, tau.i)) (x.neg.tau.h <- x.neg(tau.h, tau.i)) x.diff <- function(x,tau.e,tau.c,tau.i){ mu.g(x,1/2,tau.e,tau.c,tau.i) - mu.g(x,2,tau.e,tau.c,tau.i) } (pos.diff.tau.l <- x.diff(x.pos.tau.l, tau.e, tau.l, tau.i)) (pos.diff.tau.m <- x.diff(x.pos.tau.m, tau.e, tau.m, tau.i)) (pos.diff.tau.h <- x.diff(x.pos.tau.h, tau.e, tau.h, tau.i)) (neg.diff.tau.l <- x.diff(x.neg.tau.l, tau.e, tau.l, tau.i)) (neg.diff.tau.m <- x.diff(x.neg.tau.m, tau.e, tau.m, tau.i)) (neg.diff.tau.h <- x.diff(x.neg.tau.h, tau.e, tau.h, tau.i)) fsvals.diff.cov <- cbind(fsvals.qu.l, fsvals.qu.m, fsvals.qu.h, fsvals.ro.l, fsvals.ro.m, fsvals.ro.h) curve(mu.g(x, 1/2, tau.e, tau.l, tau.i),0,3.5, xaxt = "n", yaxt = "n", xlab = "Base Variable X", ylab = expression(paste("Membership in ", bold(C))), col = gray(0.75)) curve(mu.g(x, 1/2, tau.e, tau.m, tau.i),0,3.5, xaxt = "n", yaxt = "n", col = gray(0.45), add = TRUE) curve(mu.g(x, 1/2, tau.e, tau.h, tau.i),0,3.5, xaxt = "n", yaxt = "n", col = gray(0), add = TRUE) curve(mu.g(x, 2, tau.e, tau.l, tau.i),0,3.5, xaxt = "n", yaxt = "n", col = gray(0.75), add = TRUE) curve(mu.g(x, 2, tau.e, tau.m, tau.i),0,3.5, xaxt = "n", yaxt = "n", col = gray(0.45), add = TRUE) curve(mu.g(x, 2, tau.e, tau.h, tau.i),0,3.5, xaxt = "n", yaxt = "n", col = gray(0), add = TRUE) segments(-0.5, 0.5, max(X)+0.5, 0.5, lty = 2, col = gray(0.5)) segments(x.pos.tau.l, -0.5, x.pos.tau.l, mu.g(x.pos.tau.l, 1/2, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(x.pos.tau.m, -0.5, x.pos.tau.m, mu.g(x.pos.tau.m, 1/2, tau.e, tau.m, tau.i), lty = 2, col = gray(0.5)) segments(x.pos.tau.h, -0.5, x.pos.tau.h, mu.g(x.pos.tau.h, 1/2, tau.e, tau.h, tau.i), lty = 2, col = gray(0.5)) segments(x.neg.tau.l, 1.5, x.neg.tau.l, mu.g(x.neg.tau.l, 1/2, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(x.neg.tau.m, 1.5, x.neg.tau.m, mu.g(x.neg.tau.m, 1/2, tau.e, tau.m, tau.i), lty = 2, col = gray(0.5)) segments(x.neg.tau.h, 1.5, x.neg.tau.h, mu.g(x.neg.tau.h, 1/2, tau.e, tau.h, tau.i), lty = 2, col = gray(0.5)) segments(-0.5, mu.g(x.pos.tau.h, 1/2, tau.e, tau.h, tau.i), x.pos.tau.h, mu.g(x.pos.tau.h, 1/2, tau.e, tau.h, tau.i), lty = 2, col = gray(0.5)) segments(-0.5, mu.g(x.pos.tau.h, 2, tau.e, tau.h, tau.i), x.pos.tau.h, mu.g(x.pos.tau.h, 2, tau.e, tau.h, tau.i), lty = 2, col = gray(0.5)) segments(3.75, mu.g(x.neg.tau.l, 1/2, tau.e, tau.l, tau.i), x.neg.tau.l, mu.g(x.neg.tau.l, 1/2, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) segments(3.75, mu.g(x.neg.tau.l, 2, tau.e, tau.l, tau.i), x.neg.tau.l, mu.g(x.neg.tau.l, 2, tau.e, tau.l, tau.i), lty = 2, col = gray(0.5)) matplot(X, fsvals.diff.cov, xaxt = "n", yaxt = "n", type = c("p"), pch = 20, col = gray(c(0.75, 0.45, 0)), cex = 1.5, add = TRUE) axis(1, at=c(0,tau.e,0.57,0.85,1.14,2,2.5,3,3.5), labels = c(0,expression(tau[e[phantom(l)]]),expression(x["[+]"]^tau[l]),expression(x["[+]"]^tau[m]), expression(x["[+]"]^tau[h]),2,2.5,3,3.5)) axis(2, at=c(mu.g(x.pos.tau.l, 2, tau.e, tau.l, tau.i),mu.g(x.pos.tau.l, 1/2, tau.e, tau.l, tau.i),0.6,0.8,1), labels = c(expression(mu[quad](x["[+]"]^tau[c])),expression(mu[root](x["[+]"]^tau[c])),0.6,0.8,1)) axis(3, at=c(0,0.5,1,1.5,2.29,2.58,2.86,tau.i,3.5), labels = c(0,0.5,1,1.5,expression(x["[-]"]^tau[l]),expression(x["[-]"]^tau[m]),expression(x["[-]"]^tau[h]), expression(tau[i[phantom(l)]]),3.5)) axis(4, at=c(0,0.2,0.4,mu.g(x.neg.tau.h, 1/2, tau.e, tau.h, tau.i),mu.g(x.neg.tau.h, 2, tau.e, tau.h, tau.i)), labels = c(0,0.2,0.4,expression(mu[root](x["[-]"]^tau[c])),expression(mu[quad](x["[-]"]^tau[c])))) text(min(X) - 0.1, 0.55, labels = expression(alpha[0.5](tau[c])), pos = 4) text(tau.l, 0.47, labels = expression(tau[l]), pos = 4) text(tau.m, 0.47, labels = expression(tau[m]), pos = 4) text(tau.h, 0.47, labels = expression(tau[h]), pos = 4) arrows(0, mu.g(x.pos.tau.l, 2, tau.e, tau.l, tau.i), 0, mu.g(x.pos.tau.l, 1/2, tau.e, tau.l, tau.i), length = 0.1, angle = 30, code = 3) arrows(3.5, mu.g(x.neg.tau.l, 2, tau.e, tau.l, tau.i), 3.5, mu.g(x.neg.tau.l, 1/2, tau.e, tau.l, tau.i), length = 0.1, angle = 30, code = 3) text(0.05, 0.13, labels = expression(paste(Delta, c["[+](root)"])), pos = 4, srt = 90) text(3.3, 0.74, labels = expression(paste(Delta, c["[-](root)"])), pos = 4, srt = 90)