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

杨利霞

2020 全国大学生数学建模竞赛C题思路+代码
8 n1 G$ i% [& }/ ^, f" V7 t0 h题目链接：https://cloud.189.cn/t/ri2uUb7BRVJr


前言
: P- {5 z7 t% R, L5 ~% E
/ p1 c1 p7 Q( s/ b9 e+ V
; |4 C% ~1 v1 f6 w8 ?
  `$ L( b* K4 b
又是一年数据挖掘题型，第一次接触这种题型还是在去年的mathorcup上，这种题的难度就在于指标的建立和数据的处理上。后面会出一份关于数据挖掘题型，我的相关经验，常用的工具和代码。

- G! I9 o0 s5 r. _$ D# f. \3 }/ E
( G7 s3 j6 U: h下面的一，二问实际都在解决
' H1 h0 [! j4 {/ h( G
2 @) ]  ?0 q/ H% P: D; S  h( l
贷不贷款？
; h$ x* \: y* a% v! {4 \贷款金额多少？
数据清洗
7 r6 k0 @% K% z" x: X( g* R" q这道题的附件数据没有出现缺省或者异常数据，因此对于数据的预处理，更多的是根据问题的需求来做的。

6 y( ?- J4 R! \9 p0 S
3 z9 i; l6 Z$ Y( r8 d/ B9 J4 m$ e+ P将是否违约，违约设置为1，不违约设置为0
" u2 ?$ q$ \6 a
% n+ K4 o0 {" Q  }8 ?
信誉等级ABCD分别对应4,3,2,1


2 F  B/ h& V+ m' h  {& w; y发票状态，有效发票为a，作废发票为b
- t' g. ]2 ?+ O3 k/ W$ g/ `# f+ G
3 H" ], Y& o5 b0 N
, s" L3 k6 I4 e+ m9 k& e8 Q% L# R- `我将销项和进项所有数据，以公司代码为区别，提取到了不同的sheet当中，对于该公司有效发票数，作废发票数，负数发票数，方便对数据观察。
* b' b0 C0 w. Y4 C, {5 \7 k
' ?  `. F8 U: R5 ^1 k8 l
# 遍历所有sheet数据
( M6 V4 \& D- K" h! N! C8 @for xsn in sn.sheet_names[1:]:
7 g" E+ p3 U, d  q* W- B0 s    # 读取文件
. H9 q3 ^8 X- F2 c    datas = pd.read_excel(file_pos, sheet_name=xsn)
    datas['date']=pd.to_datetime(datas['date'],format='%Y/%m/%d')
    datas.set_index('date', drop=True)
    # 找到全部公司名称代号
, ]) W+ [/ c% w: r- R: [    code_list = list(set(list((datas['code']))))
! O* t" I5 @8 R) V( c  V: T' J! z    for name in code_list:
3 w& I9 d/ S, H" y        tmp_datas = datas[datas['code'] == name]
        tmp_datas.index = range(len(tmp_datas))
* U' R+ v. h5 l: j* [( q        # 转换日期未object类型
5 _  l: G' _2 p2 O0 [6 b  p        tmp_datas['date'] = [x.strftime('%Y/%m/%d') for x in tmp_datas['date']]
        count1 = tmp_datas['tax_status'].value_counts()
% y2 }0 ?# P3 ]$ W8 P        tmp_datas['a_count'] = list(count1)[0]
, E9 |5 x+ E; @1 X& c* ^6 x8 k        if(len(count1) > 1):
) u: {* T5 A( V! I# S            tmp_datas['b_count'] = list(count1)[1]
$ V2 I  c" _$ B. X0 A        tmp2 = tmp_datas[tmp_datas['cost'] < 0]
' d) ^2 e  p. b; ~& X5 N3 A        tmp_datas['neg_value_tax'] = len(tmp2)
        if xsn == sn.sheet_names[1]:
            tmp_datas.to_excel(writer1,sheet_name=name,index=False)
: q7 b$ `+ |. R! U        else:
            tmp_datas.to_excel(writer2,sheet_name=name,index=False)
负数发票：在之前购买的物品，并开具了相关正向发票，后来退货所以开具了值为负数的发票，抵消前面正数发票的值。

) d, m7 H+ {- h2 d
提取到信息：


6 B  I- Z8 n+ u7 y2 b: i部分公司数据记录很少，或者时间跨度大，需要综合数据指标，抵消数据数量和跨度大的影响
: Q" q0 k2 o  R- R( L有些负数发票，在之前找不到对应的正数发票，可能是因为在数据记录日期之前购买的，在之后退款，因此在附件中找不到记录。
# \- }* Q/ ?" [: d) V% y) j; Q问题一
建立指标
进项发票作废率，进项负数发票率，进项每月平均交易额，进项每月交易次数，
' C3 r. c5 w. N6 M

销项发票作废率，销项负数发票率，销项每月平均交易额，销项每月交易次数，销售收入增长率


提取出相关指标到附件
9 ?4 g4 V& c# @. n

for xsn in sn.sheet_names[1:]:
2 a3 {; u( d! ]/ T2 j4 H    # 读取文件
$ n4 E" K1 S7 D, b    datas = pd.read_excel(file_pos, sheet_name=xsn)
# E; L' H1 F) X: W0 `, L    code_list = list(set(list((datas['code']))))
    for name in code_list:
        tmp_datas = datas[datas['code'] == name]
3 z4 A& ~9 T$ x+ z# r2 `) E        tmp_datas.index = range(len(tmp_datas))
        insert_datas.append(name)
! \! R- T0 b4 ^2 p8 u* N        # 作废数
0 i5 E; E3 a0 c$ S4 E1 d: V        cacel_count = len(tmp_datas[tmp_datas['tax_status'] == 'b'])
        # 有效数
        valid_count = len(tmp_datas[tmp_datas['tax_status'] == 'a'])
* Y/ R( x  A- ~& i2 M        # 发票作废率
        count1 = (cacel_count / (cacel_count + valid_count))*100
2 L2 b) B9 X4 F8 N( l        # 负数发票数
        neg_count = len(tmp_datas[tmp_datas['cost'] < 0])
, u6 b+ M5 f& q3 C2 r1 N. }! c        # 负数发票率
4 g3 o0 u; @6 }; |        count2 = (neg_count / valid_count) * 100
        # 转换时间
) S( e' o& X& m" V        tmp_datas['date'] = [x.strftime('%Y/%m/%d') for x in tmp_datas['date']]
        # 时间最大值
! `5 ?3 d) }  U8 o+ s& t! G        max_time = tmp_datas.iloc[0:,1].max()
        # 时间最小值
        min_time = tmp_datas.iloc[0:,1].min()
6 ^- S/ A6 o( ^* n" ^, i) }1 P        # 时间差
7 i2 o8 o9 [2 H/ `        diff_time = months(max_time, min_time) + 1
        # 有效票
; t# n9 C+ ]! C  C& B6 b        valid_tax = tmp_datas[tmp_datas['tax_status'] == 'a']
2 ?: O, |% E" k8 A! x& o) z. j        # 平均月交易额
        avg_money = valid_tax['totle_cost'].sum() / diff_time
2 i. V( P9 c0 ^( e# p6 q+ n        # 平均每月交易次数
! p) d3 A8 ~" @, w' ]8 b& f        trans_count = len(tmp_datas) / diff_time
- S) D( M! S5 o. n9 L& ~& z; n        insert_datas += [count1, count2, avg_money, trans_count,]
. J, Q1 P7 d5 }" d7 g, m( r) {4 ~        if flag:
: y- A( Y! Q+ j2 @            df1.loc[len(df1)] = insert_datas
5 v- S. a, k  j/ S. [1 T) D; k            df1.to_excel(writer1,sheet_name='进项信息',index=False)
        else:
/ n* z- q+ t, G( v; R8 H9 ?) J            merge_time = tmp_datas.groupby(tmp_datas['date']).sum()
: A) m+ Q- Y8 ]. e+ R( G            # 销售收入增长率
8 W- a0 G6 [! N- k5 p            income_info= list((merge_time['cost'] - merge_time['cost'].shift(1)).fillna(1))
$ s8 b* a2 `1 T            diff_time_day = days(max_time,min_time)
            income_tax = (sum(income_info) / diff_time_day)*100
& N6 ~. y: Z0 d8 C0 _2 X; c8 R            insert_datas.append(income_tax)
. h' E' I. p, ?% D            df2.loc[len(df2)] = insert_datas
            df2.to_excel(writer1,sheet_name='销项信息',index=False)
        insert_datas = []
4 z; ^! s/ W+ c6 e    flag = False
! n' d' V: ~% G' k: E5 p. D( p
9 d  v% k& G3 A0 I3 P+ Y6 O" m9 L

# Z: b* O' x( T
1 `) P1 Z; W+ e7 U- [并将是否违约插入到最后一列


: }% ^$ ]5 D- A5 @1 J& `$ w# 提取是否违约的列表
/ M; ~/ x/ s+ Jm = []
5 L- O& |% |$ W( N4 k- i# qfor name in code_list:
% Z0 E/ m8 [- N$ Z8 t. u' C    m.append(datas[datas['code']==name]['break_contract'].tolist()[0])
    df1.loc[:,len(df1)] = m
2 t3 R5 F. Q! T  [: x5 d; T* }    df1.to_excel(writer3,sheet_name='sheet1',index=False)
5 I+ P# o/ d+ Y0 l' Y+ G# `建立模型
Logistics违约率预测模型
- _+ \+ q4 A( K6 m$ Z2 S- \3 ^. B使用Logistics违约预测模型，代入所有的指标数据为自变量，是否违约为因变量，预测出违约率。
7 g# Z! r1 b: A

8 g! {  ~+ R1 M7 M" zX=datas[['进项发票作废率','进项负数发票率','进项每月平均交易额','进项每月交易次数','销项发票作废率','销项负数发票率','销项每月平均交易额','销项每月交易次数','销售收入增长率']]
y=datas['是否违约']
X_train, X_test, y_train, y_test = train_test_split(X,y, test_size=0.2,random_state=2020)
3 `5 O! g) ^3 `. |0 b  p( `$ g) n+ Q# gX_validation, X_test, y_validation, y_test = train_test_split(X_test,y_test, test_size=0.1,random_state=2020)
model = LogisticRegression()
* L1 {' H/ b( cmodel.fit(X_train,y_train)
a=model.predict_proba(X_validation)
/ z3 D4 R" R% h! |+ ~8 ^4 g) Hresult=[]
( O  W% p, W9 C. n& m6 r5 |. `for i in range(len(a)):
0 z. X0 d  S7 _7 i* m( l4 qif a[1]>0.5:
    result.append(1)
7 w* N" J  T- w" L  selse:
    result.append(0)
% K$ ^0 k9 |& P: K, A( Z6 V3 hfrom sklearn import metrics
print('误差: %.4f' % (1-metrics.recall_score(y_validation,result,average='weighted')))
最终得到一张我们的分析表格



9 n& _/ n6 V  v


, N1 ?4 F* U# |9 `通过预测是否违约，我们就能解决贷不贷款的问题。
' W' v/ b$ l6 u+ L

9 O; T9 S2 u7 L3 H贷款金额
贷款金额的确认，根据该公司不违约率在所有公司中的权重，乘以总贷款金额确认：


$ |7 n+ w8 w2 b$ j
" e9 o- s/ c6 F/ t' q9 M

' [& J) Z* d3 D; `ri=1−Zi∑123j=1(1−Zi)×M
ri=1−Zi∑j=1123(1−Zi)×M
8 j% T- s, E9 f因此，我们得到的贷款金额是违约率和贷款总金额组成的关系式，这在第二问中能起到重要作用。
! Q" }% R/ b$ c# Z8 f

贷款年利率
6 U7 F  |6 f1 {* ?/ T绘制出年利率与客户流失率图，可以分析出两者应该是有关系的。利用SPSS拟合出不同信誉等级，年利率与客户流失率的关系式。


信誉等级        R平方        关系式
A        0.9977        y = 37.97x^3-258.57x^2+640.944*x -1.121
, J/ [7 k; w3 q: cB        0.9982        y = 33.995x^3-225.051x^2+552.829*x-1.017
C        0.9982        y = 32.157x^3-207.386x^2+504.717*x-0.973
银行获利=贷款金额x贷款年利率x(1-利率对于信誉评级客户流失率)
. p, f5 Q! _# {& q. ^! P
( C  n7 x; ?. \) a% d& v" F
在贷款金额确认，贷款年利率范围在0.4~1.5的情况下，利用上面拟合的关系式，我们能够暴力跑出最优年利率。

+ \! ~, p. r; t
double turnover_rate(double x, char ch) {
        double y = 0, result = 0;
7 t0 ~8 W. f9 _4 U. P0 ^- q        switch (ch) {
        case 'A':
% S6 x4 ]2 {1 M2 k; ^2 ]                y = 37.969520 * pow(x, 3) - 258.570452 * pow(x, 2) + 640.944427 * x - 1.121484;
                result = x * (1 - y / 100.0);
                break;
        case 'B':
                y = 33.994698 * pow(x, 3) - 225.050538 * pow(x, 2) + 552.829151 * x - 1.016503;
, G- O) O6 I) c- P                result = x * (1 - y / 100.0);
                break;
& V9 u% v: z& d8 q        case 'C':
                y = 32.156864 * pow(x, 3) - 207.385880 * pow(x, 2) + 504.716993 * x - 0.973497;
1 \3 N- b6 c- N4 b  u) J( R, C! P: r                result = x * (1 - y / 100.0);
                break;
0 r+ U; Z4 l' p! @        default:
                cout << "输出有误！"  << ch << endl;
- ?4 p2 {* a4 m( D& V        }
. S. k1 a9 n2 i2 x5 Q. ~
0 ~9 }9 ~0 T; I* |        return result;
7 q9 ?8 O& F, v# A5 Z! T. C}
8 }/ W; g" `$ }, v, w( Q
! j9 B: X2 v% j$ g# Z, ]+ c

9 Q8 F) Z( `" K8 {: S
问题二
利用代码，重新计算出各指标数据
代入Logistics违约率预测模型，预测出各公司的违约率
根据标准普尔评级建立，主标尺，对不同违约率进行A~D等级划分，信誉等级D不予贷款
将违约率代入，之前得到的公式，得到具体贷款金额
最优年利率沿用上一问
# 信用等级
' A& S# {" ^8 ?cs = []
& d, V/ Z% F; G% T- U& S# 最优年利率,客户流失率,利率值
tax = []
for i in m:
' l- R; |. [/ _9 Z0 l    if i <= 0.0069264:
) e* n0 t' p  [, ~2 d        cs.append('A')
  o# ]9 l; M0 N* t1 @; V# G) i+ S        tax.append([0.083,0.503173,0.0412366])
    elif i > 0.0069264 and i <= 0.22619:
4 R) \* D# }' Q  L        cs.append('B')
        tax.append([0.097,0.505215,0.0479942])
( S8 Q- r2 F' ]9 V" Y    elif i > 0.22619 and i <= 0.509915:
        cs.append('C')
        tax.append([0.1069,0.506501,0.052755])
    elif i > 0.509915:
( C% a9 k. j6 X3 a( |4 G( \/ ^7 V        cs.append('D')
5 D7 G" i- a. d% d& b; |1 X        tax.append([0.15,0,0])
% i% [! E! K/ k: P0 K8 Z. C    else:
        print('违规')
. o9 d( {8 k3 O& I/ {        
% h; D; Z9 G# A4 Cparr = []
8 I' g2 }# d2 G1 F4 Pfor arr in list(a):
, ^% J3 R% W8 y. I; h; H    parr.append(list(arr)[0])
: C$ g3 R) M0 w  Y' o: isum_val = sum(parr)
- G' D' Q# e/ y! r' v) u8 Samount = []
for ival in parr:
: i1 }" u' P* ]( w* A    tmp = ival / sum_val * 100000000
    if ival < 1 - 0.509915:
+ f- @. s/ I+ Y/ ~  u2 e5 L, g7 ?            amount.append(0)
! Y; `( V% {( m; b    else:
: y+ I2 O+ x8 _: g- L1 I            amount.append(tmp)


$ T! x& G2 A; X# w  s
' y1 V* N. x) }  ?
可以看到，贷款金额也都在10w~100w之内。
& ^: w4 s% c$ L) e2 I9 z& `- J8 S- U

$ y( T1 a* Y/ R7 I! a' i问题三
这一问，我们做得有些匆忙了，有其他想法的可以按照自己的想法做做，这里只拿我们的做参照。

$ z; i8 C+ c: b* G# k( y# f+ I$ g
疫情对公司影响最大的就是每月平均销售额和每月平均销售数量，因此，


对每月平均销售额和每月平均销售数量，分别取随机数，数量取10w组，其他指标数据值不变，是否违约数据根据第二问结果，设为初始值
$ j( c+ ~4 `8 P! Y代入Logistics模型中，预测出每一组的违约率
: w4 U/ C3 Q; G7 y6 }( `) ?判断每个公司违约率的变化情况，根据变化情况来增/减贷款金额和年利率。
rand_num = pro_rand()

! @6 G# S4 y. c' g* Qdata_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':[]}
8 W# F9 y$ h+ y  P5 Idf = DataFrame(data_form)

for val in rand_num:
        a = rate_func(val)
2 s8 S9 x3 p# T2 i0 D+ m        rate_list = get_rate(a)
, m+ l0 r2 l9 e8 r9 {! ]        df.loc[len(df)]  = rate_list
        df.to_excel(writer,sheet_name='违约率变化',index=False)
" V4 G; D. T! R( y* a. d; c4 w销售数量和销售金额的随机数范围是：0~MAX


; ?/ N! p; a: S7 e

————————————————
% E* |  i5 t8 R  N& Q- q- {' W版权声明：本文为CSDN博主「Hk_Mayfly」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
' Y# C+ I# V( ]2 l原文链接：https://blog.csdn.net/qq_39542714/article/details/108592087

5 r2 _& v9 |& u, S3 X5 z
: |- S( j% i0 }. l3 b# r: r# s
; W- V+ K6 j- O6 B. F" H) m