2020华数杯数学建模C题脱贫帮扶绩效评价求解全过程文档及程序

2020华数杯数学建模

C题 脱贫帮扶绩效评价

原题再现：

  消除贫困、改善民生、逐步实现共同富裕，是社会主义的本质要求，是我们党的重要使命。党的十八大以来，国家把扶贫开发工作纳入“四个全面”战略布局，作为实现第一个百年奋斗目标的重点工作，摆在更加突出的位置。
  为了更好的激励各帮扶单位提高扶贫效率，扶真贫，真扶贫。五年前，国家启动了脱贫帮扶绩效评价机制。某科研团队接受任务后，对全国 32165 个需要帮扶的贫困村进行了初步的贫困调查。从居民收入(记为 SR)、产业发展(记为 CY)、居住环境(记为 HJ)、文化教育(记为 WJ)、基础设施(记为 SS)等五个评价指标给出了评分。以此为依据，将被帮扶的村庄划分为 160 个集合，每个集合指定帮扶单位（标记为 0-159）进行帮扶。这 160 个帮扶单位按照单位属性（如国企还是民营企业等）标记为 0-5 等 6 个类型。
  2020 年，研究团队再次进行了调研，得到了被帮扶的这些村庄居民收入、产业发展、居住环境、文化教育、基础设施等五个方面的评分数据以及总分数据。为了便于比较和研究，所有数据都进行了标准化处理（标准化后的数值越大表示评分越高）。绩效评价，不能仅以最后的得分作为依据，需要考虑各个评价指标的进步幅度。因为这样才能公正的评判帮扶的效果，就能鼓励更多的帮扶单位愿意花精力去帮助非常贫困的地区。
  请你运用数学建模，解决以下问题：
  问题一：一般人的理解是，五年前的评分与 2020 年对应的各项评分有着直接的关联，如五年前的居民收入不错，现在的居民收入也会是不错的。你认为本问题有这个规律吗？请分析各个评价指标的对应关系。
  问题二：160 个帮扶单位帮扶着基础不同的村庄，帮扶单位帮扶工作的态度、目标、投入、帮扶干部素质等显然是有差异的。仅仅用 2020 年各村庄评分高低显然是无法真正有效的体现一个帮扶单位在脱贫攻坚提升方面所做出的努力。请你运用附件的数据，阐明什么类型的帮扶单位，哪些帮扶单位在脱贫帮扶上面有较高的绩效？请给不同类型的帮扶单位绩效排序，给出脱贫帮扶绩效前十名的帮扶单位编号。
  问题三：每个帮扶单位在扶贫上有不同的工作特色，如有些单位在提高居民收入上效果很好，而有些帮扶单位可能在改善基础设施上帮助的效果不错。请问，哪些帮扶单位分别在居民收入、产业发展、居住环境、文化教育、基础设施等评价指标上帮扶业绩明显？请列出各单项评价指标前五名的帮扶单位编号。
  问题四：全国计划给予 10000 个村庄“脱贫先进村庄”称号。请问，哪些因素对获得这个荣誉称号有着非常重要的影响？数据表中最后有 10 个村庄的 2020年的评价分数被删除，请你判断他们能评上“脱贫先进村庄”称号吗？如果称号分为一级和二级（一、二级称号比例为 1:3），这 10 个村庄中谁能评上“脱贫先进一级村庄”称号？
  问题五：依据你的研究成果，向国家扶贫办写一封 500 字左右的信，阐述你的观点和建议

