求区间可能度矩阵的算法(Liu 2009)

uncertain
作者

zsc

发布日期

2019年5月4日

求区间可能度矩阵的算法

0. 参考论文:

徐泽水： A consistency improving method in the analytic hierarchy process 1999年

刘芳：Acceptable consistency analysis of interval reciprocal comparison matrices 2009年

1. 主要思路：

  1. 把区间乘性互反矩阵U拆成两个正互反判断矩阵B和D,B，D都是正的互反判断矩阵。其中B的下三角元素大于D矩阵的下三角元素，B的上三角小于D的上三角元素 ，简称B的下三角大，上三角小
  2. 一致性检验，若拆分后的B，D矩阵一致性不满足条件(即\(CR<= 0.1\)) ，则用徐泽水(1999年)的文章方法进行调整，直到满足一致性条件为准(\(CR <=0.1\)).
  3. 然后分别计算矩阵\(B，D\)的权重向量\(w(B),w(D)\) ,注意这里的权重没有归一化处理.
  4. 通过公式\(w_i = [min(w_i(B),w_i(D)),max(w_i(B),w_i(D))]\),把两个权重向量组合成一个区间向量。
  5. 通过区间向量\(w\)计算出区间向量的可能度矩阵\(P\)。

2.主要函数构建：

  1. consistency(A): 求正互反判断矩阵的一致性指标，返回一个list

  2. em_get_w(A) : 特征值求权重 —— 没有归一化权重

  3. gm_get_w(A)： 几何平均求权重 — — 没有归一化权重

  4. get_w(B,D): 分别获取B，D的权重（可以指定几何平均或者特征值求权重），然后组成区间权重向量（即小的在前，大的在后），这里返回的是一个矩阵，把每一个区间数看做矩阵的一行。

  5. fenjie(U ): 把区间矩阵U分解成正互反判断矩阵B和D

  6. adjust_w(A,lambda) : 利用论文的方法进行调整，返回调整后符合一致性条件的一致性矩阵。

  7. degree_probability(a,b) 函数计算两个区间数的可能度

  8. probability_matrix(w) 给一个n*2 的区间数，求其可能度矩阵

Show the code

rm(list = ls())
# 0。 一致性指标的求解
consistency = function(A){
  lambda = Re(eigen(A)$values[1]) # 矩阵A的最大特征值
  n  = nrow(A)
  RI = c(0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45);
  CI = (lambda-n) / (n-1);
  CR = CI / RI[n];
  eig_w = eigen(A)$vectors[,1] / sum( eigen(A)$vectors[,1]);
  return(list("eig_value"=lambda,"CI"=CI,"RI"=RI[n],"CR"=CR,'eig_w'=Re( eig_w )))
}



# 1. 特征值求权重
em_get_w = function(A){
  n = nrow(A)
  stopifnot(nrow(A) ==  ncol(A))
  
  lambda = Re(eigen(A)$values[1]) # 矩阵A的最大特征值
  n = nrow(A)
  RI = c(0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45);
  CI=(lambda-n)/(n-1);
  CR=CI/RI[n];
  eig_origin = eigen(A)$vectors[,1]
  eig_w = eigen(A)$vectors[,1] # /sum(eigen(A)$vectors[,1]);
  eig_w = Re(eig_w)
  return(eig_w)
}
# 2. 几何平均求权重 -- 
gm_get_w =function(A){
  n = nrow(A)
  stopifnot(nrow(A) ==  ncol(A))
  temp = apply(A, 1, function(x) prod(x)^(1/n) )
  w = temp # /sum(temp) 
  return(w)
}
# 3. 合并权重
get_w = function(B,D,method = c('gm_get_w','em_get_w') ){
  n = nrow(B)
  stopifnot(n == nrow(D))
  method <- match.arg(method)
  f = get(method)
  w_B = f(B)
  w_D = f(D)
  w_L = rep(0,length(w_B))
  w_U = rep(0,length(w_B))
  for(i in 1:length(w_B)){
    w_L[i] = min(w_B[i],w_D[i])
    w_U[i] = max(w_B[i],w_D[i])
  }
  w = matrix(c(w_L,w_U), ncol= 2 )
  return(w)
}



