### Implementation in R of the software described in: ## Torra, V., Narukawa, Y. (2016) Numerical integration for the ## Choquet integral, Information fusion 31 137-145. ## http://dx.doi.org/10.1016/j.inffus.2016.02.007 ## Torra, V., Narukawa, Y. (2015) Distances on non-additive measures ## using the numerical Choquet integral, USB Proc. MDAI 2015. ### Examples from the above papers and also from FSS, ANOR, and EUSFLAT papers: ## Torra, V., Narukawa, Y., Sugeno, M. (2016, in press) ## On the f-divergence for non-additive measures, Fuzzy sets and ## Systems, in press. ## http://dx.doi.org/10.1016/j.fss.2015.07.006 ## Narukawa, Y., Torra, V., Sugeno, M. (2012, in press still!!) ## Choquet integral with respect to a symmetric fuzzy measure of a function ## on the real line, Annals of Operations Research, in press. ## http://dx.doi.org/10.1007/s10479-012-1166-6 ## Torra, V., Y. Narukawa, M. Sugeno, M. Carlson (2013) Hellinger distance ## for fuzzy measures, Proc. EUSFLAT Conference. ## http://dx.doi.org/10.2991/eusflat.2013.82 ### Notes about the implementation: ## 1. The function INTEGRATE needs a function f ## that must accept a vector of inputs and produce a vector of function ## evaluations at those points. ## I use the Vectorize function for this purpose. ## 2. Integration functions may give an error if the integral is divergent, ## or there are {\em similar} problems. Change limits of the integral. library(stats) ## INCREASING FUNCTIONS =============================================== ## Implemented functions to integrate increasing functions: ## ciIncreasingDLebesgueWithInverse (f, oneInvF=NULL, m, alpha0, alpha1) ## ciIncreasingDProbabilityWithInverse(f, oneInvF=NULL, dp, alpha0, alpha1) ## Function: ## Choquet integral of an increasing function f in the interval ## [alpha0, alpha1] with (or without) known inverse (oneInv) ## and a fuzzy measure that is a distorted Lebesgue with distortion m ## Input: ## f: increasing function to be integrated ## oneInvF: inverse of F (NULL if not available) in ## the interval [alpha0,alpha1] ## m: distortion of the Lebesgue measure ## alpha0, alpha1: interval where to integrate the function f ## Output: ## $value: numerical integral ## $abs.error: error of the integral (supplied by the function integrate) ## Example: (definitions below) ## ciIncreasingDLebesgueWithInverse(function1,NULL,distortion1,0,1) ciIncreasingDLebesgueWithInverse <- function (f, oneInvF=NULL, m, alpha0, alpha1) { if (is.null(oneInvF)) { oneInv <- function (f, l, u) { function (y) { uniroot(function(x) {f(x)-y}, lower=l, upper=u)$root}} oneInvF <- oneInv(f,alpha0,alpha1) } # Function to integrate, evaluable at a SINGLE point beta # Example: integrand(function1,distortion1,0,1)(0.2) integrand <- function (f, m, alpha0, alpha1) { function (beta) { m(alpha1-oneInvF(beta))}} # Vectorial function to integrate for a vector of points # Example: integrandvector(function1,distortion1, 0,1)(c(0.2,0.2)) integrandvector <- function (f, m, alpha0, alpha1) { Vectorize(integrand(f,m,alpha0,alpha1), c("beta"))} # Integral # Example: integral(function1,distortion1,0,1) integral <- integrate(integrandvector(f,m,alpha0,alpha1), lower=f(alpha0), upper=f(alpha1)) noIntegral <- m(alpha1-alpha0)*f(alpha0) res <- list(value = noIntegral + integral$value, abs.error = integral$abs.error) return (res) } ### Example 1 @ (Torra, Narukawa, 2016) ### http://dx.doi.org/10.1016/j.inffus.2016.02.007 ### ### Code R for this example: $f(x)=x^3$, $m(x)=x^2$ ### ### Result: The integral is 0.1. function1 <- function (x) { x*x*x } distortion1 <- function (x) { x*x } inverse1 <- function (x) { x^(1/3)} ciIncreasingDLebesgueWithInverse(function1,inverse1,distortion1,0,1) # $value -> [1] 0.1 # $abs.error -> [1] 2.981058e-08 ciIncreasingDLebesgueWithInverse(function1,NULL,distortion1,0,1) # $value -> [1] 0.1000001 # $abs.error -> [1] 1.042551e-05 ## Function: ## Choquet integral of an increasing function f in the interval ## [alpha0, alpha1] with (or without) known inverse (oneInv) and ## a fuzzy measure that is a distorted probability dp ## Input: ## f: increasing function to be integrated ## oneInvF: inverse of F (NULL if not available) in ## the interval [alpha0,alpha1] ## dp: distorted probability. dp(a,b) is the distorted probability ## of the interval [a,b] ## alpha0, alpha1: interval where to integrate the function f ## Output: ## $value: numerical integral ## $abs.error: error of the integral (supplied by the function integrate) ## Example: (definitions below) ## ciIncreasingDProbabilityWithInverse (function2,oneInvF=NULL,dpmuTSoTND,0,1) ciIncreasingDProbabilityWithInverse <- function(f,oneInvF=NULL,dp,alpha0,alpha1){ # Function to integrate, evaluable at a SINGLE point beta if (is.null(oneInvF)) { oneInv <- function (f, l, u) { function (y) { uniroot(function(x) {f(x)-y}, lower=l, upper=u)$root}} oneInvF <- oneInv(f,alpha0,alpha1) } # Example: integrand(function1,distortion1,0,1)(0.2) integrand <- function (f, dp, alpha0, alpha1) { function (beta) { dp(oneInvF(beta), alpha1)}} # Vectorial function to integrate for a vector of points # Example: integrandvector(function1,distortion1, 0,1)(c(0.2,0.2)) integrandvector <- function (f, dp, alpha0, alpha1) { Vectorize(integrand(f,dp,alpha0,alpha1), c("beta"))} # Integral # Example: integral(function1,distortion1,0,1) integral <- integrate(integrandvector(f,dp,alpha0,alpha1), lower=f(alpha0), upper=f(alpha1)) noIntegral <- dp(alpha0,alpha1)*f(alpha0) res <- list(value = noIntegral + integral$value, abs.error = integral$abs.error) return (res) } ### Example 2 @ (Torra, Narukawa, 2016) ### http://dx.doi.org/10.1016/j.inffus.2016.02.007 ### ### Code R for this example: $f(x)=sqrt(sin(x pi/2))$ ### $mu=m o P$ with $P$ a truncated normal distribution ### Computation of the normal distribution uses the function ### of probability of a normal distribution: ### pnorm(q, mean, sd) ### Result: The integral is 0.83897 ### Probability: Truncated normal distribution truncatedND <- function (a,b) { return ((pnorm(b,0.5,0.1) - pnorm(a,0.5,0.1))/ (pnorm(1,0.5,0.1) - pnorm(0,0.5,0.1))) } ### Distortion: Truncated sigmoid sigmoidComplete <- function (x) { 1/(1+exp(-10*(x-0.5))) } truncatedSigmoid <- function (x) { (sigmoidComplete(x)-sigmoidComplete(0))/ (sigmoidComplete(1)-sigmoidComplete(0)) } ### Distorted probability: dpmuTSoTND <- function (a,b) { truncatedSigmoid(truncatedND(a,b))} # It can be seen that this is not a distorted Lebesgue because # \mu([0,0.2])=9.15484e-05, but \mu([0.4,0.6])=0.866295. # Due to the fact that the function is increasing, any fuzzy measure will suit ### Function and inverse: ## inverse only valid in the [0,1] interval, ## function not invertible in [0,a] for a>1 function2 <- function (x) {sqrt(sin(x*(pi/2)))} inverse2 <- function (y) {(2/pi)*asin(y*y)} ### Calculation of the solution ciIncreasingDProbabilityWithInverse (function2,oneInvF=NULL,dpmuTSoTND,0,1) ## $value: [1] 0.8389763 $abs.error: [1] 6.969475e-05 ciIncreasingDProbabilityWithInverse (function2,inverse2,dpmuTSoTND,0,1) ## $value: [1] 0.8389743 $abs.error: [1] 3.613824e-08 ## Comparison of the integral of this function w.r.t. a Lebesgue measure with ## our implementation and with Mathematica dpmuId <- function (a,b) { return(b-a) } ciIncreasingDProbabilityWithInverse(function2,inverse2, dpmuId,0,1) ciIncreasingDProbabilityWithInverse(function2,NULL, dpmuId,0,1) # The Lebesgue integral w.r.t. the Lebesgue measure is 0.76276 # (using Mathematica) # This implementation in R returns: # $value -> [1] 0.7627599; $abs.error -> [1] 9.882829e-05 # $value -> [1] 0.7627577; $abs.error -> [1] 0.0001186309 ## DECREASING FUNCTIONS =============================================== ## Implemented functions to integrate decreasing functions: ## ciDecreasingDProbabilityWithInverse (f,oneInvF=NULL,dp,alpha0,alpha1) ## Function: ## Choquet integral of a decreasing function f in the interval ## [alpha0, alpha1] with (or without) known inverse (oneInv) and ## a fuzzy measure that is a distorted probability dp ## Input: ## f: decreasing function to be integrated ## oneInvF: inverse of F (NULL if not available) in ## the interval [alpha0,alpha1] ## dp: distorted probability. dp(a,b) is the distorted probability ## of the interval [a,b] ## alpha0, alpha1: interval where to integrate the function f ## Output: ## $value: numerical integral ## $abs.error: error of the integral (supplied by the function integrate) ## Example: (definitions below) ## ciIncreasingDProbabilityWithInverse (function2,oneInvF=NULL,dpmuTSoTND,0,1) ciDecreasingDProbabilityWithInverse <- function(f,oneInvF=NULL,dp,alpha0,alpha1){ # Function to integrate, evaluable at a SINGLE point beta if (is.null(oneInvF)) { oneInv <- function (f, l, u) { function (y) { uniroot(function(x) {f(x)-y}, lower=l, upper=u)$root}} oneInvF <- oneInv(f,alpha0,alpha1) } # Example: integrand(function1,distortion1,0,1)(0.2) integrand <- function (f, dp, alpha0, alpha1) { function (beta) { dp(alpha0, oneInvF(beta))}} # Vectorial function to integrate for a vector of points # Example: integrandvector(function1,distortion1, 0,1)(c(0.2,0.2)) integrandvector <- function (f, dp, alpha0, alpha1) { Vectorize(integrand(f,dp,alpha0,alpha1), c("beta"))} # Integral # Example: integral(function1,distortion1,0,1) integral <- integrate(integrandvector(f,dp,alpha0,alpha1), lower=f(alpha1), upper=f(alpha0)) noIntegral <- dp(alpha0,alpha1)*f(alpha1) res <- list(value = noIntegral + integral$value, abs.error <- integral$abs.error) return (res) } ### Example based on Example 2 @ (Torra, Narukawa, 2016) ### http://dx.doi.org/10.1016/j.inffus.2016.02.007 ### Given f the function in Example 2, ### Define f' = f(1-x) (which is decreasing) and has the same integral in [0,1] function3 <- function(x) ( function2 (1-x)) inverse3 <- function(x) ( 1-inverse2 (x)) ciDecreasingDProbabilityWithInverse (function3,oneInvF=NULL,dpmuTSoTND,0,1) ciDecreasingDProbabilityWithInverse (function3,inverse3,dpmuTSoTND,0,1) ### Example 3 @ (Torra, Narukawa, 2016) ### http://dx.doi.org/10.1016/j.inffus.2016.02.007 ### This proceeds from Example 27 of Proc. EUSFLAT 27 and ### from Example 35 of FSS paper ## Distortion for a=1 dpmuEx35a1 <- function (x) { return(x + x * x / 2) } dpmuEx35mIda1 <- function (a,b) { dpmuEx35a1(b-a)} ## Distortion generic for arbitrary a dpmuEx35 <- function (a) { function (x) { x + a * x * x / 2 }} ## Lebesgue measure (distortion = identity) dpmuId <- function (a,b) { return(b-a) } dpmuEx35mId <- function(a) { function (a1,b1) { dpmuEx35(a)(dpmuId(a1,b1)) } } ## Function to be integrated example35 <- function(a) { function (x) { (sqrt((2/a)*(1-exp(-a*x))) - sqrt(cosh(a*x)))^2} } ##plot(example35(2)) ciDecreasingDProbabilityWithInverse (example35(1), oneInvF=NULL, dpmuEx35mIda1, 0, 1) ciDecreasingDProbabilityWithInverse (example35(1), oneInvF=NULL, dpmuEx35mId(1), 0, 1) ## SINGLE PEAK FUNCTIONS =============================================== ## Implemented functions to integrate single peak functions: ## ciSinglePeakDProbabilityWithInverses (f, inverse.minus, inverse.plus, dp, gamma, alpha0,alpha1) ## Function: ## Choquet integral of a single peak function f in the interval ## [alpha0, alpha1] with (or without) known inverse ## (inverse.minus, inverse.plus) and ## a fuzzy measure that is a distorted probability dp ## gamma is the peak (the maximum) of the function (if known) ## Input: ## f: function to be integrated (with a single maximum in the interval) ## inverse.minus: the inverse of the increasing part ## inverse.plus: the inverse of the decreasing part ## dp: distorted probability. dp(a,b) is the distorted probability ## of the interval [a,b] ## gamma: the maximum in [0,1] (i.e., the peak of the function) ## alpha0, alpha1: interval where to integrate the function f ## Output: ## $value: numerical integral ## $abs.error: error of the integral (supplied by the function integrate) ## Example: (definitions below) ## ciSinglePeakDProbabilityWithInverses <- function(f, inverse.minus=NULL, inverse.plus=NULL, dp, gamma=NULL, alpha0,alpha1){ # Find peak if (is.null(gamma)) { gamma <- optimize(f, interval=c(alpha0, alpha1),maximum=TRUE)$maximum } # Find inverse left part if (is.null(inverse.minus)) { oneInv <- function (f) { function (y, l, u) { if (y<=f(l)) { l } else { uniroot(function(x) {f(x)-y}, lower=l, upper=gamma)$root}}} inverse.minus <- oneInv(f) } # Find inverse right part if (is.null(inverse.plus)) { oneInv <- function (f) { function (y, l, u) { if (y<=f(u)) { u } else { uniroot(function(x) {f(x)-y}, lower=gamma, upper=u)$root}}} inverse.plus <- oneInv(f) } # Example: integrand(function1,distortion1,0,1)(0.2) integrand <- function (f, dp, alpha0, alpha1) { function (beta) { dp(inverse.minus(beta,alpha0, gamma), inverse.plus(beta,gamma, alpha1))}} # Vectorial function to integrate for a vector of points # Example: integrandvector(function1,distortion1, 0,1)(c(0.2,0.2)) integrandvector <- function (f, dp, alpha0, alpha1) { Vectorize(integrand(f,dp,alpha0,alpha1),c("beta"))} # Integral # Example: integral(function1,distortion1,0,1) integral <- integrate(integrandvector(f,dp,alpha0,alpha1), lower=min(f(alpha0),f(alpha1)),upper=f(gamma)) noIntegral <- dp(alpha0,alpha1)*min(f(alpha0),f(alpha1)) res <- list(value = noIntegral + integral$value, abs.error = integral$abs.error) return (res) } ## Example of single peaked function (a maximum) function33 <- function (x) { 2 - (x - 0.5)*(x - 0.5) } inverse33.minus <- function (bb, alpha0, alpha1) { if (bb < function33(alpha0)) { alpha0 } else { (1-sqrt(1-4*(bb-1.75)))/2 }} inverse33.plus <- function (bb, alpha0, alpha1) { if (bb < function33(alpha1)) { alpha1 } else { (1+sqrt(1-4*(bb-1.75)))/2 }} ## It uses the fuzzy measures defined above in Example 2. ciSinglePeakDProbabilityWithInverses(function33,NULL,inverse33.plus,dpmuTSoTND,0.5,0,1) ciSinglePeakDProbabilityWithInverses(function33,inverse33.minus,NULL,dpmuTSoTND,0.5,0,1) ciSinglePeakDProbabilityWithInverses(function33,NULL,NULL,dpmuTSoTND,0.5,0,1) ciSinglePeakDProbabilityWithInverses(function33,NULL,NULL,dpmuTSoTND,NULL,0,1) ## $value --> [1] 1.994087 ; $abs.error --> [1] 5.033194e-06 ciSinglePeakDProbabilityWithInverses(function33,NULL,NULL,dpmuId,NULL,0,1) ## This integral (with Lebesgue measure) is, according to Mathematica: 1.91667 ## SINGLE VALLEY FUNCTIONS =============================================== ## Implemented functions to integrate single valley functions: ## ciSingleValleyDProbabilityWithInverses (f, inverse.minus, inverse.plus, ## dp, prob, gamma, alpha0,alpha1) ## Function: ## Choquet integral of a single valley function f in the interval ## [alpha0, alpha1] with (or without) known inverse ## (inverse.minus, inverse.plus) and ## a fuzzy measure that is a distorted probability dp ## gamma is the valley (the minimum) of the function (if known) ## Input: ## f: function to be integrated (with a single minimum in the interval) ## inverse.minus: the inverse of the decreasing part ## inverse.plus: the inverse of the increasing part ## dp: distorted probability. dp(a,b) is the distorted probability ## of the interval [a,b] ## gamma: the maximum in [0,1] (i.e., the peak of the function) ## alpha0, alpha1: interval where to integrate the function f ## Output: ## $value: numerical integral ## $abs.error: error of the integral (supplied by the function integrate) ## Example: (definitions below) ## ciSingleValleyDProbabilityWithInverses <- function(f, inverse.minus=NULL, inverse.plus=NULL, dp, prob, gamma=NULL, alpha0,alpha1){ if (is.null(gamma)) { gamma <- optimize(f,interval=c(alpha0, alpha1),maximum=FALSE)$minimum } if (is.null(inverse.minus)) { oneInv <- function (f) { function (y, l, u) { if (u<=l) { print(22222); print(l); print(u); } if (y>=f(l)) { l } else { uniroot(function(x) {f(x)-y}, lower=l, upper=u)$root}}} inverse.minus <- oneInv(f) } if (is.null(inverse.plus)) { oneInv <- function (f) { function (y, l, u) { if (u<=l) { print(22222); print(l); print(u); } if (y>=f(u)) { u } else { uniroot(function(x) {f(x)-y}, lower=l, upper=u)$root}}} inverse.plus <- oneInv(f) } # Example: integrand(function1,distortion1,0,1)(0.2) integrand <- function (f, dp, prob, alpha0, alpha1) { function (beta) { dp(prob(alpha0,inverse.minus(beta,alpha0, gamma))+ prob(inverse.plus(beta,gamma, alpha1),alpha1))}} # Vectorial function to integrate for a vector of points # Example: integrandvector(function1,distortion1, 0,1)(c(0.2,0.2)) integrandvector <- function (f, dp, prob, alpha0, alpha1) { Vectorize(integrand(f,dp, prob,alpha0,alpha1),c("beta"))} # Integral # Example: integral(function1,distortion1,0,1) integral <- integrate(integrandvector(f,dp,prob,alpha0,alpha1), lower=f(gamma), upper=max(f(alpha0),f(alpha1))) noIntegral <- dp(prob(alpha0,alpha1))* f(gamma) res <- list(value = noIntegral + integral$value, abs.error = integral$abs.error) return (res) } ### Example of single valley in Section 3.4 based on the expressions of ### Example 3 @ (Torra, Narukawa, 2016) ### Computations of the Hellinger distance for a=2, 3, 4, 5 (case a=1 considered in Example 3). ### This proceeds from Example 27 of Proc. EUSFLAT 27 and ### from Example 35 of FSS paper ## Case a=1 (decreasing function) ciDecreasingDProbabilityWithInverse (example35(1), oneInvF=NULL, dpmuEx35mIda1, 0, 1) ciDecreasingDProbabilityWithInverse (example35(1), oneInvF=NULL, dpmuEx35mId(1), 0, 1) ## Cases a=2, 3, 4, and 5 (single valley function) ciSingleValleyDProbabilityWithInverses(example35(2), inverse.minus=NULL, inverse.plus=NULL, dpmuEx35(2), dpmuId, gamma=NULL, 0,1)$value ciSingleValleyDProbabilityWithInverses(example35(3), inverse.minus=NULL, inverse.plus=NULL, dpmuEx35(3), dpmuId, gamma=NULL, 0,1)$value ciSingleValleyDProbabilityWithInverses(example35(4), inverse.minus=NULL, inverse.plus=NULL, dpmuEx35(4), dpmuId, gamma=NULL, 0,1)$value ciSingleValleyDProbabilityWithInverses(example35(5), inverse.minus=NULL, inverse.plus=NULL, dpmuEx35(5), dpmuId, gamma=NULL, 0,1)$value ## Computation of the Hellinger distance for Cases a=2, 3, 4, and 5 (single valley function) sqrt(0.5*ciSingleValleyDProbabilityWithInverses(example35(2), inverse.minus=NULL, inverse.plus=NULL, dpmuEx35(2), dpmuId, gamma=NULL, 0,1)$value) sqrt(0.5*ciSingleValleyDProbabilityWithInverses(example35(3), inverse.minus=NULL, inverse.plus=NULL, dpmuEx35(3), dpmuId, gamma=NULL, 0,1)$value) sqrt(0.5*ciSingleValleyDProbabilityWithInverses(example35(4), inverse.minus=NULL, inverse.plus=NULL, dpmuEx35(4), dpmuId, gamma=NULL, 0,1)$value) sqrt(0.5*ciSingleValleyDProbabilityWithInverses(example35(5), inverse.minus=NULL, inverse.plus=NULL, dpmuEx35(5), dpmuId, gamma=NULL, 0,1)$value) ## Computations of Figures 4 and 5 for the paper res0150 <- sapply(1:50/10.0, function(x) { sqrt(0.5*ciSingleValleyDProbabilityWithInverses(example35(x), inverse.minus=NULL, inverse.plus=NULL, dpmuEx35(5), dpmuId, gamma=NULL, 0,1)$value)}) setEPS() postscript("fig.numerical.integral.example35.a0150.eps") plot(res0150, main="Example. Distance in terms of a") dev.off() setEPS() postscript("fig.Eusflat27.FSS35.Example a=1.eps") plot(example35(1), main="Example a=1") dev.off() setEPS() postscript("fig.Eusflat27.FSS35.Example a=2.eps") plot(example35(2), main="Example a=2") dev.off() setEPS() postscript("fig.Eusflat27.FSS35.Example a=3.eps") plot(example35(3), main="Example a=3") dev.off() setEPS() postscript("fig.Eusflat27.FSS35.Example a=4.eps") plot(example35(4), main="Example a=4") dev.off() setEPS() postscript("fig.Eusflat27.FSS35.Example a=5.eps") plot(example35(5), main="Example a=5") dev.off() ### Example of single valley: ### Example 36 @ (Torra, Narukawa, Sugeno, 2016) FSS paper ### Computations for the Hellinger distance for p=1,2,3,4,5,6 example36 <- function(p) { function(x) { 1-sqrt(2*p*(p-1))*x^((p-2)/2)+(p*(p-1)/2)*x^(p-2) } } ## plot(example36(1)) # line ## plot(example36(2)) ## plot(example36(3)) # single valley function distLp <- function (p) { function (x) (x^p) } ciSingleValleyDProbabilityWithInverses(example36(1), NULL, NULL, distLp(2), dpmuId, NULL, 0,1) ## Error in integrate(integrandvector(f, dp, prob, alpha0, alpha1), lower = f(gamma), : a limit is missing ## Change lower limit to 0.01: ciSingleValleyDProbabilityWithInverses(example36(1), NULL, NULL, distLp(2), dpmuId, NULL, 0.1,1) ciSingleValleyDProbabilityWithInverses(example36(2), NULL, NULL, distLp(2), dpmuId, NULL, 0,1) # $value -> [1] 0 ; $abs.error -> [1] 0 ciSingleValleyDProbabilityWithInverses(example36(3), NULL, NULL, distLp(2), dpmuId, NULL, 0,1) # $value -> [1] 0.09379493 ; $abs.error -> [1] 9.462158e-06 ciSingleValleyDProbabilityWithInverses(example36(4), NULL, NULL, distLp(2), dpmuId, NULL, 0,1) # $value -> [1] 0.2558993 ; $abs.error -> [1] 9.565143e-06 ciSingleValleyDProbabilityWithInverses(example36(5), NULL, NULL, distLp(2), dpmuId, NULL, 0,1) # $value -> [1] 0.4066559 ; $abs.error -> [1] 0.0001048329 ciSingleValleyDProbabilityWithInverses(example36(6), NULL, NULL, distLp(2), dpmuId, NULL, 0,1) # $value -> [1] 0.5368777 ; $abs.error -> [1] 7.700711e-05 ### Example of single valley: ### Example 38 @ (Torra, Narukawa, Sugeno, 2016) FSS paper ### Computations for the Hellinger distance for p=1 example38 <- function(p) { function (x) { 2*x - 2 * sqrt (2*p*(x^p)) + p * x^(p-1)} } plot(example38(1)) ciSingleValleyDProbabilityWithInverses(example38(1), NULL, NULL, function(x) (x), dpmuId, NULL, 0,1) # $value -> [1] 0.1143923 ; $abs.error -> [1] 4.744442e-05 sqrt(0.5*ciSingleValleyDProbabilityWithInverses(example38(1), NULL, NULL, function(x) (x), dpmuId, NULL, 0,1)$value) ## POLYNOMIALS =============================================== ## Implemented functions to integrate polynomials: ## ciPositivePolynomial (vPol, alfa0, alfa1, fmdx2) ## Function: ## Function that returns a function that is the given polynomial vPol ## Input: ## vPol: vector representing the coefficients of the polynomial ## c(b0,b1,b2,.., bn) represents el polinomi es b0+b1*x+b2*x^2+b3*x^3 + .. + bn*x^n ## Output: ## A function that corresponds to the given polynomial ## Examples: ## vpol <- VectPolynomial(c(10,-60,55,-20,2.5)) ## plot(vpol,xlim=c(0.5,3.5)) ## plot(VectPolynomial(c(10,-60,55,-20,2.5)),xlim=c(0.5,3.5)) VectPolynomial <- function (vPol) { Vectorize(function (x) { sum(vPol * sapply(0:(length(vPol)-1), function(n) { x^n } )) } ,c("x"))} ## Function: ## Function to evaluate a polynomial at a point ## Input: ## vPol: vector representing the coefficients of the polynomial ## c(b0,b1,b2,.., bn) represents el polinomi es b0+b1*x+b2*x^2+b3*x^3 + .. + bn*x^n ## x: the point where to evaluate the polynomial ## Output: ## Example: ## evalPolynomial(c(1,1,2),3) # Output: 22 evalPolynomial <- function (vPol, x) { sum(vPol * sapply(0:(length(vPol)-1), function(n) { x^n })) } ## Function: ## Function to compute the derivative of the polynomial. The derivative of ## c(b0,b1,b2,…,bn) = b0+b1*x+b2*x^2+b3*x^3 + … + bn*x^n ## is ## c(b1,2b2,3b3,…,n bn) = b1 +2*b2*x^1+3*b3*x^2 + …+ n bn*x^{n-1} ## Input: ## vPol: vector representing the coefficients of the polynomial ## c(b0,b1,b2,.., bn) represents el polinomi es b0+b1*x+b2*x^2+b3*x^3 + .. + bn*x^n ## Output: ## another vector representing the derivative of the polynomial represented by vPol ## Example: ## polynomialDerivative(c(6,4,4,5,5)) # -> c(4,8,15,20) (6+4x+4x^2+5x^3+5x^4) -> (4+8x+15x^2+20x^3) ## polynomialDerivative(c(1,3,1)) # ---> c(3,2) (1+3x+x^2) -> (3+2x) ## polynomialDerivative(c(6,4,4,5,5)) ## polynomialDerivative(c(10,-60,55,-20,2.5)) polynomialDerivative <- function (vPol) { 1:(length(vPol)-1) * vPol[2: (length(vPol))] } ## Function: ## Function to return the zeros of a polynomial in a given interval union de limits of the interval ## Input: ## vPol: Polynomial (a vector of coefficients) ## c(b0,b1,b2,…,bn) = b0+b1*x+b2*x^2+b3*x^3 + … + bn*x^n ## [alfa0, alfa1]: interval ## Output: ## a vector that corresponds to the zeros in a given interval \cup the limits of the interval ## Example: ## zerosPol(c(6,4,4,5,5),-4,4) ## zerosPol(c(10,-60,55,-20,2.5),0,4) ## Notes: Polyroot calcula els zeros d’un polinomi ## Si c(1,2,1) = c(1,2,1,0,0) ## 1 + 2x^1 + 1x^2 ## Donat c(b0,b1,b2,…) el polinomi es b0+b1*x+b2*x^2+b3*x^3 + … zerosPol <- function (vPol, alfa0, alfa1) { B0 <- polyroot(vPol) # print(B0) B1 <- Re(B0[abs(Im(B0))<0.000001]) # print(B1) B2 <- (round(B1* 1000000))/1000000 # print(B2) B3 <- B2[B2>=alfa0 & B2 <= alfa1] # print(B3) if (!(is.element(alfa0, B3))) { B3 <- c(alfa0, B3) } if (!(is.element(alfa1, B3))) { B3 <- c(B3, alfa1) } return (B3) } ## Function: ## Function to return the maxima and minima of a polynomial in a given interval ## Input: ## vPol: Polynomial (a vector of coefficients) ## c(b0,b1,b2,…,bn) = b0+b1*x+b2*x^2+b3*x^3 + … + bn*x^n ## [alfa0, alfa1]: interval ## Output: ## vector with the local maxima and the local minima ## Example: ## maxminInterval(c(6,4,4,5,5),-4,4) ## maxminInterval(c(10,-60,55,-20,2.5),0,4) ## plot(VectPolynomial(c(6,4,4,5,5)),xlim=c(-4,4)) ## plot(VectPolynomial(c(10,-60,55,-20,2.5)),xlim=c(0,4)) maxminInterval <- function (vPol, alfa0, alfa1) { zerosPol (polynomialDerivative(vPol),alfa0,alfa1) } ## Function: ## max P(x) for x \in [alfa0, alfa1] ## Input: ## vPol: Polynomial (a vector of coefficients) ## c(b0,b1,b2,…,bn) = b0+b1*x+b2*x^2+b3*x^3 + … + bn*x^n ## [alfa0, alfa1]: interval ## Output: ## the global maximum in the given interval ## Example: ## maximumMaximorum(c(6,4,4,5,5),0,4) ## maximumMaximorum(c(10,-60,55,-20,2.5),0,4) maximumMaximorum <- function (vPol, alfa0, alfa1) { G <- maxminInterval(vPol, alfa0, alfa1) # print(G) max(sapply (G, VectPolynomial(vPol))) } ## Function: ## min P(x) for x \in [alfa0, alfa1] ## Input: ## vPol: Polynomial (a vector of coefficients) ## c(b0,b1,b2,…,bn) = b0+b1*x+b2*x^2+b3*x^3 + … + bn*x^n ## [alfa0, alfa1]: interval ## Output: ## the global minimum in the given interval ## Example: ## minimumMinimorum(c(6,4,4,5,5),0,4) ## minimumMinimorum(c(10,-60,55,-20,2.5),0,4) ## minimumMinimorum(c(25,-60,55,-20,2.5),0,4) ## plot(VectPolynomial(c(6,4,4,5,5)),xlim=c(0,4)) ## plot(VectPolynomial(c(10,-60,55,-20,2.5)),xlim=c(0,4)) minimumMinimorum <- function (vPol, alfa0, alfa1) { G <- maxminInterval(vPol, alfa0, alfa1) print(G) min(sapply (G, VectPolynomial(vPol))) } ## Function: ## This function defines a fuzzy measure (a distorted Lebesgue), ## applying a distortion to a set of intervals ## m(\sum (mx_i - mn_i) ) ## Input ## [intervalMin,intervalMax]: interval under consideration ## m: distortion ## Output: ## a function that receives a list of pairs of intervals and computes ## the distorted Lebesgue measure of this set ## m( [I(0)_n,I(0)_x] \cup [I(1)_n,I(1)_x] \cup ... \cup [I(n)_n,I(n)_x] ) ## m( (I(0)_x-I(0)_n) + (I(1)_x-I(1)_n) + ... (I(n)_x-I(n)_n) ) ## Example: ## ffaaux <- fmeasure(0,4,function(x) (x^2)) ## ffaaux(rbind(c(0,2),c(2,4))) ## ffaaux(rbind(c(0,1.2),c(2,4))) fmeasure <- function (intervalMin, intervalMax,m) { function (setIntervals) { vret <- 0 if (is.vector(setIntervals)) { vret <- m((setIntervals[2] - setIntervals[1])/ (intervalMax-intervalMin))} else { if (dim(setIntervals)[1]==0) { vret <- 0 } else { vret <- m(sum(apply(setIntervals,1, function(x){x[2]-x[1]}))/ (intervalMax-intervalMin)) } } return (vret) } } ## Function: ## For a given polynomial, and a value beta, this function computes mu(A_beta) ## Function that defines \mu_{\beta} = \mu(A_{\beta}) ## Input: ## vPol: The polynomial ## [alfa0,alfa1]: the interval under consideration ## fm: a fuzzy measure, it should evaluate a list of intervals and return the measure for this set ## Output: ## The fuzzy measure mu(A_beta) ## Example: ## fmdx2 <- fmeasure(0,4,function(x) (x^2)) ## muBeta<- muBetaPoly (c(25,-60,55,-20,2.5),0,4, fmdx2) ## sapply((1:250)/10, muBeta) ## muBeta(10) ## 0.03012099 ## plot(VectPolynomial(c(6,4,4,5,5)),xlim=c(-4,4)) ## plot(VectPolynomial(c(10,-60,55,-20,2.5)),xlim=c(0,4)) ## plot(VectPolynomial(c(25,-60,55,-20,2.5)),xlim=c(0,4)) muBetaPoly <- function (vPol, alfa0, alfa1,fm) { Vectorize(function (beta) { vPolN <- vPol vPolN[1] <- vPolN[1]-beta A <- zerosPol (vPolN, alfa0, alfa1) Ainterval <- cbind(A[1:length(A)-1], A[2:length(A)]) fVPol <- VectPolynomial(vPolN) selection <- apply(Ainterval,1, function(x) {(fVPol((x[1]+x[2])/2))> 0}) Ab <- Ainterval[selection,] mub <- fm (Ab) return(mub) },c("beta")) } ## Function: ## Integrates a polynomial in a given interval with respect to a fuzzy measure ## Input: ## vPol: Polynomial (a vector of coefficients) ## c(b0,b1,b2,…,bn) = b0+b1*x+b2*x^2+b3*x^3 + … + bn*x^n ## [alfa0,alfa1]: Interval under consideration ## fmdx2: a fuzzy measure, which given a set of intervals returns their measure ## Output: ## the integral of the polynomial in [alfa0,alfa1] w.r.t. fmdx2 ## $value: integral ## $abs.error: error of the integral ## Example: ## ffaaux <- fmeasure(0,4,function(x) (x^2)) ## ciPositivePolynomial(c(25,-60,55,-20,2.5),0,4, ffaaux) ciPositivePolynomial <- function (vPol, alfa0, alfa1, fmdx2) { muBeta<- muBetaPoly (vPol,alfa0,alfa1, fmdx2) mnMn <- minimumMinimorum (vPol, alfa0, alfa1) if (mnMn<0) { print ("This function is only for positive polynomials. Case not considered"); return() } I <- integrate(muBeta, lower=mnMn, upper=maximumMaximorum (vPol, alfa0, alfa1)) nI <- fmdx2(rbind(c(alfa0,alfa1)))* mnMn res<-list(value = I$value + nI, abs.error = I$abs.error) return (res) } ### Example 4 @ (Torra, Narukawa, 2016) ### http://dx.doi.org/10.1016/j.inffus.2016.02.007 ### This proceeds from Example 2 of ANOR paper ### ### Integral of the polynomials 4x-4x^2 and 2x - x^2 ### (both polynomials result into the same integral for distorted Lebesgue measures) ### w.r.t. distorted Lebesgue with m(x)=x^n fmEx4l1 <- fmeasure(0,1,function(x) (x)) fmEx4l2 <- fmeasure(0,1,function(x) (x^2)) fmEx4l3 <- fmeasure(0,1,function(x) (x^3)) fmEx4l4 <- fmeasure(0,1,function(x) (x^4)) ciPositivePolynomial(c(0,4,-4),0,1, fmEx4l1) ciPositivePolynomial(c(0,4,-4),0,1, fmEx4l2) ciPositivePolynomial(c(0,4,-4),0,1, fmEx4l3) ciPositivePolynomial(c(0,4,-4),0,1, fmEx4l4) ciPositivePolynomial(c(0,2,-1),0,1, fmEx4l1) ciPositivePolynomial(c(0,2,-1),0,1, fmEx4l2) ciPositivePolynomial(c(0,2,-1),0,1, fmEx4l3) ciPositivePolynomial(c(0,2,-1),0,1, fmEx4l4) ### To display the functions ### polynA <- VectPolynomial(c(0,4,-4)) ### plot(polynA,xlim=c(0,1)) ### polynB <- VectPolynomial(c(0,2,-1)) ### plot(polynB,xlim=c(0,1)) ### 4x-4x^2 has a single peak in [0,1]. We can use the single peaked function integral anorEx2fun <- function (x) (4*x - 4*x*x) anorEx2lam0 <- function (a) (a)^(0) anorEx2lam1 <- function (a) (a)^(1) anorEx2lam2 <- function (a) (a)^(2) anorEx2lam3 <- function (a) (a)^(3) anorEx2lam4 <- function (a) (a)^(4) anorEx2DLeb0 <- function (a,b) (anorEx2lam0(b-a)) anorEx2DLeb1 <- function (a,b) (anorEx2lam1(b-a)) anorEx2DLeb2 <- function (a,b) (anorEx2lam2(b-a)) anorEx2DLeb3 <- function (a,b) (anorEx2lam3(b-a)) anorEx2DLeb4 <- function (a,b) (anorEx2lam4(b-a)) ciSinglePeakDProbabilityWithInverses(anorEx2fun,NULL,NULL,anorEx2DLeb0,0.5,0,1) ## $value: [1] 1 $abs.error: [1] 1.110223e-14 ciSinglePeakDProbabilityWithInverses(anorEx2fun,NULL,NULL,anorEx2DLeb1,0.5,0,1) ## $value: [1] 0.6666651 $abs.error: [1] 8.686017e-05 ciSinglePeakDProbabilityWithInverses(anorEx2fun,NULL,NULL,anorEx2DLeb2,0.5,0,1) ## $value: [1] 0.5000031 $abs.error: [1] 0.0001035899 ciSinglePeakDProbabilityWithInverses(anorEx2fun,NULL,NULL,anorEx2DLeb3,0.5,0,1) ## $value: [1] 0.4000005 $abs.error: [1] 8.360559e-05 ciSinglePeakDProbabilityWithInverses(anorEx2fun,NULL,NULL,anorEx2DLeb4,0.5,0,1) ## $value: [1] 0.3333344 $abs.error: [1] 0.0001189902 ### Comparison with the analytical expression in the same example anorEx2Res <- function (n) { (2/(n+1))- (2/((n+2)*(n+1)))} ### anorEx2Res(0) ### anorEx2Res(1) ### anorEx2Res(2) ### anorEx2Res(3) ### anorEx2Res(4) ### Other examples ### Example A ## vpol <- VectPolynomial(c(25,-60,55,-20,2.5)) ## plot(vpol,xlim=c(0.0,4)) ## plot(vpol,xlim=c(0.5,3.5)) ## fmdx2 <- fmeasure(0,4,function(x) (x^2)) ## ciPositivePolynomial(c(25,-60,55,-20,2.5),0,4, fmdx2) ### $value -> [1] 3.916354 ; $abs.error -> [1] 0.0001692622 ## ### Example B ## fmdx201 <- fmeasure(0,1,function(x) (x*x)) ## vpol2 <- VectPolynomial(c(0,0,0,1)) ## plot(vpol2,xlim=c(0,1)) ## ciPositivePolynomial(c(0,0,0,1),0,1, fmdx201) ### $value -> [1] 0.1 ; $abs.error -> [1] 2.638411e-07 ## ### Comparison between the integral of a polynomial and as an increasing function ### Example C_1 ## polFunction1 <- c(0,0,0,1) # function1 <- function (x) { x*x*x } ## fmdx <- fmeasure(0,1,function(x) (x)) ## ciPositivePolynomial(c(0,0,0,1),0,1, fmdx) ## ciIncreasingDProbabilityWithInverse(function1,inverse1, dpmuId,0,1) ## ciIncreasingDProbabilityWithInverse(function1,NULL, dpmuId,0,1) ### Example C_2 ## polX <- c(0,1,2,3) ## vpolX <- VectPolynomial(polX) ## plot(vpolX,xlim=c(0,1)) ## fmdx <- fmeasure(0,1,function(x) (x)) ## ciPositivePolynomial(polX,0,1, fmdx) ## ciIncreasingDProbabilityWithInverse(vpolX,NULL, dpmuId,0,1) ### Example C_3 ## polX <- c(1,0.23,2,3,4,5) ## vpolX <- VectPolynomial(polX) ## plot(vpolX,xlim=c(0,1)) ## fmdx <- fmeasure(0,1,function(x) (x)) ## ciPositivePolynomial(polX,0,1, fmdx) ## ciIncreasingDProbabilityWithInverse(vpolX,NULL, dpmuId,0,1)