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

杨利霞

2020 全国大学生数学建模竞赛C题思路+代码
% J1 H) ^  W9 ?0 Q, p题目链接：https://cloud.189.cn/t/ri2uUb7BRVJr

+ m/ `# l3 w$ h
0 z% ?- r0 c8 S/ s: |前言
6 g4 s0 ]: b5 e" F7 z

2 x/ n- d( n4 n( Z/ E  h, f+ d

( N( z6 ?' Q4 R又是一年数据挖掘题型，第一次接触这种题型还是在去年的mathorcup上，这种题的难度就在于指标的建立和数据的处理上。后面会出一份关于数据挖掘题型，我的相关经验，常用的工具和代码。
0 X3 G$ |5 |$ a0 `
$ S# l9 Q- ]! e1 F
8 I9 f8 b$ d" }' Y; y下面的一，二问实际都在解决
# s. y# e# z( i* R# R: K9 h
! ]- `! x" D: I% W' {' L
贷不贷款？
贷款金额多少？
7 }+ b8 z6 E, o) H' h$ r数据清洗
- x6 j9 G0 e* @这道题的附件数据没有出现缺省或者异常数据，因此对于数据的预处理，更多的是根据问题的需求来做的。

& Q9 b8 t% Q8 d, j% X% x2 l
9 H1 h9 L7 p. C, A% H将是否违约，违约设置为1，不违约设置为0


* f( s: s  J5 q* K/ n* G2 @! N9 k信誉等级ABCD分别对应4,3,2,1

; O! t' X4 ^8 l2 o& C3 w6 L
发票状态，有效发票为a，作废发票为b
* i$ h; Q( s6 \) S, o

我将销项和进项所有数据，以公司代码为区别，提取到了不同的sheet当中，对于该公司有效发票数，作废发票数，负数发票数，方便对数据观察。
6 W0 @3 @  Y0 y

# 遍历所有sheet数据
; ]. l' z7 N. V& x1 I  kfor xsn in sn.sheet_names[1:]:
8 g/ O' H6 n4 T" l    # 读取文件
% D0 |5 [( J; f% r7 G    datas = pd.read_excel(file_pos, sheet_name=xsn)
    datas['date']=pd.to_datetime(datas['date'],format='%Y/%m/%d')
+ O1 X' }6 a6 h9 {    datas.set_index('date', drop=True)
    # 找到全部公司名称代号
: B9 K# @0 H+ n9 j; F) b    code_list = list(set(list((datas['code']))))
2 w6 K% h* r% h: W7 _# C    for name in code_list:
) l4 J  ~( c2 R" `6 Y2 h* t- Q* p        tmp_datas = datas[datas['code'] == name]
        tmp_datas.index = range(len(tmp_datas))
4 e; x+ M; j9 t; `7 y2 L4 T        # 转换日期未object类型
. Y" \+ K7 i7 x, @, d        tmp_datas['date'] = [x.strftime('%Y/%m/%d') for x in tmp_datas['date']]
        count1 = tmp_datas['tax_status'].value_counts()
2 Z0 M2 K) \  s, C1 {        tmp_datas['a_count'] = list(count1)[0]
        if(len(count1) > 1):
            tmp_datas['b_count'] = list(count1)[1]
( l( ?- L) ]6 T: \        tmp2 = tmp_datas[tmp_datas['cost'] < 0]
        tmp_datas['neg_value_tax'] = len(tmp2)
        if xsn == sn.sheet_names[1]:
' P& s9 W6 A7 P+ M$ g9 B9 q            tmp_datas.to_excel(writer1,sheet_name=name,index=False)
/ u5 i) v9 O" m        else:
$ O2 G, O% w7 b8 P7 z) ^# D            tmp_datas.to_excel(writer2,sheet_name=name,index=False)
负数发票：在之前购买的物品，并开具了相关正向发票，后来退货所以开具了值为负数的发票，抵消前面正数发票的值。


3 ]& [! t4 Q; p0 ~提取到信息：

- {3 L1 H! M5 Z6 z( I' \
% i! W1 w. [1 w# G; @部分公司数据记录很少，或者时间跨度大，需要综合数据指标，抵消数据数量和跨度大的影响
有些负数发票，在之前找不到对应的正数发票，可能是因为在数据记录日期之前购买的，在之后退款，因此在附件中找不到记录。
问题一
/ g( K8 l5 c$ m- `建立指标
进项发票作废率，进项负数发票率，进项每月平均交易额，进项每月交易次数，