整体求解过程概述(摘要)

  2020 年是脱贫攻坚战的决胜之年，全面建成小康社会，是中国共产党对中国人民的庄严承诺。自 2015 年以来，我国脱贫攻坚工作取得了卓越的成绩，各级扶贫单位在脱贫攻坚岗位上取得了不可磨灭的政绩。然而，在评价扶贫情况时，评价标准不够公平，评判指标不够综合等问题接踵而至，这些问题最终将影响到脱贫攻坚战的部署和扶贫工作的质量和效率。基于以上问题，本文从数据分析的角度入手，利用统计学方法，对相关扶贫单位和村庄在不同年份的数据和指标进行建模，对相关问题提出了如下解决方案：
  针对问题一，本文将 2015 年和 2020 年某一村庄的某一单项指标数据组合成为一个二维点，通过对每一个村庄的扶贫情况进行如上处理，得到一个二维点图。利用最小二乘拟合，得到了各个指标的线性方程，通过观察系数则可以直接得到该指标数据在 2015年和 2020 年的相关性情况，发现相关性均大于 0.5，进而判断出 2015 年的 2020 年对应评分存在相关性，符合现实发展规律。不仅如此，本文利用上述方法，对五项指标分别拟合得到的方程两两之间进行相关性分析，总结出各个指标之内，均存在一定关联和规律的结论。
  针对问题二，本文考虑对某一帮扶单位及某一类型的所有帮扶单位，在两年的五项评价指标分别进行特征提取，利用主成分分析（PCA）方法进行数据投影，得到其投影后的降维数据和因子得分情况，并对两组降维后的一维数据作差，分别拟合出不同的帮扶单位个体及不同类型的帮扶单位在总体绩效上的进步幅度，进而对进步幅度和评分进行加权求和，得到其最终的绩效排序。
  针对问题三，本文对各帮扶单位中，各个指标组合过后的二维数据进行 K-means 聚类分析，得到其对应指标的聚类中心点，并对应其坐标，得到该单位在 2015 年及 2020年该指标数据的新的投影坐标，利用对两坐标值作差的方法得到其当前指标的帮扶业绩，并对其进行排序，进而分析出各帮扶单位排名。
  针对问题四，本文考虑构建神经网络进行数据预测，对神经网络进行有监督训练。由于该数据集中的数据均为浮点数，数据精密，且标准化之后的数据差距较小，因此本文采用高斯-伯努利受限玻尔兹曼机（GBRBM）对一个三层神经网络进行分层初始化，再采用反向传播算法对网络进行微调，得到一个训练好的神经网络。基于以上原理，本文采用各指标数据和问题二得到的各村庄绩效训练新的神经网络，利用神经网络中的两个权值矩阵，将二者相乘，得到所有指标对绩效评分的影响权重。同时，作为一种数据预测模型，我们将最后 10 个村庄的数据输入神经网络中，将输出的数据作为其 2020 年的五项指标，并利用之前标准化的矩阵迹作为各指标权值，对所有指标进行加权求和，进而得到各村庄对应总分，混入所有村庄中依次排列，按照 1：3 的比例选取一级村庄。
  最后，针对问题五，我们总结了对五项指标的分析结论，对脱贫攻坚和扶贫工作提出建议，并向国家扶贫办写了一封信，希望以数据科学的力量，助力脱贫攻坚。

问题分析：

  本题旨在利用数学建模的方法，基于各村庄在 2015 年和 2020 年五个指标数据（居民收入、产业发展、居民环境、文化教育、基础设施），以及评估总分，对脱贫攻坚过程中的扶贫工作评估和绩效考核问题提出考核方案。
  问题一要求对单一指标寻找其中是否存在时间层面上的关联和规律，并分析五个评价指标之间是否存在对应关系。
  问题二要求对脱贫帮扶的帮扶单位和不同类型的所有扶贫单位的绩效进行评估，并对不同单位绩效进行评估排序。
  问题三要求从单一指标的角度，利用有效的数学模型，对各个指标层面上表现最佳的单位进行评估，并将对应的帮扶单位列举排序。
  问题四要求对影响“脱贫先进村庄”荣誉评选的因素进行评估，并通过数学建模，对一些数据被删除的贫困村庄的数据进行评估，以判定其是否能够获得此项荣誉。
  问题五则要求利用以上问题得到的研究结果，向国家扶贫办以书信的形式阐述自己的观点和建议，向国家提出自己的扶贫帮扶策略。