# 4. 通过U进行分解，分解出B，D矩阵，
fenjie = function(U){
  n = nrow(U)
  stopifnot(ncol(U) == 2*n)
  
  B = matrix(0,nrow = n,ncol = n)
  D = matrix(0,nrow = n,ncol = n)
  
  for(i in 1:n){
    for(j in 1:n){
      if(i<j){
        B[i,j] = U[i,j*2]
        D[i,j] = U[i,2*j-1]
      }else if(i>j){
        B[i,j] = U[i,2*j-1]
        D[i,j] = U[i,j*2]
      }else{
        B[i,j] = U[i,2*j]
        D[i,j] =U[i,2*j]
      }
    }
  }
  return(list('B'=B,'D'=D))
}


# 5. 徐泽水(1999年)的文章方法进行调整，直到满足一致性条件为准(CR <=0.1).
adjust_w <- function(A, lambda) {
  k <- 0
  n = nrow(A)
  m = ncol(A)
  stopifnot(n == m)
  temp_CR <- consistency(A)$CR
  temp_w <- consistency(A)$eig_w
  while (temp_CR >= 0.1 && k < 1000) {
    for (i in 1:n) {
      for (j in 1:n) {
          A[i, j] <- (A[i, j]^lambda) * (temp_w[i] / temp_w[j])^(1 - lambda)
      }
    }
    temp_CR <- consistency(A)$CR
    temp_w <- consistency(A)$eig_w
    k <- k + 1
  }
  return(A)
}



# 6. degree_probability 函数计算两个区间数的可能度
degree_probability <- function(a, b) {
  # 输入的a,b代表一个区间，即是一个二维向量，且小的在前面，大的元素在后面
  stopifnot(length(a) == length(b), length(a) == 2, a[1] <= a[2], b[1] <= b[2])
  temp <- 0
  if (a[1] == a[2] && b[1] == b[2]) {
    if (a[1] > b[1]) {
      temp <- 1
    } else if (a[1] == b[1]) {
      temp <- 0.5
    } else {
      temp <- 0
    }
  } else if (a[1] == a[2] && b[1] != b[2]) {
    if (a[1] > b[2]) {
      temp <- 1
    } else if (b[1] <= a[1] & a[1] <= b[2]) {
      temp <- (a[1] - b[1]) / (b[2] - b[1])
    } else {
      temp <- 0
    }
  } else if (a[1] != a[2] && b[1] == b[2]) {
    if (a[1] > b[1]) {
      temp <- 1
    } else if (a[1] <= b[1] & b[1] <= a[2]) {
      temp <- (a[2] - b[1]) / (a[2] - a[1])
    } else {
      temp <- 0
    }
  } else if (a[1] != a[2] && b[1] != b[2]) {
    if (a[1] < a[2] && a[2] <= b[1] && b[1] < b[2]) {
      temp <- 0
    } else if (a[1] <= b[1] && b[1] < a[2] && a[2] <= b[2]) {
      s_t <- (a[2] - b[1]) * (a[2] - b[1]) * 0.5
      s <- (b[2] - b[1]) * (a[2] - a[1])
      temp <- s_t / s
    } else if (a[1] <= b[1] && b[1] <= b[2] && b[2] <= a[2]) {
      s_t <- ((a[2] - b[2]) + (a[2] - b[1])) * (b[2] - b[1]) * 0.5
      s <- (b[2] - b[1]) * (a[2] - a[1])
      temp <- s_t / s
    } else if (b[1] < a[1] && a[1] < a[2] && a[2] < b[2]) { 
      # 可写等号b[1] <= a[1] && a[1]<a[2] &&a[2]<=b[2]
      s_t <- ((a[1] - b[1]) + (a[2] - b[1])) * (a[2] - a[1]) * 0.5
      s <- (b[2] - b[1]) * (a[2] - a[1])
      temp <- s_t / s
    } else if (b[1] < a[1] && a[1] < b[2] && b[2] < a[2]) {
      s_tt <- (b[2] - a[1]) * (b[2] - a[1]) * 0.5
      s <- (b[2] - b[1]) * (a[2] - a[1])
      s_t <- s - s_tt
      temp <- s_t / s
    } else if (b[1] < b[2] && b[2] <= a[1] && a[1] < a[2]) {
      temp <- 1
    } else {
      stop("运行出错")
    }
  } else {
    stop("运行出错,请检查")
  }
  return(temp)
}