0 t9 P- V# A3 m8 E, _
0 x+ O' i* L: e) _1 a9 n( k1 U销项发票作废率，销项负数发票率，销项每月平均交易额，销项每月交易次数，销售收入增长率
5 C; T6 G! H9 s5 Z4 C7 C, g0 |

9 C# J! \: {8 ]( W5 V& h! V提取出相关指标到附件
, o6 P, s7 \% {( S+ t2 w
" N+ N4 V- D: K' ~3 y
- o5 y6 a3 t6 j# ]7 }# ~for xsn in sn.sheet_names[1:]:
- M4 \6 B( e: D7 d    # 读取文件
$ V% _* Q* g+ O2 G    datas = pd.read_excel(file_pos, sheet_name=xsn)
    code_list = list(set(list((datas['code']))))
    for name in code_list:
        tmp_datas = datas[datas['code'] == name]
        tmp_datas.index = range(len(tmp_datas))
        insert_datas.append(name)
  x3 r3 f- t5 T: {: L1 T        # 作废数
. q0 g' Q: i( X* }8 i- s        cacel_count = len(tmp_datas[tmp_datas['tax_status'] == 'b'])
) r  [1 T' f. ?# W, x        # 有效数
! u3 y) i+ O$ ~        valid_count = len(tmp_datas[tmp_datas['tax_status'] == 'a'])
        # 发票作废率
        count1 = (cacel_count / (cacel_count + valid_count))*100
        # 负数发票数
        neg_count = len(tmp_datas[tmp_datas['cost'] < 0])
        # 负数发票率
5 L2 ]; c$ b2 B, {, ?( C( X, ^        count2 = (neg_count / valid_count) * 100
        # 转换时间
3 H! B; S1 S+ M! I+ Y2 u        tmp_datas['date'] = [x.strftime('%Y/%m/%d') for x in tmp_datas['date']]
        # 时间最大值
        max_time = tmp_datas.iloc[0:,1].max()
        # 时间最小值
0 Y) G* z: S4 T- p; O. o        min_time = tmp_datas.iloc[0:,1].min()
8 S  M* _5 s+ U7 x8 _/ f0 J4 W        # 时间差
        diff_time = months(max_time, min_time) + 1
) w8 \% y* q! a2 \! `2 w( i6 N        # 有效票
        valid_tax = tmp_datas[tmp_datas['tax_status'] == 'a']
. Q- h' O( s! x" s3 a4 ^        # 平均月交易额
8 W) }' M8 c& f+ A8 S3 q' _        avg_money = valid_tax['totle_cost'].sum() / diff_time
        # 平均每月交易次数
        trans_count = len(tmp_datas) / diff_time
        insert_datas += [count1, count2, avg_money, trans_count,]
  R% b) j/ u% _1 V. f        if flag:
            df1.loc[len(df1)] = insert_datas
            df1.to_excel(writer1,sheet_name='进项信息',index=False)
        else:
1 A' q. \) F3 S0 ]( z; ^( a1 @% D            merge_time = tmp_datas.groupby(tmp_datas['date']).sum()
            # 销售收入增长率
  x  Z5 |9 i5 N0 G# T0 r            income_info= list((merge_time['cost'] - merge_time['cost'].shift(1)).fillna(1))
            diff_time_day = days(max_time,min_time)
( M) S* T! f6 D( `3 N% h9 d  S            income_tax = (sum(income_info) / diff_time_day)*100
            insert_datas.append(income_tax)
            df2.loc[len(df2)] = insert_datas
' V- h8 o3 ]6 O% q            df2.to_excel(writer1,sheet_name='销项信息',index=False)
7 J+ [, i. z$ |  E+ _! Z4 w        insert_datas = []
    flag = False
4 z9 C; {! N% G


- [  a4 J0 N2 a
并将是否违约插入到最后一列
/ Q5 {) P7 V2 s4 S; m' Y* B3 _7 u; a! J
8 v5 ~6 X7 ^$ L. L
# 提取是否违约的列表
m = []
* G6 p9 J: b& O$ S. B$ y1 y2 Gfor name in code_list:
- |0 u$ i( |9 o# ]* ?6 a    m.append(datas[datas['code']==name]['break_contract'].tolist()[0])
    df1.loc[:,len(df1)] = m
    df1.to_excel(writer3,sheet_name='sheet1',index=False)
0 E& |6 k1 s; N/ J+ M3 z& B建立模型
Logistics违约率预测模型
# x: K1 V3 s, v& |" F: J使用Logistics违约预测模型，代入所有的指标数据为自变量，是否违约为因变量，预测出违约率。
6 ^% H$ l" q' m  ^
& f  k. }8 e  q5 V4 h
$ C1 Q/ n/ e$ ^( Y7 F$ L  a: JX=datas[['进项发票作废率','进项负数发票率','进项每月平均交易额','进项每月交易次数','销项发票作废率','销项负数发票率','销项每月平均交易额','销项每月交易次数','销售收入增长率']]
4 o1 K+ u$ \8 |" ~5 R% P" J6 Ey=datas['是否违约']
! g5 r( @7 H% d9 T, TX_train, X_test, y_train, y_test = train_test_split(X,y, test_size=0.2,random_state=2020)
/ c0 S8 J3 t" g. x7 eX_validation, X_test, y_validation, y_test = train_test_split(X_test,y_test, test_size=0.1,random_state=2020)
model = LogisticRegression()
model.fit(X_train,y_train)
a=model.predict_proba(X_validation)
result=[]
for i in range(len(a)):
7 S9 L! p* \: P% [! R* N2 k+ ]if a[1]>0.5:
; ]) E! {  @, r% G" N, w  w1 M    result.append(1)
% [6 E2 M  J$ G+ Selse:
+ H* Y8 e) a& i9 ]! f0 L    result.append(0)
from sklearn import metrics
- E& Y# l+ \: r7 l( E& D/ @5 kprint('误差: %.4f' % (1-metrics.recall_score(y_validation,result,average='weighted')))
最终得到一张我们的分析表格
' v; M/ q2 R% S: k

, S" t9 I7 i% s
  {0 A3 {& Z3 ?+ ?. v


通过预测是否违约，我们就能解决贷不贷款的问题。


贷款金额
贷款金额的确认，根据该公司不违约率在所有公司中的权重，乘以总贷款金额确认：

; B. }+ r: `# X6 @( `
  |, \, B% P5 ~" g8 I

) A* u2 h/ \& A9 [8 X$ `: V2 o
! p2 l! o+ ?3 P9 A) P/ m- F  ?ri=1−Zi∑123j=1(1−Zi)×M
ri=1−Zi∑j=1123(1−Zi)×M
' [2 P4 g3 l4 w因此，我们得到的贷款金额是违约率和贷款总金额组成的关系式，这在第二问中能起到重要作用。

5 i; a: h& {7 [% q- t1 }5 ~
贷款年利率
% ^0 j+ e/ Y0 h' F1 L& Z绘制出年利率与客户流失率图，可以分析出两者应该是有关系的。利用SPSS拟合出不同信誉等级，年利率与客户流失率的关系式。
7 m% c  Q! C1 n$ ~8 n! T( f
- g( z, a0 L* m% U
3 l0 {8 \5 I0 ^' F7 n+ ?$ v2 L信誉等级        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
C        0.9982        y = 32.157x^3-207.386x^2+504.717*x-0.973
银行获利=贷款金额x贷款年利率x(1-利率对于信誉评级客户流失率)


8 [8 f# K: d: K7 r在贷款金额确认，贷款年利率范围在0.4~1.5的情况下，利用上面拟合的关系式，我们能够暴力跑出最优年利率。

$ s9 D( Q+ T6 p6 n3 r) w9 [3 Y. u
double turnover_rate(double x, char ch) {
) ]. l  j6 }4 x) O2 ~8 X& E        double y = 0, result = 0;
9 a, a) A8 T3 u3 n0 b, t; b        switch (ch) {
        case 'A':
; a, e! P- T" `* z                y = 37.969520 * pow(x, 3) - 258.570452 * pow(x, 2) + 640.944427 * x - 1.121484;
                result = x * (1 - y / 100.0);
                break;
3 i  O+ V6 W9 u& j        case 'B':
  A& O" v; s2 h5 h                y = 33.994698 * pow(x, 3) - 225.050538 * pow(x, 2) + 552.829151 * x - 1.016503;
                result = x * (1 - y / 100.0);
                break;
! m. n& d% g9 F        case 'C':
+ {; T; E$ P, R; K                y = 32.156864 * pow(x, 3) - 207.385880 * pow(x, 2) + 504.716993 * x - 0.973497;
                result = x * (1 - y / 100.0);
                break;
        default:
                cout << "输出有误！"  << ch << endl;
        }
! D/ n( s8 n( g' n" {; ~
        return result;
}




问题二
7 d1 W, V# N; h利用代码，重新计算出各指标数据
代入Logistics违约率预测模型，预测出各公司的违约率
2 c% W# `/ z, i( v6 I8 |根据标准普尔评级建立，主标尺，对不同违约率进行A~D等级划分，信誉等级D不予贷款
0 Q3 l9 d( t" L( m- C) k将违约率代入，之前得到的公式，得到具体贷款金额
最优年利率沿用上一问
# 信用等级
% k; x1 C6 H6 R. W& Dcs = []
( [% [9 z  m7 B$ u' F3 X7 g# 最优年利率,客户流失率,利率值
tax = []
7 e0 b8 O4 I8 X* f. h6 ?! T5 K, [for i in m:
: O9 ?0 \6 R0 {) r& }4 H8 G! T    if i <= 0.0069264:
        cs.append('A')
        tax.append([0.083,0.503173,0.0412366])
    elif i > 0.0069264 and i <= 0.22619:
        cs.append('B')
! S. B3 m& @% z  |        tax.append([0.097,0.505215,0.0479942])
    elif i > 0.22619 and i <= 0.509915:
        cs.append('C')
  B1 L) S1 E$ w  V        tax.append([0.1069,0.506501,0.052755])
) c2 l" m( [/ I/ d) q9 F    elif i > 0.509915:
5 b' z6 A2 ]9 j! o+ d2 z        cs.append('D')
        tax.append([0.15,0,0])
# a9 w, U, b- e9 W    else:
3 u& r6 s( K8 \) V& x) K; L        print('违规')
        
parr = []
for arr in list(a):
    parr.append(list(arr)[0])
sum_val = sum(parr)
amount = []
for ival in parr:
    tmp = ival / sum_val * 100000000
    if ival < 1 - 0.509915:
- ?1 t3 n% v8 N8 s& `. K, \            amount.append(0)
    else:
6 Q. A& g. R3 S( a            amount.append(tmp)
. o7 q) r+ |" k4 a& m2 T

! S: x! j6 A- Z0 S) K5 |* B2 I
$ A9 Z! B8 J$ ~3 @) b
& K! O' Z, e& v1 B6 u, z可以看到，贷款金额也都在10w~100w之内。


+ ~& K) r( i. ^: \9 R问题三
/ n) |0 n" g0 c% f0 _这一问，我们做得有些匆忙了，有其他想法的可以按照自己的想法做做，这里只拿我们的做参照。

" y9 I8 i9 L! a; A) T2 `
3 V1 f* ~7 e- B: o' v疫情对公司影响最大的就是每月平均销售额和每月平均销售数量，因此，


对每月平均销售额和每月平均销售数量，分别取随机数，数量取10w组，其他指标数据值不变，是否违约数据根据第二问结果，设为初始值
代入Logistics模型中，预测出每一组的违约率
0 c9 K* y4 W6 a4 |判断每个公司违约率的变化情况，根据变化情况来增/减贷款金额和年利率。
rand_num = pro_rand()

data_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':[]}
- u4 a4 N/ ]" ^1 t7 b* |df = DataFrame(data_form)

  `6 l- _# k( u8 A- i4 x$ pfor val in rand_num:
$ i$ i* b5 Y$ N2 {9 W* i        a = rate_func(val)
& t5 u  G4 ]( E+ F/ N& g        rate_list = get_rate(a)
        df.loc[len(df)]  = rate_list
) P4 O( X# ^! p) k$ S        df.to_excel(writer,sheet_name='违约率变化',index=False)
1 o( l; f( Z) E) }1 I9 L  e, ^: D销售数量和销售金额的随机数范围是：0~MAX
8 [) K$ f( Y0 O$ Q4 A. X+ X



————————————————
& \* z$ K8 b+ Y! a1 U版权声明：本文为CSDN博主「Hk_Mayfly」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_39542714/article/details/108592087

: H' n  r& R1 g7 ]3 b
1 g. h; f$ z6 A  k# H$ _$ c