模型假设：

  (1) 所有标准化之后的数据之间仍可直接进行运算，即不受原始数据分布的均值和标准差影响；
  (2) 假定村庄脱贫评价仅考虑居民收入、产业发展、居住环境、文化教育和基础设施五项指标，且相对于任何其他因素独立；
  (3) 假定数据集内所有数据均真实、有效。具有统计学意义。

模型的建立与求解

单指标成对数据相关性判断

  该部分所建模型，旨在探索 2015 年与 2020 年的某项指标是否存在某种关联和规律，以及各项评价指标之间是否存在对应关系。从统计学的角度分析，考察两项指标的关联情况可以考虑采用相关性分析[1]，对这些指标进行检验，以确认其间是否存在规律。因此，如何对同一指标在两年之间进行相关性分析成为了一项亟待解决的问题。

  然而，仅仅通过画图观察其数据分布情况并不能观察出某一指标的关联和规律，更不能获得其相关性。为了获取各项指标的深层规律，模型必须对其分布的性质进行更加深刻的分析。

  基于以上定义，我们可以通过拟合某一指标分布的线性函数，观察其系数，来确定某一指标在 2015 年和 2020 年的数据之间的相关性和规律。

探究单指标评分关联性及其规律

  通过对 2015 年和 2020 年的各项指标分别进行组合，本文得到了各年各指标的二维点分布。同时利用最小二乘法，对各个分布进行线性方程拟合，得到如图 1 所示的分布和直线方程。

  可见，所有方程的相关系数均大于 0，这也就意味着，对于任意一个指标，其 2015 年的数据与 2020 年的数据呈正相关，该组数据中存在题干所述规律。