# 7. probability_matrix 给一个n*2 的区间数，求其可能度矩阵
probability_matrix <- function(Z) {
  # probability_matrix函数输入一个n*2的矩阵，每一行代表输出各个方案的综合属性值得区间数
  # 此函数输出各方案两两比较的可能度矩阵。
  # degree_probability函数求两个区间数的可能度，
  # a，b代表输入的区间数，输入这两个
  P <- matrix(0, ncol = nrow(Z), nrow = nrow(Z))
  for (i in 1:nrow(Z)) {
    for (j in 1:nrow(Z)) {
      P[i, j] <- degree_probability(Z[i, ], Z[j, ])
    }
  }
  return(P)
}

总结: 只需给出一个区间判断矩阵，返回最终的权重(前提是B和D要满足一致性条件\(CR<=0.1\))

3. 测试

3.1 例1： B 和D都满足一致性指标

Show the code

U = matrix(c(1,1,2,5,2,4,1,3,
             1/5,1/2,1,1,1,3,1,2,
             1/4,1/2,1/3,1,1,1,1/2,1,
             1/3,1,1/2,1,1,2,1,1),nrow = 4,byrow = T)
fenjie(U) 
#> $B
#>           [,1]      [,2] [,3] [,4]
#> [1,] 1.0000000 5.0000000    4    3
#> [2,] 0.2000000 1.0000000    3    2
#> [3,] 0.2500000 0.3333333    1    1
#> [4,] 0.3333333 0.5000000    1    1
#> 
#> $D
#>      [,1] [,2] [,3] [,4]
#> [1,]  1.0    2    2  1.0
#> [2,]  0.5    1    1  1.0
#> [3,]  0.5    1    1  0.5
#> [4,]  1.0    1    2  1.0

## 一致性检验
library(purrr)
fenjie(U) %>% map(function(x)consistency(x)$CR) # 求解矩阵B，D的一致性指标
#> $B
#> [1] 0.08209909
#> 
#> $D
#> [1] 0.02246186

# 直接求出区间权重，以及根据权重求出区间可能度矩阵
( w = get_w(B= fenjie(U)$B ,D = fenjie(U)$D) ) #每一行对应第i个方案的区间权重
#>           [,1]      [,2]
#> [1,] 1.4142136 2.7831577
#> [2,] 0.8408964 1.0466351
#> [3,] 0.5372850 0.7071068
#> [4,] 0.6389431 1.1892071

( P = probability_matrix(w) ) # 可能度矩阵P
#>      [,1]      [,2]      [,3]       [,4]
#> [1,]  0.5 1.0000000 1.0000000 1.00000000
#> [2,]  0.0 0.5000000 1.0000000 0.55395713
#> [3,]  0.0 0.0000000 0.5000000 0.02486059
#> [4,]  0.0 0.4460429 0.9751394 0.50000000

3.2 例2： B 和D不满足一致性指标

Show the code

U  = matrix(c(1,1,1,2,1,2,2,3,
              1/2,1,1,1,3,5,4,5,
              1/2,1,1/5,1/3,1,1,6,8,
              1/3,1/2,1/5,1/4,1/8,1/6,1,1),nrow = 4,byrow = T)
