2020 全国大学生数学建模竞赛C题思路+代码

杨利霞

2020 全国大学生数学建模竞赛C题思路+代码
题目链接：https://cloud.189.cn/t/ri2uUb7BRVJr


- l0 n( a, f$ \: s) W前言

" O% k/ A2 X, l5 W


又是一年数据挖掘题型，第一次接触这种题型还是在去年的mathorcup上，这种题的难度就在于指标的建立和数据的处理上。后面会出一份关于数据挖掘题型，我的相关经验，常用的工具和代码。
2 @" a9 O. Y; l& z- i. o2 H

下面的一，二问实际都在解决

. |- Q: S" N9 Q# B5 {+ L5 S; j/ q
, U5 S" [$ |9 c2 t6 P贷不贷款？
贷款金额多少？
数据清洗
8 Z5 q- W) {& t- O8 H2 e; D) [7 t这道题的附件数据没有出现缺省或者异常数据，因此对于数据的预处理，更多的是根据问题的需求来做的。


将是否违约，违约设置为1，不违约设置为0
9 I8 g- ^& U. x

/ Z% i  {0 c3 D  Y' `" h7 q3 M信誉等级ABCD分别对应4,3,2,1


发票状态，有效发票为a，作废发票为b

+ r, j! |8 C% H4 l7 |
我将销项和进项所有数据，以公司代码为区别，提取到了不同的sheet当中，对于该公司有效发票数，作废发票数，负数发票数，方便对数据观察。

  |# \  k: d2 `( K8 j9 T  Z
# 遍历所有sheet数据
for xsn in sn.sheet_names[1:]:
    # 读取文件
    datas = pd.read_excel(file_pos, sheet_name=xsn)
0 i6 i: b3 [1 `6 o    datas['date']=pd.to_datetime(datas['date'],format='%Y/%m/%d')
& F% O2 o  P- A8 h8 j5 \: W    datas.set_index('date', drop=True)
/ z" m! {& l; G3 x: |4 m2 w  V    # 找到全部公司名称代号
    code_list = list(set(list((datas['code']))))
    for name in code_list:
" x6 n0 G3 \6 l8 X9 r        tmp_datas = datas[datas['code'] == name]
2 e' m0 Q. z) {3 x# m, L9 Y7 }9 O        tmp_datas.index = range(len(tmp_datas))
$ Y% N  k9 D$ F2 D        # 转换日期未object类型
        tmp_datas['date'] = [x.strftime('%Y/%m/%d') for x in tmp_datas['date']]
        count1 = tmp_datas['tax_status'].value_counts()
        tmp_datas['a_count'] = list(count1)[0]
& W5 |- M) t5 g6 ?1 R        if(len(count1) > 1):
            tmp_datas['b_count'] = list(count1)[1]
        tmp2 = tmp_datas[tmp_datas['cost'] < 0]
        tmp_datas['neg_value_tax'] = len(tmp2)
9 }7 F, N% J* Z6 Z        if xsn == sn.sheet_names[1]:
            tmp_datas.to_excel(writer1,sheet_name=name,index=False)
+ v9 t) W7 r% w% ^- y5 W+ ~        else:
            tmp_datas.to_excel(writer2,sheet_name=name,index=False)
负数发票：在之前购买的物品，并开具了相关正向发票，后来退货所以开具了值为负数的发票，抵消前面正数发票的值。
( V; ]- G- `* ?& r6 s2 X& n

提取到信息：


) ^0 G$ f9 j# w. E% k部分公司数据记录很少，或者时间跨度大，需要综合数据指标，抵消数据数量和跨度大的影响
+ J( ~4 w5 K" `$ g; k# D有些负数发票，在之前找不到对应的正数发票，可能是因为在数据记录日期之前购买的，在之后退款，因此在附件中找不到记录。
& `9 }  a( w* C& S+ Y问题一
. F8 g: g. T5 @: J+ c4 H! i1 l1 u建立指标
7 @" t4 r7 g) w* [. G进项发票作废率，进项负数发票率，进项每月平均交易额，进项每月交易次数，

7 J; ~9 P1 h/ d6 I) m' E
' M# H6 Z; ]/ ?9 W- K! [0 X销项发票作废率，销项负数发票率，销项每月平均交易额，销项每月交易次数，销售收入增长率

2 {+ w3 Y, B8 ]  k" I
提取出相关指标到附件


for xsn in sn.sheet_names[1:]:
0 C5 n; I9 r5 C: p6 ~    # 读取文件
    datas = pd.read_excel(file_pos, sheet_name=xsn)
- a, Q- d( b6 D; z, x    code_list = list(set(list((datas['code']))))
    for name in code_list:
        tmp_datas = datas[datas['code'] == name]
! y- e9 l7 U, z" f7 u        tmp_datas.index = range(len(tmp_datas))
        insert_datas.append(name)
5 A. E6 A: o: ?0 |" {' S3 H3 C        # 作废数
        cacel_count = len(tmp_datas[tmp_datas['tax_status'] == 'b'])
& g9 h7 ~. Z8 Q0 w7 m* K        # 有效数
        valid_count = len(tmp_datas[tmp_datas['tax_status'] == 'a'])
        # 发票作废率
* g; n5 O' t# B' I) }7 F' u        count1 = (cacel_count / (cacel_count + valid_count))*100
        # 负数发票数
        neg_count = len(tmp_datas[tmp_datas['cost'] < 0])
3 u- }! `8 ^5 V, q4 A        # 负数发票率
' h5 W/ V/ C0 [& Y2 u! ]        count2 = (neg_count / valid_count) * 100
        # 转换时间
        tmp_datas['date'] = [x.strftime('%Y/%m/%d') for x in tmp_datas['date']]
        # 时间最大值
' V# @* M3 a! b& }5 [6 C        max_time = tmp_datas.iloc[0:,1].max()
* [, t' s! d8 r$ `5 R8 X        # 时间最小值
        min_time = tmp_datas.iloc[0:,1].min()
        # 时间差
+ k8 K( g  f, t/ K! m        diff_time = months(max_time, min_time) + 1
        # 有效票
        valid_tax = tmp_datas[tmp_datas['tax_status'] == 'a']
6 `3 }# l( ?1 j  j  r        # 平均月交易额
/ c% l- m; B* i        avg_money = valid_tax['totle_cost'].sum() / diff_time
        # 平均每月交易次数
, K8 }- _& i9 u( _5 Y        trans_count = len(tmp_datas) / diff_time
9 ?1 e* w3 H% G        insert_datas += [count1, count2, avg_money, trans_count,]
        if flag:
            df1.loc[len(df1)] = insert_datas
            df1.to_excel(writer1,sheet_name='进项信息',index=False)
, b8 r* N. l: a! d' @* [        else:
            merge_time = tmp_datas.groupby(tmp_datas['date']).sum()
+ ^2 S8 C. Z5 U6 G+ R' U5 Z            # 销售收入增长率
            income_info= list((merge_time['cost'] - merge_time['cost'].shift(1)).fillna(1))
' e" X7 b2 _  o9 G            diff_time_day = days(max_time,min_time)
            income_tax = (sum(income_info) / diff_time_day)*100
( X9 \8 N. K5 |  |            insert_datas.append(income_tax)
            df2.loc[len(df2)] = insert_datas
            df2.to_excel(writer1,sheet_name='销项信息',index=False)
        insert_datas = []
    flag = False
6 g/ u" e& ]* w
. }  m: W6 J) c$ u) P3 {
$ {4 y2 x. \; A: G9 K9 r8 l" }
3 F2 ^9 O0 U# Z- C+ v% [. o
! n7 Q! y5 v. k9 [8 \. L并将是否违约插入到最后一列
' k. Q* U5 t0 W1 J

# 提取是否违约的列表
( a' z* \8 l) J5 I+ T; rm = []
for name in code_list:
; ^  m5 t5 D* r2 [- T4 D3 w$ S    m.append(datas[datas['code']==name]['break_contract'].tolist()[0])
  a) q1 O  J% @- a- j/ l    df1.loc[:,len(df1)] = m
' \, O" ~" [! ^4 O$ \) p    df1.to_excel(writer3,sheet_name='sheet1',index=False)
6 f, [. \# S% S% X; ~9 W建立模型
8 w% M, [0 a- L, ~' MLogistics违约率预测模型
使用Logistics违约预测模型，代入所有的指标数据为自变量，是否违约为因变量，预测出违约率。
( F+ P& N1 h) ^

X=datas[['进项发票作废率','进项负数发票率','进项每月平均交易额','进项每月交易次数','销项发票作废率','销项负数发票率','销项每月平均交易额','销项每月交易次数','销售收入增长率']]
y=datas['是否违约']
+ |. i) D5 C* _$ [X_train, X_test, y_train, y_test = train_test_split(X,y, test_size=0.2,random_state=2020)
. i& G' [+ e4 ^: d! m7 f, y8 v8 ~X_validation, X_test, y_validation, y_test = train_test_split(X_test,y_test, test_size=0.1,random_state=2020)
; R5 e/ \+ \0 S  U, V; vmodel = LogisticRegression()
model.fit(X_train,y_train)
  O7 C- V4 G: e# o* I5 Ba=model.predict_proba(X_validation)
: ^6 A/ [: h; D  i6 v" _result=[]
$ O! ?$ E3 Q# x% \$ xfor i in range(len(a)):
if a[1]>0.5:
1 y- ?/ b  k3 n) x& t    result.append(1)
else:
    result.append(0)
6 Y  y0 R" P& f0 z3 e0 q3 Hfrom sklearn import metrics
% ^% E" h4 i/ a% b- ~3 E4 sprint('误差: %.4f' % (1-metrics.recall_score(y_validation,result,average='weighted')))
最终得到一张我们的分析表格
2 |8 e$ r, \% z* N


& x8 ^2 p% b* J& a5 |, |; l

- w2 h1 a, n& G8 }' _5 X
通过预测是否违约，我们就能解决贷不贷款的问题。
" V1 _( a. o. u1 ]
5 t/ t% Q9 H# R+ G0 A- F
9 B. X1 |5 G# b# \, V贷款金额
贷款金额的确认，根据该公司不违约率在所有公司中的权重，乘以总贷款金额确认：


7 K" h& b$ M! M& O0 ]* p

4 `$ M+ b+ D; Q# K  j" X
3 T- M; h; s3 k# e0 S. Sri=1−Zi∑123j=1(1−Zi)×M
1 m- c, b  ?) e- C! @! h' bri=1−Zi∑j=1123(1−Zi)×M
3 Y8 a. I1 A6 N6 c因此，我们得到的贷款金额是违约率和贷款总金额组成的关系式，这在第二问中能起到重要作用。


贷款年利率
% Z! ?9 d& W8 Z( x/ ~绘制出年利率与客户流失率图，可以分析出两者应该是有关系的。利用SPSS拟合出不同信誉等级，年利率与客户流失率的关系式。
+ m- `$ i8 E# g6 X  `! o/ Y8 m7 M

信誉等级        R平方        关系式
A        0.9977        y = 37.97x^3-258.57x^2+640.944*x -1.121
B        0.9982        y = 33.995x^3-225.051x^2+552.829*x-1.017
  G$ `: a1 z) `3 L. Z: sC        0.9982        y = 32.157x^3-207.386x^2+504.717*x-0.973
% G; I1 U- T% k9 t银行获利=贷款金额x贷款年利率x(1-利率对于信誉评级客户流失率)

; [3 E8 O( w1 l% ]/ {% N
0 \0 G2 V- p& Q/ U3 \% T9 H6 I在贷款金额确认，贷款年利率范围在0.4~1.5的情况下，利用上面拟合的关系式，我们能够暴力跑出最优年利率。
% Y+ G$ b* T# x8 \2 H
: {" }* E6 l+ ^% {; e5 T
0 b4 e; Q: F- I  R, Udouble turnover_rate(double x, char ch) {
        double y = 0, result = 0;
        switch (ch) {
9 Y2 d- j8 M# Z- z& H2 P        case 'A':
, K# o- [0 h; y9 R" |* c+ I- H6 K* \                y = 37.969520 * pow(x, 3) - 258.570452 * pow(x, 2) + 640.944427 * x - 1.121484;
                result = x * (1 - y / 100.0);
: r; ?  r7 w' J2 D9 [7 k2 P                break;
5 Z$ K! r5 o" s7 N% V3 `) A        case 'B':
                y = 33.994698 * pow(x, 3) - 225.050538 * pow(x, 2) + 552.829151 * x - 1.016503;
                result = x * (1 - y / 100.0);
                break;
        case 'C':
                y = 32.156864 * pow(x, 3) - 207.385880 * pow(x, 2) + 504.716993 * x - 0.973497;
                result = x * (1 - y / 100.0);
                break;
        default:
  K" ]3 K. d/ U                cout << "输出有误！"  << ch << endl;
        }

6 C" L- E/ r' k; c        return result;
}
5 K& d% E. W; Q6 E% `# z  j" A; K
9 q- S# K9 R5 z. S

( X2 [4 h6 V2 _) e5 B% z& ^- M
问题二
利用代码，重新计算出各指标数据
代入Logistics违约率预测模型，预测出各公司的违约率
" e2 k; e$ T0 u6 }根据标准普尔评级建立，主标尺，对不同违约率进行A~D等级划分，信誉等级D不予贷款
将违约率代入，之前得到的公式，得到具体贷款金额
4 |% @1 e9 @1 i0 q# y% E) m; Q最优年利率沿用上一问
* s7 D* F! ~7 q  s- Q, X- I# 信用等级
cs = []
# n% f" A* o+ G; u# 最优年利率,客户流失率,利率值
tax = []
for i in m:
  {7 l0 F3 j) o0 J    if i <= 0.0069264:
4 z/ p5 i7 ?8 N: R$ t        cs.append('A')
        tax.append([0.083,0.503173,0.0412366])
7 E$ W6 _# g  `    elif i > 0.0069264 and i <= 0.22619:
: Y9 j; D; N# C& c+ R+ K5 F2 _' S6 T3 |3 @        cs.append('B')
' U, i4 Q4 `5 ^* S) o4 O/ R* q        tax.append([0.097,0.505215,0.0479942])
    elif i > 0.22619 and i <= 0.509915:
        cs.append('C')
        tax.append([0.1069,0.506501,0.052755])
. @" i9 D" e" k* w+ M4 n- ?    elif i > 0.509915:
7 ~- _* R6 G5 }' [; l* }5 x2 s6 |        cs.append('D')
) |/ \# E3 [# z        tax.append([0.15,0,0])
    else:
5 P( O9 z6 p2 T3 R- H        print('违规')
        
parr = []
for arr in list(a):
    parr.append(list(arr)[0])
5 h1 m7 ]1 _& N! F. C9 v1 Fsum_val = sum(parr)
amount = []
for ival in parr:
$ u) \" W  B6 _% A: C    tmp = ival / sum_val * 100000000
    if ival < 1 - 0.509915:
            amount.append(0)
0 t% b9 A( v& z3 Z7 r* t# ?    else:
            amount.append(tmp)
: ~3 c5 h/ J0 o. ]



) k4 T7 y# D8 s3 i7 a/ F可以看到，贷款金额也都在10w~100w之内。

& T% U" S: D) i' ?% o
问题三
; x3 {% |' E6 r7 L这一问，我们做得有些匆忙了，有其他想法的可以按照自己的想法做做，这里只拿我们的做参照。

( k( ], }( h( p) ~, z
疫情对公司影响最大的就是每月平均销售额和每月平均销售数量，因此，

! H9 P3 K, z$ O& p) v
$ x4 d8 p: |" p5 Z8 T2 [4 \+ {对每月平均销售额和每月平均销售数量，分别取随机数，数量取10w组，其他指标数据值不变，是否违约数据根据第二问结果，设为初始值
代入Logistics模型中，预测出每一组的违约率
) ~3 L4 _4 _+ ^) X8 p8 _判断每个公司违约率的变化情况，根据变化情况来增/减贷款金额和年利率。
) ^1 H; h  r' n& C  R) }- S6 [rand_num = pro_rand()
# v7 r1 O% x# \+ w2 X
1 G: v% P; @9 H" L3 Rdata_form = {'E377':[], 'E311':[], 'E297':[], 'E386':[], 'E233':[], 'E367':[], 'E194':[], 'E196':[], 'E249':[], 'E205':[], 'E159':[], 'E300':[], 'E346':[], 'E395':[], 'E360':[], 'E200':[], 'E195':[], 'E247':[], 'E317':[], 'E303':[], 'E135':[], 'E217':[], 'E298':[], 'E325':[], 'E335':[], 'E157':[], 'E287':[], 'E270':[], 'E286':[], 'E212':[], 'E261':[], 'E423':[], 'E390':[], 'E189':[], 'E129':[], 'E366':[], 'E191':[], 'E246':[], 'E406':[], 'E357':[], 'E349':[], 'E387':[], 'E372':[], 'E385':[], 'E209':[], 'E267':[], 'E414':[], 'E316':[], 'E341':[], 'E281':[], 'E206':[], 'E363':[], 'E166':[], 'E190':[], 'E140':[], 'E130':[], 'E225':[], 'E347':[], 'E356':[], 'E185':[], 'E192':[], 'E207':[], 'E234':[], 'E136':[], 'E379':[], 'E274':[], 'E383':[], 'E242':[], 'E361':[], 'E408':[], 'E198':[], 'E204':[], 'E413':[], 'E389':[], 'E253':[], 'E226':[], 'E231':[], 'E182':[], 'E318':[], 'E392':[], 'E275':[], 'E425':[], 'E388':[], 'E305':[], 'E155':[], 'E348':[], 'E400':[], 'E256':[], 'E351':[], 'E201':[], 'E345':[], 'E278':[], 'E306':[], 'E308':[], 'E215':[], 'E382':[], 'E407':[], 'E376':[], 'E291':[], 'E369':[], 'E139':[], 'E296':[], 'E260':[], 'E145':[], 'E257':[], 'E216':[], 'E125':[], 'E312':[], 'E354':[], 'E380':[], 'E326':[], 'E268':[], 'E186':[], 'E179':[], 'E223':[], 'E319':[], 'E422':[], 'E224':[], 'E151':[], 'E293':[], 'E302':[], 'E309':[], 'E378':[], 'E373':[], 'E364':[], 'E144':[], 'E162':[], 'E208':[], 'E399':[], 'E355':[], 'E197':[], 'E375':[], 'E368':[], 'E334':[], 'E172':[], 'E254':[], 'E352':[], 'E230':[], 'E313':[], 'E255':[], 'E273':[], 'E220':[], 'E169':[], 'E126':[], 'E237':[], 'E153':[], 'E164':[], 'E344':[], 'E284':[], 'E320':[], 'E328':[], 'E143':[], 'E338':[], 'E290':[], 'E251':[], 'E314':[], 'E330':[], 'E424':[], 'E183':[], 'E248':[], 'E181':[], 'E271':[], 'E403':[], 'E359':[], 'E324':[], 'E304':[], 'E213':[], 'E412':[], 'E337':[], 'E235':[], 'E371':[], 'E283':[], 'E391':[], 'E370':[], 'E173':[], 'E158':[], 'E171':[], 'E299':[], 'E310':[], 'E174':[], 'E295':[], 'E402':[], 'E419':[], 'E327':[], 'E294':[], 'E339':[], 'E203':[], 'E285':[], 'E410':[], 'E241':[], 'E152':[], 'E228':[], 'E245':[], 'E263':[], 'E404':[], 'E160':[], 'E292':[], 'E397':[], 'E329':[], 'E163':[], 'E301':[], 'E365':[], 'E178':[], 'E288':[], 'E175':[], 'E211':[], 'E252':[], 'E276':[], 'E156':[], 'E232':[], 'E321':[], 'E128':[], 'E150':[], 'E148':[], 'E161':[], 'E222':[], 'E142':[], 'E134':[], 'E238':[], 'E396':[], 'E239':[], 'E227':[], 'E147':[], 'E188':[], 'E243':[], 'E401':[], 'E362':[], 'E418':[], 'E272':[], 'E333':[], 'E405':[], 'E416':[], 'E358':[], 'E202':[], 'E259':[], 'E180':[], 'E409':[], 'E244':[], 'E394':[], 'E282':[], 'E210':[], 'E124':[], 'E398':[], 'E374':[], 'E277':[], 'E265':[], 'E421':[], 'E342':[], 'E322':[], 'E165':[], 'E340':[], 'E381':[], 'E331':[], 'E280':[], 'E393':[], 'E353':[], 'E221':[], 'E219':[], 'E323':[], 'E336':[], 'E199':[], 'E170':[], 'E137':[], 'E168':[], 'E350':[], 'E141':[], 'E214':[], 'E177':[], 'E250':[], 'E269':[], 'E167':[], 'E307':[], 'E420':[], 'E411':[], 'E262':[], 'E176':[], 'E218':[], 'E229':[], 'E315':[], 'E264':[], 'E127':[], 'E132':[], 'E154':[], 'E133':[], 'E332':[], 'E138':[], 'E258':[], 'E266':[], 'E236':[], 'E415':[], 'E146':[], 'E193':[], 'E187':[], 'E384':[], 'E417':[], 'E289':[], 'E279':[], 'E131':[], 'E184':[], 'E240':[], 'E149':[], 'E343':[]}
4 s% \2 r  ]# y8 t9 i8 B' ^df = DataFrame(data_form)

for val in rand_num:
' Z" I: q- K2 K. `7 W        a = rate_func(val)
/ z. w% y0 y) Q) q        rate_list = get_rate(a)
& N) O" i, \5 @9 ?- t7 b        df.loc[len(df)]  = rate_list
* c# S' C# l; y9 n5 y+ S        df.to_excel(writer,sheet_name='违约率变化',index=False)
销售数量和销售金额的随机数范围是：0~MAX

( A$ b) C& w; X3 O
; |4 b# v) {9 S0 {, b1 Z9 {
8 j8 B8 N/ I: ^+ a9 L$ J1 q, {
————————————————
版权声明：本文为CSDN博主「Hk_Mayfly」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_39542714/article/details/108592087
6 q5 H. K: z. |' l, r8 [
2 y( Q4 u$ E. r8 m

' R1 V$ W. p( s, Q% K/ F5 ]