探究各指标之间对应关系

  接下来，本文将通过对数据进行相关性分析，探索五个指标之内是否存在对应关系。通过对五个指标在各年的数据进行相关性分析，确认各个变量之间是否存在关联。如表1，表 2 所示，本文分析了各年内所有指标之间的相关性。

  根据表 1，表 2，所有的相关性均大于 0.5，可知各因素之间存在相关性，且均为中度相关，通过以上分析可知，本模型给定的五个指标之间并非完全相互独立，两两之间均存在相关性。

  对于扶贫工作的绩效评价，只参考总分是不合理的，显然，将所有评价指标综合起来考量的模型，相较于单一参考总分，能够更加全面、客观、综合的反应扶贫工作进展和效果。这就要求建立一个能够将五项评价指标进行有效结合，并与总分进行统一的模型，来更好的对帮扶单位成绩进行有效评估[2]。
  首先考虑五项指标的结合问题，由于绩效评价更倾向于考虑村庄的进步幅度，总分的进步幅度可以轻松的用两年的差值𝛥𝑄表示，而此五项指标的进步幅度数据较为分散，但却互有相关。对于这些并非各自独立的指标，本文首先采用主成分分析对各年的五项指标进行数据特征提取，得到因子得分最高的一维数据。




  通过以上两个表格可知，类型 4 的帮扶单位，在脱贫攻坚中成果显著，具有较高的绩效，并且，在全国 160 个帮扶单位中，编号为 149、113、158、115、131、155、54、132 和 45 的帮扶单位具有较好的绩效，排名在全国前十。

  不同于上一个模型，该问题要求模型针对不同帮扶单位的各项指标分别进行业绩分析，并对其帮扶业绩做出评估和排序，因此，如何针对某一帮扶单位同一指标下的所有数据进行特征提取，成为了建立模型的主要方向。本文采用第一问的数据构建方法，将2015 年和 2020 年同一指标下的数据两两组合，构成新的二维数据分布，并利用 K-means聚类分析方法，对各个帮扶单位下的所有村庄情况进行聚类分析，得到各个指标下所有帮扶单位的聚类中心点，通过反向对应数据点，得到其大致的 2015 年和 2020 年的总体指标数据，通过作差的形式，得到其对应的帮扶业绩，并依据帮扶业绩情况进行排序，得出前五名帮扶单位编号。



  首先，如图 2 所示，本文对所有 160 个帮扶单位在五项指标上的聚类中心点进行了绘制。这些数据点的横纵坐标做差后能够直接看出其帮扶业绩的进步情况。对各个指标中所有的帮扶单位两年间的差值进行排序，得到了各指标排名前五的单位编号及其进步幅度，如表格 5 所示。
  根据表 3 可知，在居民收入指标上，帮扶单位 38，25，159，155 和 92 表现最佳，在产业发展指标上，帮扶单位 151，138，23，10 和 131 表现最佳，在居民环境指标上，32，139，45，37 和 102 表现最佳，在文化教育指标上，帮扶单位 14，43，65，42 和148 表现最佳，最后，在基础设施指标上，帮扶单位 151，139，23，10 和 131 的表现最佳，进步幅度较大。
  问题四要求对评选先进的重要影响因素进行分析，并对几个村庄的数据进行预测，以确定其是否能够评上荣誉称号。关于对重要影响因素的分析问题上，本文需要对数据的深层次特征进行提取，同时，更要利用已有的三万五千条数据，对给出的十个村庄的未知数据进行准确预测。显然，就本题而言，每一条数据都包含 2015 年和 2020 年的具体信息，属于有监督学习，在这种情况下进行数据预测，神经网络是不二之选。然而，常规的神经网络对于深层数据挖掘的效果十分有限，且深层神经网络的训练代价较高，对于本文来说从时间角度并不适用。


  值得一提的是，通过观察数据发现，本数据集当中的数据多为离散的浮点数，不适用于传统的二值受限玻尔兹曼机模型，因此，本文对该模型进行改进，将原有隐层后验概率分布——伯努利分布改为高斯分布，称为高斯-伯努利受限玻尔兹曼机[7]，适应了数据特征。
  尔后，本文将预测得到的十个村庄的数据，以第二题模型中的绩效评价方法进行判定，将其得分与全国所有村庄进行比较，排位在前 10000 名之内的，自然为脱贫先进村庄，尔后在其中，按照 1：3 的比例筛选等级。

  首先，考虑 2015 年的五项指标和总分数据为训练集数据，2020 年的五项指标和总分数据作为目标数据，利用高斯-伯努利受限玻尔兹曼机模型，对输入层和隐层之间的权值矩阵进行初步训练，再建立神经网络模型，利用反向传播算法对初始化后的多层感知机进行训练。本文采用了输入层隐层输出层=6 × 12 × 6的网络结构，学习率设为 0.2，收敛目标值设为10的-5次幂，预设迭代次数为 5000 次，训练过程中的误差如图 5 所示，梯度变化情况如图 6 所示。


  神经网络训练完成后，将数据被删除的十个村庄在 2015 年的数据输入神经网络中，得到预测数据如下：

  最后，将以上数据代入到公式(12)，得到十个被删除数据的村庄的预测得分和在全国过所有村庄的排名如表 7 所示：

  根据题目给定条件，排名前 10000 名的村庄可以获得荣誉称号，且在这 10000 个村庄内按照 1：3 的比例划分一等和二等，即 2500：7500，因此，编号为 12916、21570、22096、34208、42883 的村庄将获得二级荣誉称号。由于这些村庄里没有进入前 2500 名的村庄，故这十个村庄均不能获得一等荣誉称号。
  接下来，本文将对影响荣誉称号评定的因素进行评价。在这里，本文建立了一个结构为10 × 5 × 1的神经网络，学习率设为 0.2，收敛目标值设为10的-5次幂，预设迭代次数为 5000次，训练过程中的误差如图 7 所示，梯度变化情况如图 8 所示。



  通过观察可知，对于最终的绩效评价，2015 年的所有指标所占权重均为负值，而2020 年的所有指标所占权重均为正值，这意味着 2015 年各项指标的数据与最终绩效评价呈负相关，而 2020 年各项指标的数据与最终绩效评价呈正相关，其中，2015 年的指标里，居民收入对最终绩效的影响最大，而对于 2020 年的指标来讲，产业发展则对最终绩效评价具有较大影响。

