library(RODPS) library(rattle) library(stringr) library(sqldf) library(Hmisc) library(reshape) load("d:/temp/GeShuiDaShuJu2.RData") #rodps.init("D:/R/R-3.2.4revised/library/RODPS/odps_config.ini") #原始指标数据 risk<-rodps.read.table("udaftest3") #平均工资上限 pjgz_sh<-as.character(round(mean(risk$pjgz)+20*sd(risk$pjgz),2)) #计算下降总量,排除异常点 sql<-paste("select djxh,(case when risk>=6 then risk-5 else risk end) as risk, perc,pjrs,amount,round(pjrs*amount,2) as zl from risk where pjgz<=",pjgz_sh,sep="") riskindexs<-sqldf(sql) #下降总量取自然对数,用经过自然对数变换的总量计算风险指数更合理 riskindexs$lnzl<-round(log(riskindexs$zl),2) #e^0=1,log(1)=0,有10条 sqldf("select count(*) as c from riskindexs where lnzl=0") #log(0)= -Inf sqldf("select count(*) as c from riskindexs where lnzl<0") #特殊处理, -Inf的全部置0 riskindexs$lnzl[which(riskindexs$lnzl<0)]<-0 #各指标变量作0~1标准化处理,去除量纲影响 riskindexs$risk<-rescaler(riskindexs$risk, "range") riskindexs$perc<-rescaler(riskindexs$perc, "range") riskindexs$pjrs<-rescaler(riskindexs$pjrs, "range") riskindexs$amount<-rescaler(riskindexs$amount, "range") riskindexs$zl<-rescaler(riskindexs$zl, "range") riskindexs$lnzl<-rescaler(riskindexs$lnzl, "range") #按2个不同标准计算风险指数,区别是总量是否经过对数变换,这是未作对数变换的 riskindexs$riskindex<-riskindexs$zl*0.7+riskindexs$risk*0.2+riskindexs$perc*0.1 #0~1标准化 riskindexs$riskindex<-rescaler(riskindexs$riskindex, "range") #按2个不同标准计算风险指数,区别是总量是否经过对数变换,这是经过对数变换的 riskindexs$riskindex2<-riskindexs$lnzl*0.7+riskindexs$risk*0.2+riskindexs$perc*0.1 #0~1标准化 riskindexs$riskindex2<-rescaler(riskindexs$riskindex2, "range") #二分类目标变量,打上标签,是否有风险 riskindexs$tag<-ifelse(riskindexs$risk>0,"1","0") riskindexs$tag<-as.factor(riskindexs$tag) names(riskindexs)<-c("djxh","risk_01","perc_01","pjrs_01","amount_01","zl_01","lnzl_01","riskindex","riskindex2","tag") riskindexs<-merge(risk,riskindexs,by="djxh",all.Y=TRUE) riskindexs$hy_dm<-substr(riskindexs$hy_dm,1,2) riskindexs<-riskindexs[order(-riskindexs$riskindex2),] #只有1条,不便于划分训练集、验证集、测试集,予以排除。 sqldf("select count(*) from riskindexs where kzztdjlx_dm='1100'") #sqldf("select * from kzztdjlx where kzztdjlx_dm='1100'") riskindexs<-riskindexs[which(riskindexs$kzztdjlx_dm!="1100"),] #解决rattle部分汉字显示为乱码的问题,设定R语言的提示信息为英文。 #Sys.setenv(LANGUAGE="zh") Sys.setenv(LANGUAGE="en") rattle(dataset="riskindexs") #对模式识别算法结果总体情况观察 #所有有风险的记录 sqldf("select count(*) from riskindexs where riskindex>0") tmp<-sqldf("select * from riskindexs where riskindex>0") #画下降总量概率密度图看看,有指数分布的形状 p1<-ggplot(data=riskindexs,aes(x=amount_01))+ geom_density(colour="green")+ stat_density(position='identity',fill='darkblue')+ ylab("概率密度")+xlab("下降总量")+ ggtitle("下降总量=(平均人数*下降额度)")+ theme(text=element_text(size=20)) print(p1) #画风险指数1概率密度图看看,有正态分布的形状 p1<-ggplot(data=tmp,aes(x=riskindex))+ geom_density(colour="green")+ stat_density(position='identity',fill='darkblue')+ ylab("概率密度")+xlab("风险指数")+ ggtitle("风险指数1 概率密度图\n=(平均人数*下降额度)*70%+年份*20%+下降幅度*10%")+ theme(text=element_text(size=20)) print(p1) sqldf("select count(*) from riskindexs where riskindex2>0") tmp<-sqldf("select * from riskindexs where riskindex2>0") #画风险指数2概率密度图看看,有正态分布的形状 p1<-ggplot(data=tmp,aes(x=riskindex2))+ geom_density(colour="green")+ stat_density(position='identity',fill='darkblue')+ ylab("概率密度")+xlab("风险指数")+ ggtitle("风险指数2 概率密度图\n=ln(平均人数*下降额度)*70%+年份*20%+下降幅度*10%")+ theme(text=element_text(size=20)) print(p1) #单样本KS测试,测试是否服从正态分布, p-value < 2.2e-16,<0.05,拒绝服从正态分布的原假设 #注意Rattle的K-S测试是比较两个样本是否来自同一个分布,用 library(fBasics)的ks2Test()函数,p值的判断不同。 ks.test(tmp$riskindex2,"pnorm",alternative ="two.sided") #画Q-Q图进一步确认,注意纵横坐标都是分位数,但横坐标分位数的50%分位坐标值为0,所以有正有负。 #非标准正态分布的斜率为样本标准差,截距为样本均值。 #在Q-Q图的首尾两段偏离直线,正态分布的Q-Q图是一条直线。但这个分布在中间一段接近正态分布,指数公式效果不错。 p1<-ggplot(data=tmp, aes(sample = riskindex2)) + stat_qq(colour="darkblue") + stat_qq_line(colour="red")+ ylab("样本累积概率密度")+xlab("理论累积概率密度")+ ggtitle("风险指数2 Q-Q图\n=ln(平均人数*下降额度)*70%+年份*20%+下降幅度*10%")+ theme(text=element_text(size=20)) print(p1) #这是基础作图工具的画法 #qqnorm(tmp$riskindex2) #qqline(tmp$riskindex2) #代缴单位统计,需要加上行业代码,经济类型代码及纳税人单位类型代码,以作微观分析。 #这个表是上传大集中金三个税代缴单位对照表到阿里云端后重新生成的。后面的微观分析都基于该表。 dwtj2<-rodps.table.read("dwrstj3") # 选择要查看的单位 test<-sqldf("select * from riskindexs where gzze>=10000000 and pjgz>=150000 and pjgz<=250000 and pjrs>=100 and risk>='14' order by riskindex2 desc") test<-head(test,6) # 选出这些单位历年的数据,画折线图看看 fxqy<-sqldf("select * from dwtj2 where djxh in(select distinct djxh from test) order by djxh,yy") fxqy$yy<-as.factor(fxqy$yy) fxqy$djxh<-substr(fxqy$djxh,11,18) #画折线图 key <- list(title = "登记序号", x = .01, y = .99,cex=1.5,lines = TRUE,points=TRUE,border=FALSE,columns=4) xyplot(pjgz/10000~yy,data=fxqy,type="b",groups=djxh, auto.key=key, #par.settings = parSettings, scales=list(cex=1.5), main=list("部分代缴单位历年平均年薪折线图",cex=2),xlab=list("年",cex=1.5),ylab =list("平均年薪(万)",cex=1.5)) # 读入企业所得税年报,只看看营业收入、营业成本、营业收入3列 library(data.table) selectedCols<-c(1,2,5,6,10,11,19) qysdsnb<-fread("D:/ThinkPad_X1_Carbon/temp/企业所得税年报.csv",header=TRUE,sep=",", select=selectedCols,stringsAsFactors = FALSE,showProgress = TRUE,data.table = FALSE) # 选出单位的年报 fxqynb<-qysdsnb[which(qysdsnb$DJXHDS %in% test$djxh | qysdsnb$DJXHGS %in% test$djxh ),] fxqynb$DJXHDS<-substr(fxqynb$DJXHDS,11,18) fxqynb$DJXHGS<-substr(fxqynb$DJXHGS,11,18) #读取随机森林模型False Negative误判输出,再抽取其历年平均工资明细,画出其变化图,分析具体情况 #可以看到,并非全部误判,随机森林算法找出了模式匹配算法没有发现的一些有风险的变化, #可以补充模式匹配算法受参数所限的一些不足。 forest_fn<-rodps.sql("select * from forest_big where tag=0 and prediction_result=1 order by prediction_score desc limit 10000;") forest_fn_detail<-sqldf("select * from dwtj2 where djxh in(select distinct djxh from forest_fn) order by djxh,yy") color<-c("Darkblue","Darkgreen","Blue","darkred","darkorchid1","coral","pink") p1<-ggplot(data=forest_fn_detail, aes(x=yy, y=pjgz/10000,group=djxh,colour=djxh)) + geom_line(size=1)+geom_point(colour="Black")+ xlab("年度")+ylab("平均工资(万)")+ ggtitle('随机森林False-Negative误判企业历年平均工资折线图')+ theme(legend.position="top",text=element_text(size=20))+ #scale_x_discrete(labels=levels(hstj7$rsfd))+ scale_colour_manual(name = "登记序号",values=color) print(p1) forest_fn_detail2<-forest_fn_detail[which(forest_fn_detail$djxh!='10124404000003614930'),] p1<-ggplot(data=forest_fn_detail2, aes(x=yy, y=pjgz/10000,group=djxh,colour=djxh)) + geom_line(size=1)+geom_point(colour="Black")+ xlab("年度")+ylab("平均工资(万)")+ ggtitle('随机森林False-Negative误判企业历年平均工资折线图')+ theme(legend.position="top",text=element_text(size=20))+ #scale_x_discrete(labels=levels(hstj7$rsfd))+ scale_colour_manual(name = "登记序号",values=color) print(p1) #保存工作区 save.image("d:/temp/GeShuiDaShuJu2.RData",compress=TRUE)