业务临时需要,需要确定多因素影响下的综合权重值,现使用层次分析法和拉格朗日多项式插值算法做简易值计算
1、建立层次分析结构模型,分析影响综合指标的各个因素,分层级,上层受下层影响,而同层各因素之间基本上相对独立(本实例只有两层)
2、连接oracle数据库,读取各个因素实际指标值
library(DBI)
library(ROracle)
drv=dbDriver('Oracle')
conn=dbConnect(drv,'AQTS_ZHZX','AQTS_ZHZX','192.160.0.0:1521/iims')
rs=dbSendQuery(conn,"with weather as
(
select * from (
select * from weather_fact_aqts w
where w.cjrq=(select max(cjrq) from weather_fact_aqts)
order by w.cjsd desc) where rownum=1
)
select v.yfbm,v.dlbm,v.ldbm,weather.cjrq,weather.cjsd,v.wfshb,a.sgshb,
weather.sbbh,weather.wd,weather.sd,weather.fl,weather.njd
,weather.fog_rank,weather.wind_rank from VIO_FACT_AQTS v
join ACD_FACT_AQTS a on v.yfbm=a.sgyfbm and v.dlbm=a.lhbm and v.ldbm=a.glsbm
join weather on 1=1")
data=fetch(rs)
data=as.matrix(data)
3、影响综合因素的指标分两级:
#FIRST,A1:weather,A2:Flux,A3:acd,A4:vio
#SECOND1,B1:fog,B2:rain,B3:wind,B4:snow
#SECOND2,B2:FLUX1,FLUX2
4、构造第一层成对比较矩阵(采用成对比较法和1~9尺度)
#FIRST,A1:weather,A2:Flux,A3:acd,A4:vio
A=matrix(0,4,4)
diag(A)=1
A[1,2]=1/2;A[2,1]=2
A[1,3]=1/4;A[3,1]=4
A[1,4]=1/4;A[4,1]=4
A[2,4]=1/2;A[4,2]=2
A[2,3]=1/2;A[3,2]=2
A[3,4]=A[4,3]=1
5、构造第二层成对比较矩阵(构造针对上一层成对比较矩阵)
#SECOND,B1:fog,B2:rain,B3:wind,B4:snow
B=matrix(0,4,4)
diag(B)=1
B[1,2]=2;B[2,1]=1/2
B[1,3]=6;B[3,1]=1/6
B[1,4]=2;B[4,1]=1/2
B[2,3]=3;B[3,2]=1/3
B[2,4]=1;B[4,2]=1
B[3,4]=1/3;B[4,3]=3
#SECOND2,B2:FLUX1,FLUX2
C=matrix(1,2,2)
C[1,2]=3
C[2,1]=1/3
6、计算权向量并一致性校验
求权向量以及特征向量,进行一致性检验(因时间问题,这块没有加校验,后续修改)
7、求各层级因素权重,这里使用规范列平均法
Ac=colSums(A)
for (i in 1:nrow(A)) {
A[,i]=A[,i]/Ac[i]
}
weigA=rowMeans(A)
#B
Bc=colSums(B)
for (i in 1:nrow(B)) {
B[,i]=B[,i]/Bc[i]
}
weigB=rowMeans(B)
#C 1:flux,2:vehicle
Cc=colSums(C)
for (i in 1:nrow(C)) {
C[,i]=C[,i]/Cc[i]
}
weigC=rowMeans(C)
########row 1:fog,2:rain,3:wind,4:snow,5:flux,6:vehicle,7:acd,8:vio
########new row 1:fog,2:rain,3:wind,4:snow,5:even,6:acd,7:vio
weit=matrix(0,8,4)
weit[1:4,1]=weigB
weit[5:6,2]=weigC
weit[7,3]=1
weit[8,4]=1
######
w=as.matrix(weigA)
p=weit
we=p%*%w
wenew=c(we[1:4],sum(we[5:6]),we[7],we[8])
8、构造各个因素的函数
#get the function of accident
acdx<-c(0,0.1,0.2,0.3,0.4)
acdy<-c(5,4,3,2,1)
acdf <- LagrangePolynomial(acdx,acdy)
acdv <- function(x,f){
x=as.numeric(x)
ranka=f(x)
ranka[ranka>5]=5
ranka[ranka<0]=0
return (ranka)
}
#v=c(-0.1,0.2)
#vv=acdv(v,acdf)
#get the function of violation
viox<-c(0,0.1,0.2,0.3,0.4)
vioy<-c(5,4,3,2,1)
viof <- LagrangePolynomial(viox,vioy)
#get the function of flux flux
fluxfx<-c(0,0.15,0.3,0.45,0.6)
fluxfy<-c(5,4,3,2,1)
fluxff <- LagrangePolynomial(fluxfx,fluxfy)
#get the function of flux vehicel
fluxvx<-c(0,0.2,0.35,0.5,0.65)
fluxvy<-c(5,4,3,2,1)
fluxvf <- LagrangePolynomial(fluxvx,fluxvy)
#get the function of fog rain snow
fogx<-c(0,50,100,200,500)
fogy<-c(1,2,3,4,5)
fogf <- LagrangePolynomial(fogx,fogy)
#get the function of wind
windx<-c(0,13.9,20.7,28.5,32.7)
windy<-c(5,4,3,2,1)
windf <- LagrangePolynomial(windx,windy)
9、对应各因素指标值与权重相乘求取综合指标级别,并写入数据库
rankc=cbind(data[,c(2:5)],as.matrix(acdv(data[,12],fogf))*wenew[1],as.matrix(acdv(data[,12],fogf))*wenew[2],
as.matrix(acdv(data[,11],fogf))*wenew[3],as.matrix(acdv(data[,12],fogf))*wenew[4],as.matrix(rep(3,nrow(data)))*wenew[5],
as.matrix(acdv(data[,7],acdf))*wenew[6],as.matrix(acdv(data[,6],viof))*wenew[7])
#comprehensive rank
alv=as.numeric(rankc[,5:11])
dim(alv)=c(nrow(rankc),7)
RANKOVER<-cbind(rankc[,1:4],ceiling(rowSums(alv)))
colnames(RANKOVER)<-c('DLBM','LDBM','CJRQ','CJSD','RANK')
RANKOVER=as.data.frame(RANKOVER)
#dbRemoveTable(conn, 'ZHZX_FACT_AQTS')
dbWriteTable(conn,'ZHZX_FACT_AQTS',RANKOVER, row.names = F, append = TRUE)
附:拉格朗日多项式插值算法
#LagrangePolynomial
LagrangePolynomial <- function(x,y) {
len = length(x)
if(len != length(y))
stop("length not equal!")
if(len < 2)
stop("dim size must more than 1")
#pretreat data abd alloc memery
xx <- paste("(","a -",x,")")
m <- c(rep(0,len))
#combin express
for(i in 1:len) {
td <- 1
tm <- "1"
for(j in 1:len) {
if(i != j) {
td <- td*(x[i] - x[j])
tm <- paste(tm,"*",xx[j])
}
}
tm <- paste(tm,"/",td)
m[i]<-tm #m[i] <- parse(text=tm)
}
#combin the exrpession
m <- paste(m,"*",y)
r <- paste(m,collapse="+")
#combin the function
fbody <- paste("{ return(",r,")}")
f <- function(a) {}
#fill the function's body
body(f) <- parse(text=fbody)
return(f)
}
(记录日常点滴)对算法以及R使用还不是很熟悉,后续完善!如有错误,希望指点>>>OK