fenjie(U)
#> $B
#>           [,1] [,2]  [,3] [,4]
#> [1,] 1.0000000  2.0 2.000    3
#> [2,] 0.5000000  1.0 5.000    5
#> [3,] 0.5000000  0.2 1.000    8
#> [4,] 0.3333333  0.2 0.125    1
#> 
#> $D
#>      [,1]      [,2]      [,3] [,4]
#> [1,]  1.0 1.0000000 1.0000000    2
#> [2,]  1.0 1.0000000 3.0000000    4
#> [3,]  1.0 0.3333333 1.0000000    6
#> [4,]  0.5 0.2500000 0.1666667    1
fenjie(U) %>% map(function(x)consistency(x)$CR) # 可以发现B和D的一致性不满足条件
#> $B
#> [1] 0.2752221
#> 
#> $D
#> [1] 0.1232627

## 方法一： 进行调整
adjust_w(fenjie(U)$B,lambda = 0.6) %>% round(4)
#>        [,1]   [,2]   [,3]   [,4]
#> [1,] 1.0000 1.4696 1.9339 3.9751
#> [2,] 0.6805 1.0000 3.4564 5.5703
#> [3,] 0.5171 0.2893 1.0000 5.6119
#> [4,] 0.2516 0.1795 0.1782 1.0000
adjust_w(fenjie(U)$D,lambda = 0.88) %>% round(4)
#>        [,1]   [,2]   [,3]   [,4]
#> [1,] 1.0000 0.9428 0.9943 2.1061
#> [2,] 1.0607 1.0000 2.7732 4.1112
#> [3,] 1.0057 0.3606 1.0000 5.5695
#> [4,] 0.4748 0.2432 0.1795 1.0000

## 方法二:进行调整 --  多参数映射
temp = fenjie(U) %>% map2(.,list(0.6,0.88),function(x,y)round( adjust_w(x,y),4) )
temp
#> $B
#>        [,1]   [,2]   [,3]   [,4]
#> [1,] 1.0000 1.4696 1.9339 3.9751
#> [2,] 0.6805 1.0000 3.4564 5.5703
#> [3,] 0.5171 0.2893 1.0000 5.6119
#> [4,] 0.2516 0.1795 0.1782 1.0000
#> 
#> $D
#>        [,1]   [,2]   [,3]   [,4]
#> [1,] 1.0000 0.9428 0.9943 2.1061
#> [2,] 1.0607 1.0000 2.7732 4.1112
#> [3,] 1.0057 0.3606 1.0000 5.5695
#> [4,] 0.4748 0.2432 0.1795 1.0000
temp %>% map(function(x)consistency(x)$CR)#检验权重情况
#> $B
#> [1] 0.09525181
#> 
#> $D
#> [1] 0.09474896


# 求出调整后的B和D的权重
( w = get_w(B= temp$B ,D = temp$D) ) #每一行对应第i个方案的区间权重
#>           [,1]      [,2]
#> [1,] 1.1853702 1.8333497
#> [2,] 1.8648143 1.9025351
#> [3,] 0.9572122 1.1921409
#> [4,] 0.2995165 0.3794326
( P = probability_matrix(w) ) # 可能度矩阵P
#>              [,1] [,2]      [,3] [,4]
#> [1,] 0.5000000000  0.0 0.9998494  1.0
#> [2,] 1.0000000000  0.5 1.0000000  1.0
#> [3,] 0.0001505719  0.0 0.5000000  1.0
#> [4,] 0.0000000000  0.0 0.0000000  0.5

4. 总结:

无论什么样的区间矩阵(一致性是否满足)，都可以用步骤

# list(0.6,0.88)中的0.6和0.88 为调整adjust_w()函数参数中的lambda值对应
temp = fenjie(U) %>% map2(., list(0.6,0.88),function(x,y)round( adjust_w(x,y),4) )
w = get_w(B= temp$B ,D = temp$D)#每一行对应第i个方案的区间权重
w

P = probability_matrix(w)
P