给国家扶贫办的一封信

尊敬的国家扶贫办工作人员：
   您好！2020 年是全面建成小康社会目标实现之年，是脱贫攻坚收官之年，自 2015年以来，我国脱贫攻坚工作取得了卓越的成绩，各级扶贫单位在脱贫攻坚岗位上取得了不可磨灭的政绩。
  我们的研究根据 2015 年和 2020 年的统计数据从居民收入、产业发展、居住环境、文化教育、基础设施五个指标入手分析得出指标之间都具有相关性，且在时间序列上都具有正相关性。故我们建议在评价帮扶情况时要综合考虑上述多种指标，并且扶贫是长效工作，在评价时需要考虑较长一段时间内的扶贫效果。同时，仅仅通过 2020 年村庄评分数据无法有效的体现帮扶单位在扶贫工作上所做的努力。我们通过统计学方法对不同帮扶单位各项指标的进步幅度进行加权融合，分析得出了类型 4 的帮扶单位的帮扶效果最为显著，故建议可以通过研讨学习的方式，分享类型 4 的帮扶单位的成功经验，使整体的脱贫攻坚工作取得更大的进步。不同的帮扶单位都有着各自的工作特色，在不同方面的工作成效也有所差异。根据我们的研究结果，我们建议让帮扶单位 38，151，32，14，151 分别从居民收入、产业发展、居民环境、文化教育、基础设施五个方面介绍自己的扶贫工作，只有发扬每个扶贫单位的特色才能帮助不同类型的扶贫单位的工作越做越好。通过授予先进称号的激励方法不仅可以表彰先进，更可以起到激励民众在脱贫攻坚过程中的信心，利用高斯-伯努利 RBM 神经网络的方法，我们研究得出 2015 年各项指标的数据与最终绩效评价呈负相关，而 2020 年各项指标的数据与最终绩效评价呈正相关的结论，且在 2015 年的指标里，居民收入对最终绩效的影响最大，而对于 2020 年的指标来讲，产业发展则对最终绩效评价具有较大影响。从结果可以看出，从时间角度上，影响绩效的因素是变化的，故我们建议在不同时间及客观条件下，先进称号的评价体系应该需要灵活调整，或设置多种特色奖项作为奖励。
   脱贫攻坚工作作为我国全面建成小康社会的重点工作，是实现中华民族伟大复兴之路上的重要阶段。希望我们的研究结论可以为脱贫攻坚工作更好的开展贡献一份力量，
谢谢！

论文缩略图：

程序代码：

%%此程序为改进的 BP 神经网络评分程序%%输入目标矩阵
close all
 
clear echo on clc
% NEWFF——生成一个新的前向神经网络	% TRAIN——对 BP 神经网络进行训练
% SIM——对 BP 神经网络进行仿真
pause
%	敲任意键开始
clc
%	定义训练样本矢量	% P 为输入矢量
P = [-1:0.05:1];
% T 为目标矢量
randn('seed'，78341223); T = sin(2*pi*P)+0.1*randn(size(P));
%	绘制样本数据点
plot(P，T，'+'); echo off
hold on; plot(P，sin(2*pi*P)，':');
%	绘制不含噪声的正弦曲线
echo on clc pause clc
%	创建一个新的前向神经网络net=newff(minmax(P)，[20，1]，{'tansig'，'purelin'}); pause
clc  echo off clc
disp('1.	L-M  优化算法 TRAINLM'); disp('2.	贝叶斯正则化算法 TRAINBR');
choice=input('请选择训练算法(1，2):'); figure(gcf);
if(choice==1)
echo on clc
%	采用 L-M 优化算法 TRAINLM net.trainFcn='trainlm';
pause clc
%	设置训练参数
net.trainParam.epochs = 500; net.trainParam.goal = 1e-6; net=init(net);
%	重新初始化
pause clc
elseif(choice==2) echo on
clc
%	采用贝叶斯正则化算法 TRAINBR net.trainFcn='trainbr';
pause clc
%	设置训练参数
net.trainParam.epochs = 500; randn('seed'，192736547); net = init(net);
%	重新初始化
pause clc
end

