8.20
昨天比赛出成绩了,一等奖,哈哈哈和队友还是很开心的,希望国赛能好好准备发挥
拿到三个题以后,各自看了半个小时,和队友开了会。因为我们三个都没有偏物理专业的,所以A题就没有考虑,B题C题感觉都可以处理,C题更加开放好写。
整个题目解决比较顺畅,大概花了2天,最后一天对论文进行了润色。比赛完了负责国赛带我们组的老师点评这个题有更好的解法,三个队员都有一些问题,以后继续努力吧
1
所属类别
2021 年“华数杯”全国大学生数学建模竞赛
参赛编号
本科组 CM 211297
基于 GA-BP 模型的客户销售方案
摘要
对电动汽车客户挖掘和销售决策研究是十分有意义的问题,这对于电动汽车行业
快速发展有着重要影响。公司不仅需要提高客户对产品的满意度,对性能做出改变,
更要对汽车的销售,对客户做出更加个性化的销售方案。因此,本文的主要目标就是
对清洗后的数据进行全面分析,从中提炼出最影响客户购买电动汽车的因素,并利用
这些因素建立一套合理有效的针对不同品牌电动汽车的客户挖掘模型。
针对问题一,本文从 a1~a8 指标中找到了四组异常数据进行剔除,并借助箱型图
法进行了检验。借助 B5,B6,B7 以及 B13,B15,B16,B17 之间的逻辑关系对异常数据进行
修改与剔除。得到正常数据后主要针对 a1~a8 指标进行均值、方差、偏度、峰度上的
描述性统计分析。数据处理的最终结果见附件附录一整理.xlsx。
针对问题二,首先对第一问清洗后的客户数据进行训练,以 90%作为训练集,10%
作为数据集进行求解。由于 BP 网络随机生成初始参数易造成局部最优和收敛速度慢,
而后使用遗传算法对其进行优化。最终归一化得到该指标的贡献率,分析即可得出各
因素的影响程度高低,对不同品牌电动汽车的销售影响因素最大的是 a1、B1、B3、B5 和 B9 项。本文将权值贡献率低于 2%的因素划分为影响最小的因素,有 a2、a4、a7、 B4、B6、B11、B14、B15 八项因素。
针对问题三,通过主成分分析法和 GA-BP 神经网络建立客户挖掘模型。利用主成
分分析法对第二问结果中大于 2%的因素集进行降维处理,得到累计贡献率达 85%的客
户评分模型。基于问题一中筛选后的数据进行训练,网络初始化后使用遗传算法优化
神经网络预测,输入需要预测的 15 条数据后得到各位顾客的购买意愿,将 15 名目标
客户购买电动车的可能性填入表格,对模型的优良性进行分析。
针对问题四,本文利用第三问的主成分分析法所得的评分模型对附件三 15 名客户
进行评分,基于客户 ABC 分类管理法筛选出三名未购车的客户。基于 GA-BP 神经网络
的潜在客户挖掘模型进行控制变量法分析,将测试集带入模型进行预测识别,得到各
指标的满意度增长表,见附录满意度增长表.xlsx。用效益难度比来刻画销售者应当在
哪方面提升服务力度,将通过 BP 神经网络得到的结果利用 Matlab 绘图工具绘制折线
图,最终得到在短时间内,提升与 a1 指标相关的服务力度将获得最好的效果。
针对问题五,需要基于上述四问的分析结果和客户生命周期理论给出销售策略建
议。根据问题四应该根据电池技术性能等方面对客户进行销售。结合问题一的描述性
分析与客户挖掘模型得出,针对品牌 1,品牌 2 可以走大众路线,在价格,性价比等方
面可以出台优惠政策;针对品牌 3,则可以走小众化路线,一是可以调查获取数据利用
客户挖掘模型,二是可以在宣传时着重凸显品牌的特点。具体的销售策略见正文。
最后,我们分析了模型的优缺点,在对缺点提出相应的改进方式的基础之上,将
该模型推广到其他领域。
关键词:描述性统计 数据清洗 GA-BP 神经网络 客户 ABC 分类管理 主成分分析
2
-
问题重述 1.1 问题背景
近十几年来,中国的汽车竞争激烈,客户需求快速变化,整车厂必须迅速准确地
认识到自身发展所需的必要条件,才能在高速前进的市场中保持蓬勃发展[1]。 而新能源汽车产业是战略性新兴产业,大力发展以电动汽车为代表的新能源汽车
是解决能源环境问题的有效途径。然而,电动汽车相较于传统汽车,消费者在一些领
域,如电池问题,也存在着一些疑惑。
本文需要根据题目给出的各品牌各项购买意愿来制定出各品牌相应的购买策略,
如电池性能、舒适度、经济性、安全性、动力性等等,同时还需要结合目标客户体验
者个人特征的信息建立不同品牌电动汽车的客户挖掘模型,根据模型判断出目标客户
购买电动车的可能性与销售策略。并对销售部门给出销售策略建议。
1.2 附件分析
本题一共给出了三个附件数据,本文将对这些附件进行初步的解读。
附件一:给出了 1964 位目标客户的体验数据,由三部分组成。第一部分是对各品
牌的各项指标的满意度的打分,包括电池技术性能(电池耐用和充电方便)满意度得
分(满分 100 分,下同)a1、舒适性(环保与空间座椅)整体表现满意度得分 a2、经
济性(耗能与保值率)整体满意度得分 a3、安全性表现(刹车和行车视野)整体满意
度得分 a4、动力性表现(爬坡和加速)整体满意度得分 a5、驾驶操控性表现(转弯和
高速的稳定性)整体满意度得分 a6、外观内饰整体表现满意度得分 a7、配置与质量品
质整体满意度得分 a8。第二部分是目标客户个人特征统计数据,第三部分是各企业的
购买意愿统计。
附件二:给出了调查者针对 1964 位目标客户所做的个人特征调查表,户口、所在
地、驾龄、家庭人数、家庭收入丰方面统计客户数据。
附件三:给出了待判定的数据,问题三需要我们对附件 3 中 15 名目标客户购买电
动车的可能性做出判断。其中 B7 列即客户孩子数部分未知,需要对此数据分析。
1.3 问题要求
现某汽车公司最新推出三款品牌电动汽车,包含合资品牌、自主品牌和新势力品
牌。需要根据消费者对电动汽车的购买意愿制定相应的销售策略。公司统计来了 1964
位目标客户体验三种品牌汽车后给出的满意度得分。此外,题目还提供了目标客户体
验个人特征的信息。
基于上述背景和附件信息我们需要建立数学模型解决以下问题: (1)从 1964 位目标顾客的满意度与特征信息数据统计中提取信息,进行数据清
洗,指出异常值和数据以及处理方法,并且对数据进行描述性统计分析。分析目标客
户对不同品牌汽车满意度的差异比较。 (2)在上述给出的影响客户购买体验的电动车的因素很多,分析电动汽车本身的
因素和目标客户个人特征的因素判断主要影响不同品牌汽车的销售的影响因素。
(3)基于(1)(2)问的数据分析后的结果建立不同品牌电动汽车的客户挖掘模
型并要求模型进行分析。同时分析 15 名目标客户对电动汽车的购买可能性。
(4)综合考虑品牌满意度的提高和服务难度的提高来加大销售力度。基于前面的
研究成果和提高满意度的思路,在挑选的 1 名没有购买电动汽车的客户,给出销售策
略。
(5)基于上述研究结果,向销售部分提出提高销售的销售策划建议。 -
问题分析
3
本文要解决的是基于三种不同品牌的销售策略问题。无论是合资品牌、自主品牌
还是新势力品牌,针对每一品牌,最终都要给出其最合适的销售策略,探究哪些因素可
能会对不同品牌电动汽车的销售有影响并且建立客户挖掘模型。问题一是先要对数据
进行清洗与分析,处理好异常值和缺失值并对目标客户对于不同品牌汽车满意度的比
较分析;而问题二是针对影响销售因素的分析,本文采用 BP 神经网络方法对经过问题
一处理后的数据进行分析,BP 网络求解的具体算法思路如下:
图 1 问题二的解决流程图
问题三是在前两问的基础上建立客户挖掘模型并用该模型判断客户购买电动车的
可能性。针对问题四基于神经网络建立潜在客户挖掘模型。在该模型过程中,需确定
学习样本集,对预测模型进行训练,确定网络参数,进行预测仿真,检验网络;还需测
试样本集,分别进行预测潜在购机人群的综合属性。问题五则是针对上述研究结果对
相关部门提出销售建议,五个问题是一脉相承,层层递进的。 -
模型假设与符号系统
3.1 模型的假设 (1)未婚情况下不会孕育小孩
(2)假设客户在填写调研表时,B7 问题回答不准确以致出现误差数据的可能性
远大于 B6,大于 B5 (3)与父母同住情况时家中的人口数包括父亲、母亲其中一位或两位同居以及与
兄弟姐妹同住的情况
(4)各品牌竞争相互独立
(5)营销者在短时间内只能针对 a1 到 a8 中的某一个指标最多提升 5%的服务满
意度
(6)服务难度和满意度的提升只需考虑为整数的情况,即 1%,2%,3%,4%,5%
3.2 符号系统
本文部分符号说明如下:
符号 意义
因素集 x 矩阵每一列的均值
j
因素集 x 矩阵每一列的标准差
购买电动汽车的客户数
W 初始权值
第 i 项指标提升 n%的效益
Din
第 i 项指标提升 n%的难度
Pin
第 i 项指标提升 n%后客户购车概率
j Ein
4 -
问题一数据处理与分析
4.1 问题一数据的清洗 (1)对 a1~a8 进行异常数据筛选
观察表格数据,发现 a1,a3,a5,a7 列存在异常数据,本文使用 excel 筛选出
a1,a3,a5 超过 100 的数据,筛选出 a7 中过小的数据,将其剔除。
(2)a 项数据使用箱图法筛选数据
使用 matlab 箱图法寻找差距过大的数据,验证 excel 分析正确。Matlab 所得结果
如下图(图中已经标注出了异常的四个数据):
图 2 箱图法筛选 a 项异常数据结果
结果图中展示横坐标为 a1~a8 项,纵坐标为附件 1 中数据,可以看出 a1,a3,
a5,a7 项存在明显的异常数据。
(3)B 项数据的筛选
针对 B1~B17 组数据,本文通过 B5,B6,B7 的关系以及对 B13,B15,B16,B17 的逻 辑关系进行分析筛选并修改出异常数据。
根据附件二可以得知,B5 等项代表的含义如下表:
表 1:B 项主要代码所表示含义
指标代码 含义
B5 家庭人口数
B6 婚姻家庭情况代号
B7 家庭孩子个数
B13 家庭年收入
B15 家庭可支配年收入
B16 全年房贷支出占比
B17 全年车贷支出占比
其中 B6 值为代码,有以下含义:1 表示“未婚,单独居住”,2 表示“未婚,与
父母同住”,3 表示“已婚/同居无子女(两人世界)”,4 表示“已婚/同居无子女
(与父母同住)”,5 表示“已婚,有小孩,不与父母同住”,6 表示“已婚,有小
孩,与父母同住”,7 表示“离异/丧偶”,8 表示“其他”。
5
易知正常数据应该满足:
由于题目中并未给出对年收入更为详细的相关数据,所以本文直接利用 matlab 将
不满足上述不等式的数据剔除。
B5,B6,B7 三组数据应满足某种数学关系,基于假设 1,假设 2,本文利用 matlab
将不缺失的 B7 数据作为准确数据分类讨论,如下: (1)当 B7=1 时,应用 EXCEL 进行数据的筛选可以得到 B6 可取 1,2,3,4,5,6。 若 B6=1,2,3,4 即分别代表“未婚,单独居住”,“未婚,与父母居住”,“已婚/
同居无子女(两人世界)”“已婚/无子女(与父母居住)”。证明是不可能有小孩
的,所以不成立。此时,本文根据 B6,B7 的约束对表中数据进行修改。
若 B6=5,即代表已婚,有小孩,不与父母同住。则可知 B5 为 3,即家庭 3 口
人。据此再次对附件 1 中数据进行筛选。
若 B6=6,即代表已婚,有小孩,与父母同住,则 B5(家庭人口数)可以为
4,5,6,7,8,9。由附件 1 中数据可知当 B7 为 1,B6 为 6 时表格中 B5 只有 3,4,5,6 各
值,可以得知 3 为异常值。本文把 B5 为 3 的数据改为 5,即家庭情况为一对夫妻孕育
有一个孩子,孩子爷爷奶奶一起居住。
(2)当 B7=2 时,应用 EXCEL 进行数据的筛选可以得到 B6 可取 1,2,3, 5,6。 若 B6=1,2,3 即分别代表“未婚,单独居住”,“未婚,与父母居住”,“已婚/
同居无子女(两人世界)”。证明是不可能有小孩的,所以不成立。此时,本文根据
B6,B7 的约束对表中数据进行修改。
若 B6=5,即代表已婚,有小孩,不与父母同住。则可知 B5 为 4,即家庭 4 口
人。据此再次对附件 1 中数据进行筛选。
若 B6=6,即代表已婚,有小孩,与父母同住,则 B5(家庭人口数)可以为
5,6,7,8,9,10。由附件 1 中数据可知当 B7 为 2,B6 为 6 时表格中 B5 只有 4,5,6 各
值,可以得知异常值。本文把 B5 为 4 的数据改为 6,即家庭情况为一对夫妻孕育有两
个孩子,孩子爷爷奶奶一起居住,据此修改附件。
(3)当 B7=3 时,应用 EXCEL 进行数据的筛选 B6 仅可取 5 或 6 且仅只有 6 个
数据。
若 B6=5,即分别代表已婚,有小孩,不与父母同住。所以家庭人数(B5)只能
为 5。此时,本文根据 B6,B7 的约束对表中数据进行修改。
若 B6=6,即代表已婚,有小孩,与父母同住。附件 1 中 2 个数据均合理,此情况
没有异常值。 (4)B7 值为未知时,根据附件 1 可知此时 B6 有近 99%的数据为未婚或没有子
女的情况,由此可知缺失数据应该全为 0。此时将 B6,B7 作为标准数据填入表格,根
据约束条件来调整 B5 的值,修改后的表格程序太长不做赘述,详见支撑文件。
4.2 数据描述性统计分析
本文主要针对满意度进行描述性分析,将产品 1,2,3 独立分析,分别求解品牌 1,
品牌 2,品牌 3 以及 3 种综合的 a1~a8 项数据分析。本文主要从均值、方差、偏度、
峰度四个方面进行描述性统计分析。
4.2.1 均值的分析与评价
均值是表示一组数据集中趋势的量数,它是反映数据集中趋势的一项指标。本文
对 a 项满意度使用 matlab 进行均值处理后得到均值,通过 excel 可得均值如下图:
16 17 13 15+ 13
100
B B B B B
图 3a 项各类满意度均值统计
由上图各均值具体数据可知 3 种品牌综合得分为 77.5 左右,以及各项得分差距并
不大,同时由图可知 a3 与 a5 的得分偏低。a3 反映了经济性(耗能与保值率)整体满
意度得分,a5 为动力性表现(爬坡和加速)整体满意度得分。人们对电动汽车的经济
性与动力性表现满意度不是很高,为提高销售策略,可以从这两方面入手。
4.2.2 方差分析与评价
方差可以衡量数据的离散程度,也可用来度量随机变量和其数学期望(即均值)之
间的偏离程度。本文对 a 项满意度使用 matlab 进行均值处理后得到如下方差,通过 excel
可得下图:
图 4 a 项满意度的各品牌方差值
由图可知,3 种品牌中品牌 3 的方差最大,说明 1964 位体验者对品牌 3 的评价褒贬
不一,但是 3 种品牌中购买意愿的均值品牌 3 最大(为 9.4%),这说明品牌 3 的受众群
体比较专一,若可以为品牌建立良好的客户挖掘模型,则品牌 3 可为公司带来巨大的利
润收益。
4.2.3 偏度的分析与评价
偏度这一指标,又称偏斜系数、偏态系数,是用来帮助判断数据序列的分布规律性
的指标。在偏度系数的绝对值较大的时候,最有可能的含义是“离群”数据离群的程度
很高(很大或很小),亦即分布曲线某侧的拖尾很长。本文通过 matlab 得到满意度的各项
数据的偏度,并绘制如下表格
6
图 5 各品牌满意度的偏度
由图可知 3 种品牌及其总和的偏度分布在 0~-0.7 之间,品牌 1,2 的偏度较为集中,
品牌 3 的偏度较为分散。偏度的分布说明还存在相当多的人对电动汽车存在顾虑,公司
在挖掘客户的同时也应当着力消除人们的顾虑。
4.2.4 峰度的分析与评价
与偏度一样,峰度也是一个用于评价数据系列分布特征的指标。峰度越大,说明该
数据系列中的极端值越多。在正态分布情况下,峰度系数值是 3。峰度的绝对值越大,
说明数据越陡峭,峰度的绝对值大于 3,意味着数据严重不正态。对附件 1 修改后的数
据进行分析,通过 matlab 可得到个品牌的峰度大概情况如下图
图 6 各品牌满意度的峰度
峰度的分布约有 84.3%的数据分布在 3.5~4 之间,比较接近 3,即偏离正态分布的
程度不是很大,也就是说明了数据处理的有效性。 -
问题二的分析与求解
5.1 问题二的分析思路
问题二基于第一问清洗后的附件 1 给出的客户数据进行训练,随机选取数据,以
90%作为训练集,10%作为数据集进行求解。使用 BP 神经网络对导入的数据进行训
练,最终绝对值后归一化就是该指标的贡献率,分析即可得出各因素的影响程度高
低。
使用 MIV 算法即可计算各指标影响占比,BP 网络随机生成初始参数易造成局部最
优和收敛速度慢,故利用遗传算法对其进行优化。
5.2 BP 神经网络的建立与求解
5.2.1 求解 BP 神经网络基本原理 7
8
本文采用的 BP 神经网络由三层组成: 输入层,隐含层与输出层。顾名思义:输入
单元接受外部给的信号与数据;输出单元实现系统处理结果的输出;隐含单元处在输入
和输出单元之间,从网络系统外部是无法观测到隐含单元的结构的。如下图。 图 7 BP 人工神经网络构造示意图
除了上述三个处理信息的单元之外,神经元间的连接强度大小由权值等参数来决定。
在处理时需要先让每个神经元代表对数据进行一次处理: 图 8 神经元处理数据示意图
每个隐含层和输出层神经元输出与输入的函数关系为:
其中 表示神经元 i 与神经元 j 之间连接的权重, 代表神经元 j 的输出, sigmoid
是一个特殊的函数用于将任意实数映射到(0,1)区间。
本文用一个完成训练的神经网络处理回归问题,每个样本拥有 n 个输入。相应地,
神经网络拥有 n 个输入神经元和 1 个输出神经元。
接着将 n 个特征依次送入输入神经元, 隐含层神经元获得输入层的输出并计算自
己输出值, 输出层的神经元根据隐含层输出计算出回归值。
然后开始训练神经网络,首先我们随机初始化连接权重 ,对某一训练样本进行一
次前馈过程得到各神经元的输出。
首先计算输出层的误差: Wij Oj Wij l j 1
( )
1 l I O sigmoid I
e
9
其中 代表神经元 j 的误差, 表示神经元 j 的输出, 表示当前训练样本的参考
输出, sigmod′(x)是上文 sigmod 函数的一阶导数。
计算隐含层误差:
隐含层输出不存在参考值, 使用下一层误差的加权和代替 。计算完误差后
就可以更新 和 j :
其中λ是一个称为学习率的参数,一般在(0,0.1)区间上取值。每一个训练样本都
会更新一次整个网络的参数。再额外设置训练终止的条件。
最简单的训练终止条件为设置最大迭代次数,如将数据集迭代 1000 次后终止训练,
本文此题解答时迭代了八百多次,准确度较高。
单纯的设置最大迭代次数不能保证训练结果的精确度,更好的办法是使用损失函数
作为终止训练的依据。损失函数可以选用输出层各节点的方差:
为了避免神经网络进行无意义的迭代,我们通常在训练数据集中抽出一部分用作校
验。当预测误差高于阈值时提前终止训练。
5.2.2 BP 神经网络的实现
对于 BP 神经网络的构建,本文在对 ANN 进行分析是从其三要素入手,即从网络拓
扑结构、传递函数、学习算法三方面对此题开始研究。通过调整网络的权重和偏置这
两个参数对 BP 网络进行训练分解,即 BP 神经网络的训练过程分两部分:
(1)前向传输,逐层波浪式的传递输出值
正向传播时,输入样本从输入层传入,经各隐层逐层处理后,传向输出层。若输出
层的实际输出与期望的输出不符,则转入误差的反向传播阶段。
(2)逆向反馈,反向逐层调整权重和偏置
反向传播时,将输出以某种形式通过隐层向输入层逐层反传,并将误差分摊给各层
的所有单元,从而获得各层单元的误差信号,此误差信号即作为修正各单元权值的依
据。
至此,通过上述方法完成了一次神经网络的训练过程,通过不断的使用所有数据
记录进行训练,从而得到一个分类模型。每一轮训练都使用数据集的所有记录,遇见
一下两种条件停止:设置最大迭代次数,比如使用数据集迭代 100 次后停止训练,或
者计算训练集在网络上的预测准确率,达到一定门限值后停止训练。 Ej Oj Tj
j j – (T O) Wj
10
图 9 BP 算法的信号流向图
为了提高网络训练的精度,本文采用了一个隐含层,而增加其神经元个数的方法
来获得,这在结构实现上要比增加网络层数简单得多。一般而言,用精度和训练网络
的时间来恒量一个神经网络设计的好坏:
(1)神经元数太少时,网络不能很好的学习,训练迭代的次数也比较多,训练精
度也不高。
(2)神经元数太多时,网络的功能越强大,精确度也更高,训练迭代的次数也
大,可能会出现过拟合现象。
由此,神经网络隐层神经元个数的选取原则是:在能够解决问题的前提下,再加
上一两个神经元,以加快误差下降速度即可。
5.3 GA-BP 神经网络模型
传统的 BP 具有容易陷入局部极小值、学习速度慢、对初始值要求严格等缺点,G
A算法是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模
型,是一种通过模拟自然进化过程搜索最优解的方法[2]。本文将遗传算法应用于 BP
中,利用遗传算法来全局寻优 BP 网络的初始权值与阈值。通过优化神经网络初始权值
和阈值提高学习速度,避免陷入局部极小值,提高模型的速度与准确性。GA-BP 客户预
测模型实施步骤如下:
1)确定 GA-BP 各层神经元个数,在 GA-BP 模型的输入中共 25 个因素,隐含层选
一层,构建网络结构。进行训练集和测试集的划分。
2)归一化处理。由于不同的输入参数具有不同的量纲,不同的样本数据评价标准
不一样,因此需要对其量纲化统一标准,归一化处理也能避免吞噬小值和加速网络收
敛。
3)种群初始化并计算个体适应度。网络的权值和阈值为 GA 中要优化的参数,由
于归一化之后的数据都是介于 -1到 1 之间的小数,如果采用二进制编码容易导致溢
出等问题,故本文采用实数编码,通过编码操作,先随机生成一定数量的种群代表
GA-BP 模型中权值和阈值的集合[4]。 4)迭代求解最佳权值和阈值。迭代即进化,上一代进化到下一代产生新的种群,
计算新种群的适应度,并记录好当前代的最佳染色体个体。当算法终止时,可以获得
全局最优解,此时再通过解码操作得到对应的权值和阈值。
5)神经网络模型训练。利用上面把最优初始阈值权值赋予网络预测。BP 网络训练
后马上对其预测,最后对 GA 和 BP 网络预测结果分析。
11
5.4 模型结果分析
经过八百多次迭代后得到最佳的结果,下面对相关结果进行分析。
下面左图通过均方差来衡量网络的性能,可以看出迭代次数越多性能越好。右图
为内容为记录 Gradient 和 Validation Checks,训练状态的跟踪。
图 10 左图为 Performance 结果,右图为 Training tate 结果
下图左图通过绘制回归线来测量神经网络对应数据的拟合程度。可以看出拟合程
度相对较高。右图为平均适应度,适应度函数的选取会直接影响到本文遗传算法的收
敛速度以及能否找到最优解,因为遗传算法在进化搜索中基本不利用外部信息,仅以
适应度函数为依据,利用种群每个个体的适应度来进行搜索。
图 11 拟合度(左)和适应度(右)
基于遗传算法“优胜劣汰”的特点,遗传算法不需要适应度函数满足连续可微等
条件,唯一要求是针对输入可计算出能加以比较的非负结果。目标函数总取非负值,
并且是以求函数最大值为优化目标,故可直接利用目标函数值作为个体的适应度。 5.4.2 结果分析
通过上述模型可得到三种不同品牌的权重贡献率,如下图所示:
图 12 品牌 1(左)、品牌 2(中)、品牌 3(右)的权重贡献率
从图中可以发现三种品牌的权重贡献率大致相同,数值上仅有细微差别。例如在
a8 指标下,对于品牌 1 的权值贡献率为 5.418%,品牌 2 的权值贡献率为 5.421%,而
品牌 3 的权值贡献率为 5.42%。 对各品牌的销售不同的影响因素可以发现,在权重贡献率大于 2%以上的因素较为
明显,所以取 2%以上的权重为主要影响因素,即 a1, a3, a5, a6, a8, B1, B2, B3,
B5,B7,B8,B9,B10,B12,B13, B16,B17,以上因素会对不同品牌电动汽车的销售有影
响。
5.5 模型的优化分析对比 BP 神经网络通过梯度下降法训练网络,而其初始参数随机生成,若选择不当,会
造成局部最优和收敛速度慢等问题。故采用遗传算法[5]
,对其进行优化,本文的预测都
是使用的 GA-BP 模型。
利用遗传算法优化 BP 神经网络不再简单的随机生成初始参数,而是对种群进行多
次迭代初始化,淘汰较差的参数,最终产生合适的网络初始参数。
使用 GA 优化,主要是对网络进化参数部分调整,优化 BP 神经网络参数的设置,
这样会使预测结果更为准确。本文对纯 BP 模型与 GA-BP 模型进行分析后得到其误差对
比图如下:
图 13 纯 BP 模型与 GA-BP 模型的误差对比
此处由图可以看出,使用遗传算法后的 BP 神经网络模型的误差更小,数值更集中
了,本例表明,用遗传算法优化 BP 神经网络可加快训练收敛速度,提升校准结果精
度。
6. 问题三的建模与求解 6.1 主成分分析模型
6.1.1 模型的建立
当研究的问题涉及到多个具有错综复杂关系的变量时,为简化问题解答,常突出
问题的主要矛盾,抓住主要变量,可使用主成分分析的方法。主成分分析法即是通过
对某事物样本的表象分析,实现揭露使得该事物具有这种表象的主成分,从而实现降
维,简化分析[6]。
针对第二问所得出的柱状图,我们对贡献率再 2%以上的因素进行筛选,得到第三
问的初始因素集:1,3,5,6,8,9,10,11,13,15,16,17,18,20,21,24,
12
13
25。
运用主成分分析法对客户能对公司带来利益的可能性进行分析,从而建立客户挖
掘模型
6.1.2 模型的求解 (1)建立因素集
上述可知有 17 个指标,用字母 p 表示,用 表示购买电动汽车的客户数
由于不同指标的量纲有所不同,为避免量纲错误对客户挖掘模型产生影响,接下
来需要将 x 标准化处理。
, ;
其中 为因素集 x 矩阵每一列的均值, 为因素集 x 矩阵每一列的标准差。 (2)求解相关系数矩阵
构造相关矩阵求解系数矩阵,如下:
(3)求解 R 的特征值及其对应的特征向量
假设 R 矩阵的特征值为: ,特征向量为:
所以由特征向量构成了 17 个新指标:
1 11 1 21 2 171 17
2 12 1 22 2 172 17
17 117 1 217 2 1717 17
y …
y …
…
y …
u x u x u x
u x u x u x
u x u x u x
通过计算累积贡献率寻找累积贡献率达到 85%时的指标 ,选取 为主成
分,求:
q1 q 17
1
( )=
k k k k 累积贡献率
计算附件 3 中的综合得分
11 12 1
1 2
…
… … … …
…
pp
x x x
x
x x x
ij
ij j
j x x
j 1,2 .
1
,. .,17
2 i
,, , j j =
ij p p R (r ) 1
r ,
1
ki kj
ij
k
a a
n
i, 1, 2,…,17 j
1 2 3 4 17 , , , ,…, m m m m m yq
1 11 21 31, 41 171
2 12 22 32, 42 172
17 117 217 317, 417 1717
u [ , , ,…, ]
u [ , , ,…, ]
…
u [ , , ,…, ]
u u u u u
u u u u u
u u u u u
1, 2 q ,…, y y y
14
选取每个品牌得分最高的两名客户为最可能购车的客户。带入数据即可利用
matlab 编程求解各成分贡献值。
6.2 基于 BP 神经网络的潜在客户挖掘模型
由于影响购买电动汽车的因素很多,且影响因素与购买电动车之间通常存在复杂的
非线性关系,所以无法直接求得购买电动汽车与影响因素的函数关系,也就无法通过一
般方法预测函数关系直接得出购买电动汽车的人群特征情况。
所以本文采用误差反传算法或其变化形式的网络模型即 BP 网络,在第二问模型上
进行优化。利用第一问清洗好的数据,处理后的数据划分明显,第二问最后得到的影
响显著的指标也比较符合实际等等,直接利用对附件 3 通过机器学习训练并识别了,
同时因为品牌 1、品牌 2、品牌 3 是单独的,所以将 3 个分别训练,算法流程图如下:
图 14 问题 3 BP 分类模型流程图
6.2.1 样本的选择与预处理
对删除和修改异常值后的数据进行训练,之前处理后的数据(附件 1 中清洗后的
客户的各项指标和满意度)等进行划分。为了更好的进行训练和测试,获得输入和输
出变量,在试验前将试验训练样本集和测试样本集进行数据归一化处理。
6.2.2 网络初始化
根据问题二所建立的的 GA-BP 神经网络模型求解该问,输入层为 25 个节对应 625
影响因素,即 a1,…,a8 各项满意度指标和 b1,…,b17 客户信息指标;输出层为客户的
购买意愿。隐层设计一般为 4-16 个神经元节点,若神经元个数过多,会加大计算量,
导致过度拟合;反之若隐层神经元个数过少,会很大程度上影响该性能。所以本题隐
层设为 8 个神经元节点。
6.2.3 使用遗传算法优化网络预测
带入第二问所建立的 GA-BP 模型进行预测。通过遗传算法优化网络预测后可以得
到训练后的网络,即输入 3 题所需要预测的 15 条数据。得到每位客户的购买意愿率,
对 3 个品牌进行排序,得到每个品牌中意愿率最大的 2 人,本文将这 6 人作为购买电
动车的客户。
6.3 结果分析 表 2: 15 位客户的购买意愿预估值
客户编号 品牌编号 是否购买
1 1 1
2 1 1
3 1 0
4 1 0
5 1 0
6 2 1
7 2 1
8 2 0
9 2 0 q
15
1 1
= ( ) j j j k k
Z y
15
10 2 0
11 3 1
12 3 1
13 3 0
14 3 0
15 3 0
6.4 模型优良性的评价
建立客户挖掘模型,即公司需要针对客户推出更个性化的销售方案。著名学者李
小庆(2018)在金融行业中指出,通过对客户行为的大数据分析,向客户推荐更多定制
化的金融产品以及服务渠道, 体现个性化和差异化服务,从而增强客户的忠诚度和粘
性。于苧(2018)在酒店行业中指出,个性化的服务可以从两个层面理解:(1)在企
业制定的服务标准化的基础上,依照客户自身的需求,为客户提供个性化的服务,达
到使客户满意的结果,从而增加客户的粘性。(2)公司在为客户提供个性化服务的同
时,也在不断改进公司的流程和服务[3]。而本文针对问题三就是根据客户独立的个体进
行个性化的分析。
而本文采用主成分分析可以排除主观因素的干扰,得到客观的结果,而且可以将
附件中相关性强的指标如 B13,B16,B17 等通过主成分分析降维变为相互独立的因素。
累计贡献率的临界值取得为 85%,可以反映原因素集 85%以上的信息,以品牌三为例,
第一主成分中第 11 项指标所占比重最大,说明其在所有因素中的影响最大。
主成分分析后可以发现,以品牌 3 为例,第一主成分中第十一项指标的系数最
大,说明第十一项指标对顾客评分影响最大,与第二问所的品牌三的贡献率柱状图相
照应,说明了模型的准确性。
7. 问题四的分析与求解 7.1 问题分析
题目中要求通过提高服务力度来提升客户满意度增加电动车销量,而在短的时间
内提高 a1-a8 五个百分点的满意度是有可能的,但服务难度与提高的满意度百分点是
成正比的,即提高体验满意度 5%的服务难度是提高体验满意度 1%服务难度的 5 倍。即
是说短时间内最多能对 a1-a8 中某一项提升 5%。
首先通过客户 ABC 分类管理办法选出没有购买电动汽车的客户,基于上述 GA-BP
模型的研究,针对 a1~a8 中第二问求出的显著影响的指标,进行控制变量法分析,可
按 10%、20%、等提高指标值,将其作为测试集带入模型进行预测识别。本问的 GA-BP
模型算法流程图如下:
图 15 问题四 GA-BP 预测算法流程图
7.2 基于客户 ABC 分类管理办法筛选目标客户
客户 ABC 分类法最早由意大利经济学家维尔弗雷多帕累创造。客户 ABC 分类法源自
于帕累托法则,区别在于法则强调关键,而分类法强调主要和次要之间的区别,并将客
户分为 A 类、B 类和 C 类[3]。
帕累托原理是一种定量的研究方法。通常二八定律讨论的是金字塔顶端的 20%的客
16
户,而不是数量最多的 80%位于底端的客户。A 类客户只有 20%,但是给企业贡献了 80%
的收益。不同级别的客户给公司带来的贡献也不一样。企业根据客户的价值,进行差异
化服务策略,合理分配资源,提高企业的竞争力。A、B、C 三类客户分别为关键客户,
约占 10%-20%;重要客户,约占 15%-25%;低价值客户,约占 60-75%。对于企业来说,找
出占比 20%的关键客户,提供更优质的服务,将会带来业绩上的提升和更好的发展[3]。
客户 ABC 分类管理办法是将主要与次要区别开,将客户分为 A 类,B 类和 C 类。二
八定律则是说明公司 80%的利润由 20%的 A 类客户提供。而分类的标准则是根据客户的
价值,此处可以用主成分分析的得分来衡量。下表为根据评分选取的客户信息(红色
为选取客户): 表 4:利用主成分分析法筛选的客户
7.2 依据控制变量法对数据进行 BP 算法处理 (1)对数据进行训练得客户初始购买意愿
基于上述建立的 GA-BP 神经网络模型求解该问,由于 a1,…,a8 指标与购买意愿
有直接关联,所以只需对该问将输入层设置为 8 个节,即 a1,…,a8 各项满意度指标;
隐层设计一般为 4 个神经元节点。于是本文对第 1 问清洗后的数据 a1-a8 列所有数据
进行训练得到训练网。带入上述主成分分析得到的 3 种品牌选择的 3 位客户(编号分
别为 5,10,15)对电动汽车的满意度指标相关数据预测出 3 位客户当前的购买意愿,将
此数据作为初始意愿。
(2)控制变量法分析各项指标
控制变量是在蒙特卡洛方法中用于减少方差的一种技术方法。本题使用该方法
以通过对已知量的了解来减少对未知量估计的误差。
首先控制 a2,…,a8 各项指标数据不变,即仅提高 a1 项电池技术性能(电池耐
用和充电方便)满意度得分,3 位客户 a1 项满意度均增加 1%接得到预测当前购买意
愿,记录;若 3 位客户 a1 项满意度均增加 2%接得到预测当前购买意愿,记录;依次
递增,直到 a1 项满意度增加 5%。 表 5:第一轮计算后客户购买意愿度预测表
客户编号 初始 a1 增长 1% 2% 3% 4% 5%
5 0.2118 0.2214 0.2316 0.2423 0.2537 0.2657
10 -0.0257 -0.0203 -0.015 -0.0097 -0.0047 0.0002
15 0.05 0.0662 0.0833 0.1013 0.12 0.1391
接着将 a1 项数据复原,控制 a1,a3,…,a8 保持不变,将 3 位客户 a2 项满意度
均增加 1%接得到预测当前购买意愿,记录;同理递推,直至 3 位客户满意度均增加
5%接记录得到预测当前购买意愿。
最后得到增长率和购买意愿的数学关系式,选取短时间内提升最高的因素。最
后按该项因素提高购买意愿的方向提升。
17
7.3 预测识别结果及增长率分析
上述每做一次预测,就可以得到三个客户的购买意愿率,控制变量后统计出每增
长 1%服务满意度后的预测值。再使用 matlab 计算出每位客户在不同满意度得分指标
变化下购买意愿的增长速率。
对 BP 模型分析预测后的数据进行处理
为使销售额提高明显,营销者在短时间内应该针对 a1 到 a8 中的某一个指标进行
服务难度的加大,即只能最多提升某一指标的五个百分点。由于题目中说提高体验满
意度 5%的服务难度是提高体验满意度 1%服务难度的 5 倍,为了在销售策略中体现对服
务难度的考量,基于假设(6)我们用效益/服务难度来衡量某一指标提升的效果,此
值越大,说明对此项指标的投入产出比最小。于是我们可以建立如下模型
由于效益可以用 BP 所得的满意度增长率来表示,而服务难度与满意度提升的百分
比成正比:
其中 表示第 i 项指标提升 n 个百分点的效益; 表示第 i 项指标提升 n 个百分
点的难度; 表示第 i 项指标提升 n 个百分点后客户购买该品牌电动汽车的满意度增
长率。
利用 matlab 绘图工具做出以提升的百分比为自变量,以 BP 输出结果为因变量的
折线图,如下图所示:
图 16 各品牌顾客的购买意愿增长情况
上图表示 a 项 8 类满意度随服务难度提升时购买满意度增长率的变化情况。以针对
in i 1,…,8
% 1 . ,. .,5
in
in
E P n D n
, Ein Din Pin
10 号顾客提升五个百分满意度为例,得到的效益难度比如下表所示:
表 5:10 号客户提高体验满意度到 5%后的增益表
由图可以看出基于 a1 指标建立营销策略可以得到最高的效益难度比。不同顾客区
别在于第二选项的不同,例如针对 10 号顾客 a1 与 a3 的差别并不大,所以若有客观因
素的限制,也可以针对 a3 指标建立营销策略。a1 反应的指标是电池技术性能,这与当
前大对数人们对电动汽车的顾虑是相一致的。
7.4 销售策略
据上述结果可知,销售者在向顾客推销产品时,应当侧重对产品电池性能的描
述,着重突出电池性能好的地方,以电池为驱动的汽车作为一项新兴技术成果,公司
应当基于电池性能对大众进行宣传。例如扩大广告投放,并且广告中应尽可能体现电
池驱动车辆给人民生活出行带来的便捷处。同时也应为顾客提供尽可能周到的服务,
基于客户 ABC 分类管理法可知 A 类客户能给公司带来更多的收益,所以在有限条件下
应该首先为 VIP 客人提供更加贴心的服务。
8. 问题五的建议的提出
8.1 分析
本题需基于上述 4 问的分析并结合客户生命周期理论给出销售策略建议。
客户生命周期理论是对时间因素的考量,此理论将客户与公司的关系划分为四个
时期:考察期、形成期、稳定期、退化期。附件一中的数据可认为是在考察期获得的
数据,此阶段仅有少量人购车,对公司的贡献有限;形成期客户开始为公司制造利
益;稳定期客户为公司带来较高利益;退化期客户与公司的交易逐渐减少[3]。 图 17 客户关系生命周期对企业的贡献变化趋势图
第一问本团队从均值、方差、偏度和峰度四个方面分析客户对各种品牌电动汽车
的满意度;第二问分析了需要提高的主要服务指标;第 3 问客户挖掘模型可以帮助公
司筛选出最优潜力购买的 A 类客户;第四问给出了具体需要提高多少的策略。
8.2 提出建议
首先要注重新技术的研发与宣传销售。本团队对各项满意度指标均值进行分析发
现 a3,a5 的得分最低。人们对电动汽车的经济性与动力性表现满意度不是很高,需要
加大相关部分对新技术的研究,做出经济性和动力性能都高的电动汽车,销售时可通
过宣传该品牌电动汽车的性价比来吸引客户等。
其次需要针对不同品牌分别制定销售策略,从各数据的稳定性看,针对品牌 1 和
品牌 2 可以走大众路线,在性价比方面可以出台优惠政策;针对品牌 3,则可以走小众
化路线,一是可以调查获取数据利用客户挖掘模型,二是可以在宣传时着重凸显品牌
18
的特点。
另一方面注重线下体验感销售,从偏度分析可以发现仍有很多人对购买电动汽车
存在顾虑,所以可以从打消客户的购车顾虑出发,例如建立电动汽车体验区,让更多
的人了解到电动汽车的优势,拉进与电动汽车的距离等让电动汽车走进大众的视野。
同时结合客户满意理论,调查客户对电动汽车的期望值,看与实际满意度评分的
差距。对于有购车意愿的人,先让客户做附件 2 中的问卷调查,利用本问所建立的客
户挖掘模型得出结果后,对客户进行 ABC 类划分,提供个性化高质量服务。
最后在向客户推销时,注意侧重对电动汽车电池性能优越性的描述。
9. 模型的评价 9.1 客户挖据模型的评价
基于神经网络的客户挖掘模型具有很多优点,主要有两大突出优点。一是自学习
和自适应能力极强,经过训练,神经网络可自动提取出输入数据与输出数据的“合理
规则”,并将其记忆于权值中,本题中经训练与学习,神经网络高度地学习了已经购买
了电动汽车的客户的各项指标,能够准确的根据客户数据为不同品牌的电动汽车挖掘
出可能会购买的客户。二是神经网络模型具有较强的非映射能力,能够以任意精度逼
近任何非线性的连续函数。
但该法也有着一些缺点。主要是收敛速度慢,本题中主要由学习率决定。如果学
习率过大,将使网络在训练过程中发生振荡,乃至不收敛现象;如果学习率过小,影响
权值和阈值的更新速度从而需要更多的迭代次数,减小了收敛速度。
基于此种情况,本文加入遗传算法以改善 BP 模型的不足,并对两张模型的误差进
行分析。基于此模型选出影响客户决策的主要因素、判断客户购买电动车的潜力并制
定销售策略。
9.2 模型的改进
1.本文训练神经网络需要较长时间,主要是由于学习速率太小所造成的,可采用
变化的或自适应的学习速率来加以改进。可以采用 P 算法改进来加快训练速度,避免
陷入局部极小值等,常见的改进方法有带动量因子算法、自适应学习速率、变化的学
习速率以及作用函数后缩法等。动量因子法的基本思想是在反向传播的基础上,在每
一个权值的变化上加上一项正比于前次权值变化的值,并根据反向传播法来产生新的
权值变化。而自适应学习 速率的方法则是针对一些特定的问题的。改变学习速率的方
法的原则是,若连续几次迭代中,若目标函数对某个权倒数的符号相同,则这个权的
学习速率增加, 反之若符号相反则减小它的学习速率。而作用函数后缩法则是将作用
函数进行平移,即加上一个常数。
2.对于局部最小量方面,可以采用梯度下降法来尽可能收敛到局部最小值,采用
多层网络或较多的神经元,有可能得到更好的结果。
9.3 模型推广
主成分分析模型适用于探究主要影响因素的多变量问题,可将多变量简化为少数
几个变量,简化分析过程,且可明确各因素的影响程度。基于神经网络的客户挖掘模
型在加入遗传算法后适用于数据的规律探索,并对某现象进行预测的问题,其通过其
自学习等能力提出数据的内在规律,从而确定其他发展性内容。
19
20
10. 参考文献
[1] Florian Frederik Deutgen,费蓓娅.中国汽车销售的未来(上)[J].汽车与配
件,2021(06):48-53.
[2] 刘浩然、赵翠香、李轩等.一种基于改进遗传算法的神经网络优化算法研究[J]、仪
器仪表学报,2016,37(7):1573-1580. [3]刘倩. B 公司客户价值管理研究[D].上海外国语大学,2020.
[4]王倚文,许承东,彭雅奇,牛飞.基于 GA-BP 神经网络的 BDS 轨道误差模型研究[J].计
算机仿真,2020,37(02):82-86.
[5] 刘鹏,赵言正,闫维新.基于改进遗传算法的管道机器人摩擦参数辨识[J].中国民
航大学学报,2018,36(06):54-58. [6] 司守奎,孙兆亮.数学建模算法与应用(第 2 版)[M].北京:国防工业出版社, 2015.
(4):231-240.
附录 第一问:数据处理的 matlab 代码
% 数据处理
a=xlsread(‘附录 1 目标客户体验数据.xlsx’,‘C2:J1965’);
ab=xlsread(‘附录 1 目标客户体验数据.xlsx’,‘A2:AB1965’);
boxplot(a);
for i=1:1964
if ab(i,17)==1
if ab(i,15)==3
ab(i,16)=5;
elseif ab(i,15)==4||ab(i,15)==5||ab(i,15)==6
ab(i,16)=6;
elseif ab(i,15)==6
if ab(i,16)==6
ab(i,15)=5;
end
end
elseif ab(i,17)==2
if ab(i,15)==3
ab(i,15)=4;
elseif ab(i,15)==4
ab(i,16)=5;
elseif ab(i,15)==5||ab(i,15)==6
ab(i,16)=6;
end
elseif ab(i,17)==3
if ab(i,15)==3
ab(i,15)=5;
end
elseif ab(i,17)
=1||ab(i,17)
=2||ab(i,17)~=3
if ab(i,16)==5||ab(i,16)==7||ab(i,16)==8
ab(i,:)=0;
end
if ab(i,15)==4
if ab(i,16)==3
ab(i,:)=0;
end
end
if ab(i,1)~=0
ab(i,17)=0;
end
end
end
for i=1:1964
if ab(i,26)+ab(i,27)>100
ab(i,:)=0;
end
if (ab(i,23)
ab(i,26)+ab(i,23)
ab(i,27))/100+ab(i,25)>ab(i,23)
ab(i,:)=0;
end
if ab(i,24)>ab(i,23)
ab(i,:)=0;
end
21
end
xlswrite(‘附录一整理.xlsx’,ab);
% 统计分析
ab1=xlsread(‘附录一整理.xlsx’,‘品牌 1’,‘A2:AB463’);
ab2=xlsread(‘附录一整理.xlsx’,‘品牌 2’,‘A2:AB967’);
ab3=xlsread(‘附录一整理.xlsx’,‘品牌 3’,‘A2:AB107’);
abnew=xlsread(‘附录一整理.xlsx’,‘品牌 123’,‘A2:AB1535’);
mean1=mean(ab1);
mean2=mean(ab2);
mean3=mean(ab3);
meannew=mean(abnew);
xlswrite(‘附录一整理.xlsx’,mean1,2,‘A465:AB465’);
xlswrite(‘附录一整理.xlsx’,mean2,3,‘A969:AB969’);
xlswrite(‘附录一整理.xlsx’,mean3,4,‘A109:AB109’);
xlswrite(‘附录一整理.xlsx’,meannew,5,‘A1537:AB1537’);
var1=var(ab1);
var2=var(ab2);
var3=var(ab3);
varnew=var(abnew);
xlswrite(‘附录一整理.xlsx’,var1,2,‘A467:AB467’);
xlswrite(‘附录一整理.xlsx’,var2,3,‘A971:AB971’);
xlswrite(‘附录一整理.xlsx’,var3,4,‘A111:AB111’);
xlswrite(‘附录一整理.xlsx’,varnew,5,‘A1539:AB1539’);
skewness1=skewness(ab1);
skewness2=skewness(ab2);
skewness3=skewness(ab3);
skewnessnew=skewness(abnew);
xlswrite(‘附录一整理.xlsx’,skewness1,2,‘A469:AB469’);
xlswrite(‘附录一整理.xlsx’,skewness2,3,‘A973:AB973’);
xlswrite(‘附录一整理.xlsx’,skewness3,4,‘A113:AB113’);
xlswrite(‘附录一整理.xlsx’,skewnessnew,5,‘A1541:AB1541’);
kurtosis1=kurtosis(ab1);
kurtosis2=kurtosis(ab2);
kurtosis3=kurtosis(ab3);
kurtosisnew=kurtosis(abnew);
xlswrite(‘附录一整理.xlsx’,kurtosis1,2,‘A471:AB471’);
xlswrite(‘附录一整理.xlsx’,kurtosis2,3,‘A975:AB975’);
xlswrite(‘附录一整理.xlsx’,kurtosis3,4,‘A115:AB115’);
xlswrite(‘附录一整理.xlsx’,kurtosisnew,5,‘A1543:AB1543’);
第二问:基于遗传算法神经网络的预测代码
% 清空环境变量
clc
clear
%
%% 网络结构建立
a=xlsread(‘2.xlsx’);
input_format=a(:,1:25)’;
output=a(:,26)’;
%% 划分训练集和测试集
22
tempp = randperm(size(input_format,2));
P_train = input_format(:,tempp(1:round(0.9
size(input_format,2))));
P_test = input_format(:,tempp(round(0.9
size(input_format,2)):end));
T_train = output(:,tempp(1:round(0.9
size(input_format,2))));
T_test = output(:,tempp(round(0.9
size(input_format,2)):end));
input_train=P_train;
output_train=T_train;
input_test=P_test;
output_test=T_test;
%% 基本参数
inputnum = 25;
hiddennum = 8;
outputnum = 1;
N = size(P_test,2);
%训练样本输入输出数据归一化
[inputn,inputps]=mapminmax(input_train);
[outputn,outputps]=mapminmax(output_train);
%构建网络
net=newff(inputn,outputn,hiddennum);
net.trainParam.epochs = 1000;
net.trainParam.max_fail=1000;
net.trainParam.goal = 1e-6;
net.trainParam.min_grad = 1e-10;
net.trainParam.lr = 0.1;
[net,tr] = train(net,inputn,outputn);
inputn_test=mapminmax(‘apply’,input_test,inputps);
t_sim = sim(net,inputn_test);
t_sim=mapminmax(‘reverse’,t_sim,outputps);
r2 = 1 – (sum((t_sim – output_test).^2) / sum((output_test – mean(output_test)).^2));
%% 遗传算法参数初始化
maxgen=50; %进化代数,即迭代次数
sizepop=50; %种群规模
pcross=[0.4]; %交叉概率选择,0 和 1 之间
pmutation=[0.2]; %变异概率选择,0 和 1 之间
%节点总数
numsum=inputnum
hiddennum+hiddennum+hiddennum
outputnum+outputnum;
lenchrom=ones(1,numsum); %个体长度
bound=[-3
ones(numsum,1) 3
ones(numsum,1)]; %个体范围
%将种群信息定义为一个结构体
individuals=struct(‘fitness’,zeros(1,sizepop), ‘chrom’,[]);
avgfitness=[]; %每一代种群的平均适应度
bestfitness=[]; %每一代种群的最佳适应度
bestchrom=[]; %适应度最好的染色体
%计算个体适应度值
for i=1:sizepop
%随机产生一个种群
individuals.chrom(i,:)=code(lenchrom,bound); %编码(binary 和 grey 的编码结果为一个实
23
数,float 的编码结果为一个实数向量)
x=individuals.chrom(i,:);
%计算适应度
w1=x(1:inputnum
hiddennum);
B1=x(inputnum
hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum
hiddennum+hiddennum+1:inputnum
hiddennum+hiddennum+hiddennum*outputnum)
;
B2=x(inputnum
hiddennum+hiddennum+hiddennum
outputnum+1:inputnum
hiddennum+hiddennum
+hiddennum
outputnum+outputnum);
%网络进化参数
net.trainParam.epochs=20;
net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
net.trainParam.show=100;
net.trainParam.showWindow=0;
%网络权值赋值
net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=reshape(B2,outputnum,1);
%网络训练
net=train(net,inputn,outputn);
an=sim(net,inputn);
error=sum(sum(abs(an-outputn)));
individuals.fitness(i)=error; %染色体的适应度
end
%FitRecord=[];
%找最好的染色体
[bestfitness bestindex]=min(individuals.fitness);
bestchrom=individuals.chrom(bestindex,:); %最好的染色体
avgfitness=sum(individuals.fitness)/sizepop; %染色体的平均适应度
%记录每一代进化中最好的适应度和平均适应度
trace=[avgfitness bestfitness];
%% 迭代求解最佳初始阀值和权值
%进化开始
for i=1:maxgen
disp([num2str(i)])
%选择
individuals=Select(individuals,sizepop);
avgfitness=sum(individuals.fitness)/sizepop;
%交叉
individuals.chrom=cross(pcross,lenchrom,individuals.chrom,sizepop,bound);
%变异
individuals.chrom=Mutation(pmutation,lenchrom,individuals.chrom,sizepop,i,maxgen,bound);
%计算适应度
24
for j=1:sizepop
x=individuals.chrom(j,:); %个体
w1=x(1:inputnum
hiddennum);
B1=x(inputnum
hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum
hiddennum+hiddennum+1:inputnum
hiddennum+hiddennum+hiddennum*outputnum)
;
B2=x(inputnum
hiddennum+hiddennum+hiddennum
outputnum+1:inputnum
hiddennum+hiddennum
+hiddennum
outputnum+outputnum);
%网络进化参数
net.trainParam.epochs=20;
net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
net.trainParam.show=100;
net.trainParam.showWindow=0;
%网络权值赋值
net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=reshape(B2,outputnum,1);
%网络训练
net=train(net,inputn,outputn);
an=sim(net,inputn);
error=sum(sum(abs(an-outputn)));
individuals.fitness(i)=error; %染色体的适应度
individuals.fitness(j)=error;
end
%找到最小和最大适应度的染色体及它们在种群中的位置
[newbestfitness,newbestindex]=min(individuals.fitness);
[worestfitness,worestindex]=max(individuals.fitness);
%最优个体更新
if bestfitness>newbestfitness
bestfitness=newbestfitness;
bestchrom=individuals.chrom(newbestindex,:);
end
individuals.chrom(worestindex,:)=bestchrom;
individuals.fitness(worestindex)=bestfitness;
%记录每一代进化中最好的适应度和平均适应度
avgfitness=sum(individuals.fitness)/sizepop;
trace=[trace;avgfitness bestfitness];
end
%% 遗传算法结果分析
figure(1)
[r c]=size(trace);
plot([1:r]’,trace(:,2),‘b–’);
25
title([‘适应度曲线 ’ ‘终止代数=’ num2str(maxgen)]);
xlabel(‘进化代数’);ylabel(‘适应度’);
legend(‘平均适应度’,‘最佳适应度’);
disp(‘适应度 变量’);
%% 把最优初始阀值权值赋予网络预测
% %用遗传算法优化的 BP 网络进行值预测
w1=x(1:inputnum
hiddennum);
B1=x(inputnum
hiddennum+1:inputnum
hiddennum+hiddennum);
w2=x(inputnum
hiddennum+hiddennum+1:inputnum
hiddennum+hiddennum+hiddennum
outputnum)
;
B2=x(inputnum
hiddennum+hiddennum+hiddennum
outputnum+1:inputnum
hiddennum+hiddennum
+hiddennum
outputnum+outputnum);
net1=newff(inputn,outputn,hiddennum);
net1.iw{1,1}=reshape(w1,hiddennum,inputnum);
net1.lw{2,1}=reshape(w2,outputnum,hiddennum);
net1.b{1}=reshape(B1,hiddennum,1);
net1.b{2}=reshape(B2,outputnum,1);
%% BP 网络训练
%网络进化参数
net1.trainParam.epochs=1000;
%net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
net1.trainParam.max_fail=1000;
net.trainParam.max_fail=1000;
net1.trainParam.goal = 1e-20;
%网络训练
[net1,per2]=train(net1,inputn,outputn);
%% BP 网络预测
%数据归一化
inputn_test=mapminmax(‘apply’,input_test,inputps);
an=sim(net1,inputn_test);
%预测结果反归一化
test_simu=mapminmax(‘reverse’,an,outputps);
r21 = 1 – (sum((test_simu – output_test).^2) / sum((output_test – mean(output_test)).^2));
%% GA 优化 BP 网络预测结果分析
error1=t_sim-output_test
error2=test_simu-output_test
figure(4)
plot(1:N,error1,‘k-
’,1:N,error2,‘r–h’)
legend(‘GA-BPNN’,‘BPNN’);
title(‘两种模型误差对比’);
第二问:权重贡献率分析代码
%% 清空环境变量
clc
clear
load best
26
a=xlsread(‘12.xlsx’) %导入 excell 文件
a=mapminmax(a)
p=a(:,1:25)’ %输入值 a:b,是从 a 行或列到 b 行或列
t=a(:,26)’ %输出值
%最终绝对值后归一化就是你的指标的贡献率
%% 变量筛选 MIV 算法的初步实现(增加或者减少自变量
p=p’;
[m,n]=size§;
yy_temp=p;
% p_increase 为增加 10%的矩阵 p_decrease 为减少 10%的矩阵
for i=1:n
p=yy_temp;
pX=p(:,i);
pa=pX
1.1;
p(:,i)=pa;
aa=[‘p_increase’ int2str(i) ‘=p’];
eval(aa);
end
for i=1:n
p=yy_temp;
pX=p(:,i);
pa=pX
0.9;
p(:,i)=pa;
aa=[‘p_decrease’ int2str(i) ‘=p’];
eval(aa);
end
%% 利用原始数据训练一个正确的神经网络
nntwarn off;
p=p’;
net=net1;
%% 变量筛选 MIV 算法的后续实现(差值计算)
% 转置后 sim
for i=1:n
eval([‘p_increase’,num2str(i),’=transpose(p_increase’,num2str(i),’)’])
end
for i=1:n
eval([‘p_decrease’,num2str(i),’=transpose(p_decrease’,num2str(i),’)’])
end
% result_in 为增加 10%后的输出 result_de 为减少 10%后的输出
for i=1:n
eval([‘result_in’,num2str(i),’=sim(net,’,‘p_increase’,num2str(i),’)’])
end
for i=1:n
eval([‘result_de’,num2str(i),’=sim(net,’,‘p_decrease’,num2str(i),’)’])
end
for i=1:n
eval([‘result_in’,num2str(i),’=transpose(result_in’,num2str(i),’)’])
end
for i=1:n
eval([‘result_de’,num2str(i),’=transpose(result_de’,num2str(i),’)’])
end
%% MIV 的值为各个项网络输出的 MIV 值 MIV 被认为是在神经网络中评价变量相关的最好指标
之一,其符号代表相关的方向,绝对值大小代表影响的相对重要性。
for i=1:n
IV= [‘result_in’,num2str(i), ‘-result_de’,num2str(i)];
27
eval([‘MIV_’,num2str(i) ,’=mean(’,IV,’)’])
end
M(1,:)=abs(MIV_1)
M(2,:)=abs(MIV_2)
M(3,:)=abs(MIV_3)
M(4,:)=abs(MIV_4)
M(5,:)=abs(MIV_5)
M(6,:)=abs(MIV_6)
M(7,:)=abs(MIV_7)
M(8,:)=abs(MIV_8)
M(9,:)=abs(MIV_9)
M(10,:)=abs(MIV_10)
M(11,:)=abs(MIV_11)
M(12,:)=abs(MIV_12)
M(13,:)=abs(MIV_13)
M(14,:)=abs(MIV_14)
M(15,:)=abs(MIV_15)
M(16,:)=abs(MIV_16)
M(17,:)=abs(MIV_17)
M(18,:)=abs(MIV_18)
M(19,:)=abs(MIV_19)
M(20,:)=abs(MIV_20)
M(21,:)=abs(MIV_21)
M(22,:)=abs(MIV_22)
M(23,:)=abs(MIV_23)
M(24,:)=abs(MIV_24)
M(25,:)=abs(MIV_25)
MM=sum(M);
for i=1:25
rate(i)=M(i,:)/MM;
end
%% 作图
figure=bar(rate
100)
xlabel(‘因素集’);
ylabel(‘权值贡献率’);
title(‘MIV 权值贡献率分析’);
第三问:主成分分析代码
%以品牌三为例,品牌 1 品牌 2 仅选取的表格数据不同
AB3=xlsread(‘第三问数据处理.xlsx’,‘品牌 3’,‘C2:S11’);
ab3=zscore(AB3);%数据标准化
std=corrcoef(ab3);%求解相关系数矩阵
[vec,val]=eig(std);%求解特征值,特征向量
newval=diag(val);%将特征值作为新向量
[y,n]=sort(newval);%对新向量排序
rate=y/sum(y);%计算贡献率
sumrate=0;
newi=[];
for k=length(y)👎1
sumrate=sumrate+rate(k);
newi(length(y)+1-k)=n(k);
if sumrate>0.85
28
break;
end
end %记下累积贡献率大 85%的特征值的序号放入 newi 中
fprintf(‘主成分数:%g\n\n’,length(newi));
for i=1:1:length(newi)
for j=1:1:length(y)
aa(i,j)=sqrt(newval(newi(i)))
vec(j,newi(i));
end
end %计算载荷
aaa=aa.
aa;
for i=1:1:length(newi)
for j=1:1:length(y)
zcfzh(i,j)=aa(i,j)/sqrt(sum(aaa(i,:)));
end
end%载荷归一化
zcfz=zcfzh
ab3’;
m=[0.4008 0.1667 0.1320 0.1162 0.0901];%主成分的贡献率
df=m
zcfz;%计算附件一品牌 3 买主得分
xlswrite(‘第三问数据处理.xlsx’,df’,2,‘U2:U11’);
m3=xlsread(‘fulu3.xlsx’,‘C12:S16’);
M3=zscore(m3);
zcfz3=zcfzh
M3’;
df3=m
zcfz3;
xlswrite(‘fulu3.xlsx’,df3’,‘T12:T16’);%计算附件 3 品牌 3 客户得分
第三问:GA-BP 模型预测的代码
%% 该代码为基于遗传算法神经网络的预测代码
% 清空环境变量
clc
clear
%
%% 网络结构建立
a=xlsread(‘2.xlsx’);
input_format=a(:,1:25)’;
output=a(:,26)’;
%% 划分训练集和测试集
tempp = randperm(size(input_format,2));
P_train = input_format(:,tempp(1:round(0.9
size(input_format,2))));
P_test = input_format(:,tempp(round(0.9
size(input_format,2)):end));
T_train = output(:,tempp(1:round(0.9
size(input_format,2))));
T_test = output(:,tempp(round(0.9
size(input_format,2)):end));
input_train=P_train;
output_train=T_train;
input_test=P_test;
output_test=T_test;
29
%% 基本参数
inputnum = 25;
hiddennum = 8;
outputnum = 1;
N = size(P_test,2);
%训练样本输入输出数据归一化
[inputn,inputps]=mapminmax(input_train);
[outputn,outputps]=mapminmax(output_train);
%构建网络
net=newff(inputn,outputn,hiddennum);
net.trainParam.epochs = 1000;
net.trainParam.max_fail=1000;
net.trainParam.goal = 1e-6;
net.trainParam.min_grad = 1e-10;
net.trainParam.lr = 0.1;
[net,tr] = train(net,inputn,outputn);
inputn_test=mapminmax(‘apply’,input_test,inputps);
t_sim = sim(net,inputn_test);
t_sim=mapminmax(‘reverse’,t_sim,outputps);
r2 = 1 – (sum((t_sim – output_test).^2) / sum((output_test –
mean(output_test)).^2));
%% 遗传算法参数初始化
maxgen=50; %进化代数,即迭代次数
sizepop=50; %种群规模
pcross=[0.4]; %交叉概率选择,0 和 1 之间
pmutation=[0.2]; %变异概率选择,0 和 1 之间
%节点总数
numsum=inputnum
hiddennum+hiddennum+hiddennum
outputnum+outputnum;
lenchrom=ones(1,numsum); %个体长度
bound=[-3
ones(numsum,1) 3
ones(numsum,1)]; %个体范围
%将种群信息定义为一个结构体
individuals=struct(‘fitness’,zeros(1,sizepop), ‘chrom’,[]);
avgfitness=[]; %每一代种群的平均适应度
bestfitness=[]; %每一代种群的最佳适应度
bestchrom=[]; %适应度最好的染色体
%计算个体适应度值
for i=1:sizepop
30
%随机产生一个种群
individuals.chrom(i,:)=code(lenchrom,bound); %编码(binary 和 grey 的
编码结果为一个实数,float 的编码结果为一个实数向量)
x=individuals.chrom(i,:);
%计算适应度
w1=x(1:inputnum
hiddennum);
B1=x(inputnum
hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum
hiddennum+hiddennum+1:inputnum
hiddennum+hiddennum+hiddennu
m*outputnum);
B2=x(inputnum
hiddennum+hiddennum+hiddennum
outputnum+1:inputnum
hiddennu
m+hiddennum+hiddennum
outputnum+outputnum);
%网络进化参数
net.trainParam.epochs=20;
net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
net.trainParam.show=100;
net.trainParam.showWindow=0;
%网络权值赋值
net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=reshape(B2,outputnum,1);
%网络训练
net=train(net,inputn,outputn);
an=sim(net,inputn);
error=sum(sum(abs(an-outputn)));
individuals.fitness(i)=error; %染色体的适应度
end
%FitRecord=[];
%找最好的染色体
[bestfitness bestindex]=min(individuals.fitness);
bestchrom=individuals.chrom(bestindex,:); %最好的染色体
avgfitness=sum(individuals.fitness)/sizepop; %染色体的平均适应度
%记录每一代进化中最好的适应度和平均适应度
trace=[avgfitness bestfitness];
31
%% 迭代求解最佳初始阀值和权值
%进化开始
for i=1:maxgen
disp([num2str(i)])
%选择
individuals=Select(individuals,sizepop);
avgfitness=sum(individuals.fitness)/sizepop;
%交叉
individuals.chrom=cross(pcross,lenchrom,individuals.chrom,sizepop,bound);
%变异
individuals.chrom=Mutation(pmutation,lenchrom,individuals.chrom,sizepop,i
,maxgen,bound);
%计算适应度
for j=1:sizepop
x=individuals.chrom(j,:); %个体
w1=x(1:inputnum
hiddennum);
B1=x(inputnum
hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum
hiddennum+hiddennum+1:inputnum
hiddennum+hiddennum+hiddennu
m*outputnum);
B2=x(inputnum
hiddennum+hiddennum+hiddennum
outputnum+1:inputnum
hiddennu
m+hiddennum+hiddennum
outputnum+outputnum);
%网络进化参数
net.trainParam.epochs=20;
net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
net.trainParam.show=100;
net.trainParam.showWindow=0;
%网络权值赋值
net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=reshape(B2,outputnum,1);
%网络训练
net=train(net,inputn,outputn);
32
an=sim(net,inputn);
error=sum(sum(abs(an-outputn)));
individuals.fitness(i)=error; %染色体的适应度
individuals.fitness(j)=error;
end
%找到最小和最大适应度的染色体及它们在种群中的位置
[newbestfitness,newbestindex]=min(individuals.fitness);
[worestfitness,worestindex]=max(individuals.fitness);
%最优个体更新
if bestfitness>newbestfitness
bestfitness=newbestfitness;
bestchrom=individuals.chrom(newbestindex,:);
end
individuals.chrom(worestindex,:)=bestchrom;
individuals.fitness(worestindex)=bestfitness;
%记录每一代进化中最好的适应度和平均适应度
avgfitness=sum(individuals.fitness)/sizepop;
trace=[trace;avgfitness bestfitness];
end
%% 遗传算法结果分析
figure(1)
[r c]=size(trace);
plot([1:r]’,trace(:,2),‘b–’);
title([‘适应度曲线 ’ ‘终止代数=’ num2str(maxgen)]);
xlabel(‘进化代数’);ylabel(‘适应度’);
legend(‘平均适应度’,‘最佳适应度’);
disp(‘适应度 变量’);
%% 把最优初始阀值权值赋予网络预测
% %用遗传算法优化的 BP 网络进行值预测
w1=x(1:inputnum
hiddennum);
B1=x(inputnum
hiddennum+1:inputnum
hiddennum+hiddennum);
w2=x(inputnum
hiddennum+hiddennum+1:inputnum
hiddennum+hiddennum+hiddennu
m
outputnum);
B2=x(inputnum
hiddennum+hiddennum+hiddennum
outputnum+1:inputnum
hiddennu
m+hiddennum+hiddennum
outputnum+outputnum);
net1=newff(inputn,outputn,hiddennum);
net1.iw{1,1}=reshape(w1,hiddennum,inputnum);
33
34
net1.lw{2,1}=reshape(w2,outputnum,hiddennum);
net1.b{1}=reshape(B1,hiddennum,1);
net1.b{2}=reshape(B2,outputnum,1);
%% BP 网络训练
%网络进化参数
net1.trainParam.epochs=1000;
%net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
net1.trainParam.max_fail=1000;
net.trainParam.max_fail=1000;
net1.trainParam.goal = 1e-20;
%网络训练
[net1,per2]=train(net1,inputn,outputn);
%% BP 网络预测
%数据归一化
inputn_test=mapminmax(‘apply’,input_test,inputps);
an=sim(net1,inputn_test);
%预测结果反归一化
test_simu=mapminmax(‘reverse’,an,outputps);
r21 = 1 – (sum((test_simu – output_test).^2) / sum((output_test –
mean(output_test)).^2));
%% GA 优化 BP 网络预测结果分析
error1=t_sim-output_test
error2=test_simu-output_test
figure(4)
plot(1:N,error1,‘k-
’,1:N,error2,‘r–h’)
legend(‘GA-BPNN’,‘BPNN’);
title(‘两种模型误差对比’);
aa=xlsread(‘31.xlsx’);
aa=mapminmax(‘apply’,aa’,inputps);
aa_output=sim(net,aa)
aa_output=mapminmax(‘reverse’,aa_output,outputps)
第四问:效益增度代码
huatu5.m
%第五名客户
x=0:1:5;
a1=[0.2118 0.2214 0.2316 0.2423 0.2537 0.2657];
a2=[0.2118 0.2149 0.218 0.2211 0.2242 0.2274];
a3=[0.2118 0.2178 0.224 0.2304 0.237 0.2438];
35
a4=[0.2118 0.2199 0.2281 0.2362 0.2442 0.252];
a5=[0.2118 0.2095 0.2071 0.2046 0.2019 0.1992];
a6=[0.2118 0.1982 0.1856 0.1738 0.1629 0.1529];
a7=[0.2118 0.2131 0.2145 0.2159 0.2174 0.2189];
a8=[0.2118 0.1968 0.1831 0.1705 0.1591 0.1487];
plot(x,a1,’-ro’,x,a2,’-go’,x,a3,’-bo’,x,a4,’-yo’,x,a5,’-r
’,x,a6,’-
g*’,x,a7,’-b*’,x,a8,’-y*’);
set(gca,‘XTick’,[0:1:6])
legend(‘a1’,‘a2’,‘a3’,‘a4’,‘a5’,‘a6’,‘a7’,‘a8’);
xlabel(‘百分点’);
ylabel(‘购买概率’);
huatu10.m
%第十名客户
x=0:1:5;
a1=[-0.0257 -0.0203 -0.015 -0.0097 -0.0047 0.0002];
a2=[-0.0257 -0.029 -0.0322 -0.0354 -0.0384 -0.0412];
a3=[-0.0257 -0.0207 -0.0158 -0.0108 -0.0059 -0.0011];
a4=[-0.0257 -0.0311 -0.0359 -0.04 -0.0435 -0.0463];
a5=[-0.0257 -0.0218 -0.0179 -0.0138 -0.0096 -0.0053];
a6=[-0.0257 -0.0287 -0.0312 -0.0331 -0.0343 -0.0349];
a7=[-0.0257 -0.0224 -0.0191 -0.0158 -0.0125 -0.0092];
a8=[-0.0257 -0.026 -0.0255 -0.0242 -0.0219 -0.0189];
plot(x,a1,’-ro’,x,a2,’-go’,x,a3,’-bo’,x,a4,’-yo’,x,a5,’-r*’,x,a6,’-
g*’,x,a7,’-b*’,x,a8,’-y*’);
set(gca,‘XTick’,[0:1:6])
legend(‘a1’,‘a2’,‘a3’,‘a4’,‘a5’,‘a6’,‘a7’,‘a8’);
xlabel(‘百分点’);
ylabel(‘购买概率’);
huatu15.m
%第十五名客户
x=0:1:5;
a1=[0.05 0.0662 0.0833 0.1013 0.12 0.1391];
a2=[0.05 0.0542 0.0586 0.0633 0.0681 0.0731];
a3=[0.05 0.0541 0.0582 0.0625 0.0669 0.0713];
a4=[0.05 0.0561 0.0623 0.0686 0.0748 0.0808];
a5=[0.05 0.0651 0.0813 0.0981 0.1154 0.1327];
36
a6=[0.05 0.0486 0.0473 0.046 0.0449 0.044];
a7=[0.05 0.0529 0.056 0.059 0.0621 0.0652];
a8=[0.05 0.0432 0.0361 0.029 0.022 0.0151];
plot(x,a1,’-ro’,x,a2,’-go’,x,a3,’-bo’,x,a4,’-yo’,x,a5,’-r*’,x,a6,’-
g*’,x,a7,’-b*’,x,a8,’-y*’);
set(gca,‘XTick’,[0:1:6])
legend(‘a1’,‘a2’,‘a3’,‘a4’,‘a5’,‘a6’,‘a7’,‘a8’);
xlabel(‘百分点’);
ylabel(‘购买概率’);
第 4 问:基于控制变量法的 BP 预测代码
%% 该代码为基于遗传算法神经网络的预测代码
% 清空环境变量
clc
clear
%
%% 网络结构建立
a=xlsread(‘12.xlsx’);
input_format=a(:,1:8)’;
output=a(:,26)’;
%% 划分训练集和测试集
tempp = randperm(size(input_format,2));
P_train = input_format(:,tempp(1:round(0.9
size(input_format,2))));
P_test = input_format(:,tempp(round(0.9
size(input_format,2)):end));
T_train = output(:,tempp(1:round(0.9
size(input_format,2))));
T_test = output(:,tempp(round(0.9
size(input_format,2)):end));
input_train=P_train;
output_train=T_train;
input_test=P_test;
output_test=T_test;
%% 基本参数
inputnum = 8;
hiddennum = 4;
outputnum = 1;
N = size(P_test,2);
%训练样本输入输出数据归一化
[inputn,inputps]=mapminmax(input_train);
[outputn,outputps]=mapminmax(output_train);
%构建网络
net=newff(inputn,outputn,hiddennum);
net.trainParam.epochs = 1000;
net.trainParam.max_fail=1000;
net.trainParam.goal = 1e-6;
net.trainParam.min_grad = 1e-10;
net.trainParam.lr = 0.1;
[net,tr] = train(net,inputn,outputn);
inputn_test=mapminmax(‘apply’,input_test,inputps);
t_sim = sim(net,inputn_test);
t_sim=mapminmax(‘reverse’,t_sim,outputps);
r2 = 1 – (sum((t_sim – output_test).^2) / sum((output_test –
mean(output_test)).^2));
%% 遗传算法参数初始化
maxgen=50; %进化代数,即迭代次数
sizepop=50; %种群规模
pcross=[0.4]; %交叉概率选择,0 和 1 之间
pmutation=[0.2]; %变异概率选择,0 和 1 之间
%节点总数
numsum=inputnum
hiddennum+hiddennum+hiddennum
outputnum+outputnum;
lenchrom=ones(1,numsum); %个体长度
bound=[-3
ones(numsum,1) 3
ones(numsum,1)]; %个体范围
%将种群信息定义为一个结构体
individuals=struct(‘fitness’,zeros(1,sizepop), ‘chrom’,[]);
avgfitness=[]; %每一代种群的平均适应度
bestfitness=[]; %每一代种群的最佳适应度
bestchrom=[]; %适应度最好的染色体
%计算个体适应度值
for i=1:sizepop
%随机产生一个种群
individuals.chrom(i,:)=code(lenchrom,bound); %编码(binary 和 grey 的
编码结果为一个实数,float 的编码结果为一个实数向量)
x=individuals.chrom(i,:);
%计算适应度
w1=x(1:inputnum
hiddennum);
B1=x(inputnum
hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum
hiddennum+hiddennum+1:inputnum
hiddennum+hiddennum+hiddennu
m*outputnum);
B2=x(inputnum
hiddennum+hiddennum+hiddennum
outputnum+1:inputnum
hiddennu
37
m+hiddennum+hiddennum
outputnum+outputnum);
%网络进化参数
net.trainParam.epochs=20;
net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
net.trainParam.show=100;
net.trainParam.showWindow=0;
%网络权值赋值
net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=reshape(B2,outputnum,1);
%网络训练
net=train(net,inputn,outputn);
an=sim(net,inputn);
error=sum(sum(abs(an-outputn)));
individuals.fitness(i)=error; %染色体的适应度
end
%FitRecord=[];
%找最好的染色体
[bestfitness bestindex]=min(individuals.fitness);
bestchrom=individuals.chrom(bestindex,:); %最好的染色体
avgfitness=sum(individuals.fitness)/sizepop; %染色体的平均适应度
%记录每一代进化中最好的适应度和平均适应度
trace=[avgfitness bestfitness];
%% 迭代求解最佳初始阀值和权值
%进化开始
for i=1:maxgen
disp([num2str(i)])
%选择
individuals=Select(individuals,sizepop);
avgfitness=sum(individuals.fitness)/sizepop;
%交叉
individuals.chrom=cross(pcross,lenchrom,individuals.chrom,sizepop,bound);
%变异
38
individuals.chrom=Mutation(pmutation,lenchrom,individuals.chrom,sizepop,i
,maxgen,bound);
%计算适应度
for j=1:sizepop
x=individuals.chrom(j,:); %个体
w1=x(1:inputnum
hiddennum);
B1=x(inputnum
hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum
hiddennum+hiddennum+1:inputnum
hiddennum+hiddennum+hiddennu
m*outputnum);
B2=x(inputnum
hiddennum+hiddennum+hiddennum
outputnum+1:inputnum
hiddennu
m+hiddennum+hiddennum
outputnum+outputnum);
%网络进化参数
net.trainParam.epochs=20;
net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
net.trainParam.show=100;
net.trainParam.showWindow=0;
%网络权值赋值
net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=reshape(B2,outputnum,1);
%网络训练
net=train(net,inputn,outputn);
an=sim(net,inputn);
error=sum(sum(abs(an-outputn)));
individuals.fitness(i)=error; %染色体的适应度
individuals.fitness(j)=error;
end
%找到最小和最大适应度的染色体及它们在种群中的位置
[newbestfitness,newbestindex]=min(individuals.fitness);
[worestfitness,worestindex]=max(individuals.fitness);
%最优个体更新
39
if bestfitness>newbestfitness
bestfitness=newbestfitness;
bestchrom=individuals.chrom(newbestindex,:);
end
individuals.chrom(worestindex,:)=bestchrom;
individuals.fitness(worestindex)=bestfitness;
%记录每一代进化中最好的适应度和平均适应度
avgfitness=sum(individuals.fitness)/sizepop;
trace=[trace;avgfitness bestfitness];
end
%% 遗传算法结果分析
figure(1)
[r c]=size(trace);
plot([1:r]’,trace(:,2),‘b–’);
title([‘适应度曲线 ’ ‘终止代数=’ num2str(maxgen)]);
xlabel(‘进化代数’);ylabel(‘适应度’);
legend(‘平均适应度’,‘最佳适应度’);
disp(‘适应度 变量’);
%% 把最优初始阀值权值赋予网络预测
% %用遗传算法优化的 BP 网络进行值预测
w1=x(1:inputnum
hiddennum);
B1=x(inputnum
hiddennum+1:inputnum
hiddennum+hiddennum);
w2=x(inputnum
hiddennum+hiddennum+1:inputnum
hiddennum+hiddennum+hiddennu
m
outputnum);
B2=x(inputnum
hiddennum+hiddennum+hiddennum
outputnum+1:inputnum
hiddennu
m+hiddennum+hiddennum
outputnum+outputnum);
net1=newff(inputn,outputn,hiddennum);
net1.iw{1,1}=reshape(w1,hiddennum,inputnum);
net1.lw{2,1}=reshape(w2,outputnum,hiddennum);
net1.b{1}=reshape(B1,hiddennum,1);
net1.b{2}=reshape(B2,outputnum,1);
%% BP 网络训练
%网络进化参数
net1.trainParam.epochs=1000;
%net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
net1.trainParam.max_fail=1000;
net.trainParam.max_fail=1000;
net1.trainParam.goal = 1e-20;
40
%网络训练
[net1,per2]=train(net1,inputn,outputn);
%% BP 网络预测
%数据归一化
inputn_test=mapminmax(‘apply’,input_test,inputps);
an=sim(net1,inputn_test);
%预测结果反归一化
test_simu=mapminmax(‘reverse’,an,outputps);
r21 = 1 – (sum((test_simu – output_test).^2) / sum((output_test –
mean(output_test)).^2));
%% GA 优化 BP 网络预测结果分析
error1=t_sim-output_test
error2=test_simu-output_test
figure(4)
plot(1:N,error1,‘k-*’,1:N,error2,‘r–h’)
legend(‘GA-BPNN’,‘BPNN’);
title(‘两种模型误差对比’);
aa=xlsread(‘41.xlsx’);
aa=mapminmax(‘apply’,aa’,inputps);
aa_output=sim(net,aa)
aa_output=mapminmax(‘reverse’,aa_output,outputps)
41