clear
clc
close all
[a,b,raw] = xlsread('D:\A 论文集\华教杯\赛题\2020 年“华数杯”全国大学生数学建模竞赛
赛题\2020 年“华数杯”全国大学生数学建模竞赛 C 题\附件 C 题数据.xlsx');
data15 = a(:,4:8);
data20 = a(:,10:14);
colorr = {'r.','g.','b.','m.','k.'};
dt2015 = a(:,4:8);
dt2015_zong = a(:,9);
dt2020 = a(:,10:end-1);
dt2020_zong = a(:,end);
[~,pca_dt2015] = pca(dt2015);
pca_dt2015 = pca_dt2015(:,1);
[~,pca_dt2020] = pca(dt2020);
pca_dt2020 = pca_dt2020(:,1);
index_diff = pca_dt2020-pca_dt2015;
zong_diff = dt2020_zong-dt2015_zong;
score = 0.4649*index_diff + (1-0.4649)*zong_diff;
score(:,2) = a(:,1);
[vhW, vb, hb, fvar, errs] = GaussianRBM([data15,data20],params);
%% 神经网络部分
trainset = a(:,4:9)';
targetset = a(:,10:15)';
testset = a(end-9:end,4:9)';
net=newff(trainset,targetset,12,{'tansig','tansig','purelin'},'trainlm');
net.trainParam.epochs = 5000;
net.trainParam.goal=0.0000001;
[net,tr]=train(net,trainset,targetset);
testres=sim(net,testset);
 
%% RBM 部分
function [vhW, vb, hb, fvar, errs] = GaussianRBM(batchdata, params)
[n d nBatches]=size(batchdata);
assert(params.v_var > 0);
fstd = ones(1,d)*sqrt(params.v_var);
params.v_var=[];
r = params.epislonw_vng;
std_rate = linspace(0, params.std_rate, params.maxepoch);
std_rate(:) = params.std_rate;
std_rate(1:min(30, params.maxepoch/2)) = 0; %learning schedule for 
variances
assert( all(size(params.PreWts.vhW) == [d params.nHidNodes]));
assert( all(size(params.PreWts.hb) == [1 params.nHidNodes]));
assert( all(size(params.PreWts.vb) == [1 d]));
vhW = params.PreWts.vhW;
vb = params.PreWts.vb;
hb = params.PreWts.hb;
vhWInc = zeros( d, params.nHidNodes);
hbInc = zeros( 1, params.nHidNodes);
vbInc = zeros( 1, d);
invfstdInc = zeros(1,d);
Ones = ones(n,1);
q=zeros(1, params.nHidNodes); %keep track of average activations
errs = zeros(1, params.maxepoch);
fprintf('\rTraining Learning v_var Gaussian-Binary RBM %d-%d epochs:%d 
r:%f',...
 d, params.nHidNodes, params.maxepoch, r);
for epoch = 1:params.maxepoch
 
 if rem(epoch, int32(params.maxepoch/20)) == 0 || epoch < 30
 fprintf('\repoch %d',epoch);
 end
 
 errsum=0;
 ptot = 0;
 for batch = 1:nBatches
 Fstd = Ones*fstd;
 %%%%%%%%% START POSITIVE PHASE 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 data = batchdata(:,:,batch); %nxd 
 pos_hidprobs = 1./(1 + exp(-(data./Fstd)*vhW - Ones*hb)); %p(h_j 
=1|data) 
 pos_hidstates = pos_hidprobs > rand( size(pos_hidprobs) );
 
 pos_prods = (data./Fstd)'* pos_hidprobs;
 pos_hid_act = sum(pos_hidprobs);
 pos_vis_act = sum(data)./(fstd.^2); %see notes on this derivation
%%%%%%%%% END OF POSITIVE PHASE %%%%%%%%%
 for iterCD = 1:params.nCD
 
 %%%%%%%%% START NEGATIVE PHASE 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 negdataprobs = pos_hidstates*vhW'.*Fstd+Ones*vb;
 negdata = negdataprobs + randn(n, d).*Fstd;
 neg_hidprobs = 1./(1 + exp(-(negdata./Fstd)*vhW - Ones*hb )); 
%updating hidden nodes again
 pos_hidstates = neg_hidprobs > rand( size(neg_hidprobs) );
 
 end %end CD iterations
 
 neg_prods = (negdata./Fstd)'*neg_hidprobs;
 neg_hid_act = sum(neg_hidprobs);
 neg_vis_act = sum(negdata)./(fstd.^2); %see notes for details
 
 %%%%%%%%% END OF NEGATIVE PHASE 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 errsum = errsum + sum(sum( (data-negdata).^2 ));
 
 if epoch > params.init_final_momen_iter,
 momentum=params.final_momen;
 else
 momentum=params.init_momen;
 end
 
 %%%%%%%%% UPDATE WEIGHTS AND BIASES 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 vhWInc = momentum*vhWInc + r/n*(pos_prods-neg_prods) - 
r*params.wtcost*vhW;
 vbInc = momentum*vbInc + (r/n)*(pos_vis_act-neg_vis_act);
 hbInc = momentum*hbInc + (r/n)*(pos_hid_act-neg_hid_act);
 
 invfstd_grad = sum(2*data.*(Ones*vb-data/2)./Fstd,1) + sum(data' .* 
(vhW*pos_hidprobs') ,2)';
 invfstd_grad = invfstd_grad - ( sum(2*negdata.*(Ones*vbnegdata/2)./Fstd,1) + ...
 sum( negdata'.*(vhW*neg_hidprobs') ,2 )' );
 
 invfstdInc = momentum*invfstdInc + std_rate(epoch)/n*invfstd_grad;
 if params.SPARSE == 1 %nair's paper on 3D object recognition 
 %update q
 if batch==1 && epoch == 1
 q = mean(pos_hidprobs);
 else
 q_prev = q;
 q = 0.9*q_prev+0.1*mean(pos_hidprobs);
 end 
 
 p = params.sparse_p;
 grad = 0.1*params.sparse_lambda/n*sum(pos_hidprobs.*(1-
pos_hidprobs)).*(p-q)./(q.*(1-q));
 gradW 
=0.1*params.sparse_lambda/n*(data'./Fstd'*(pos_hidprobs.*(1-
pos_hidprobs))).*repmat((p-q)./(q.*(1-q)), d,1);
 
 hbInc = hbInc + r*grad;
 vhWInc = vhWInc + r*gradW;
 end
 
 ptot = ptot+mean(pos_hidprobs(:));
 
 vhW = vhW + vhWInc;
 vb = vb + vbInc;
 hb = hb + hbInc; 
 
 invfstd = 1./fstd;
 invfstd = invfstd + invfstdInc;
 fstd = 1./invfstd;
 fstd = max(fstd, 0.005); %have a lower bound! 
 %%%%%%%%%%%%%%%% END OF UPDATES 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 end
 if rem(epoch, int32(params.maxepoch/20)) == 0 || epoch < 30
 fprintf(1, ' p%1.2f ', ptot/nBatches );
 fprintf(1, ' error %6.2f stdr:%.5f fstd(x,y): [%2.3f %2.3f] mm:%.2f 
', errsum, std_rate(epoch), fstd(1), fstd(2), momentum);
 fprintf(1, 'vh_W min %2.4f max %2.4f ', min(min(vhW)), 
max(max(vhW)));
 end
 errs(epoch) = errsum; 
end
fvar = fstd.^2；