本文目录一览：
- 1、如何使用python做统计分析
- 2、python使用hist画频率直方图时，怎样修改填
- 3、怎样用python的matplotlib模块画累积分布图
- 4、图像处理的Python问题，怎么解决
- 5、利用数据和专家模型实现对企业的信用评级
如何使用python做统计分析
Shape Parameters
形态参数
While a general continuous random variable can be shifted and scaled
with the loc and scale parameters, some distributions require additional
shape parameters. For instance, the gamma distribution, with density
γ(x,a)=λ(λx)a−1Γ(a)e−λx,
requires the shape parameter a. Observe that setting λ can be obtained by setting the scale keyword to 1/λ.
虽然一个一般的连续随机变量可以被位移和伸缩通过loc和scale参数，但一些分布还需要额外的形态参数。作为例子，看到这个伽马分布，这是它的密度函数
γ(x,a)=λ(λx)a−1Γ(a)e−λx,
要求一个形态参数a。注意到λ的设置可以通过设置scale关键字为1/λ进行。
Let’s check the number and name of the shape parameters of the gamma
distribution. (We know from the above that this should be 1.)
让我们检查伽马分布的形态参数的名字的数量。（我们知道从上面知道其应该为1）
from scipy.stats import gamma
gamma.numargs
1
gamma.shapes
'a'
Now we set the value of the shape variable to 1 to obtain the
exponential distribution, so that we compare easily whether we get the
results we expect.
现在我们设置形态变量的值为1以变成指数分布。所以我们可以容易的比较是否得到了我们所期望的结果。
gamma(1, scale=2.).stats(moments="mv")
(array(2.0), array(4.0))
Notice that we can also specify shape parameters as keywords:
注意我们也可以以关键字的方式指定形态参数：
gamma(a=1, scale=2.).stats(moments="mv")
(array(2.0), array(4.0))
Freezing a Distribution
冻结分布
Passing the loc and scale keywords time and again can become quite
bothersome. The concept of freezing a RV is used to solve such problems.
不断地传递loc与scale关键字最终会让人厌烦。而冻结RV的概念被用来解决这个问题。
rv = gamma(1, scale=2.)
By using rv we no longer have to include the scale or the shape
parameters anymore. Thus, distributions can be used in one of two ways,
either by passing all distribution parameters to each method call (such
as we did earlier) or by freezing the parameters for the instance of the
distribution. Let us check this:
通过使用rv我们不用再更多的包含scale与形态参数在任何情况下。显然，分布可以被多种方式使用，我们可以通过传递所有分布参数给对方法的每次调用（像我们之前做的那样）或者可以对一个分布对象冻结参数。让我们看看是怎么回事：
rv.mean(), rv.std()
(2.0, 2.0)
This is indeed what we should get.
这正是我们应该得到的。
Broadcasting
广播
The basic methods pdf and so on satisfy the usual numpy broadcasting
rules. For example, we can calculate the critical values for the upper
tail of the t distribution for different probabilites and degrees of
freedom.
像pdf这样的简单方法满足numpy的广播规则。作为例子，我们可以计算t分布的右尾分布的临界值对于不同的概率值以及自由度。
stats.t.isf([0.1, 0.05, 0.01], [[10], [11]])
array([[ 1.37218364, 1.81246112, 2.76376946],
[ 1.36343032, 1.79588482, 2.71807918]])
Here, the first row are the critical values for 10 degrees of freedom
and the second row for 11 degrees of freedom (d.o.f.). Thus, the
broadcasting rules give the same result of calling isf twice:
这里，第一行是以10自由度的临界值，而第二行是以11为自由度的临界值。所以，广播规则与下面调用了两次isf产生的结果相同。
stats.t.isf([0.1, 0.05, 0.01], 10)
array([ 1.37218364, 1.81246112, 2.76376946])
stats.t.isf([0.1, 0.05, 0.01], 11)
array([ 1.36343032, 1.79588482, 2.71807918])
If the array with probabilities, i.e, [0.1, 0.05, 0.01] and the array of
degrees of freedom i.e., [10, 11, 12], have the same array shape, then
element wise matching is used. As an example, we can obtain the 10% tail
for 10 d.o.f., the 5% tail for 11 d.o.f. and the 1% tail for 12 d.o.f.
by calling
但是如果概率数组，如[0.1,0.05,0.01]与自由度数组,如[10,11,12]具有相同的数组形态，则元素对应捕捉被作用，我们可以分别得到10%，5%，1%尾的临界值对于10，11,12的自由度。
stats.t.isf([0.1, 0.05, 0.01], [10, 11, 12])
array([ 1.37218364, 1.79588482, 2.68099799])
Specific Points for Discrete Distributions
离散分布的特殊之处
Discrete distribution have mostly the same basic methods as the
continuous distributions. However pdf is replaced the probability mass
function pmf, no estimation methods, such as fit, are available, and
scale is not a valid keyword parameter. The location parameter, keyword
loc can still be used to shift the distribution.
离散分布的简单方法大多数与连续分布很类似。当然像pdf被更换为密度函数pmf，没有估计方法，像fit是可用的。而scale不是一个合法的关键字参数。Location参数，关键字loc则仍然可以使用用于位移。
The computation of the cdf requires some extra attention. In the case of
continuous distribution the cumulative distribution function is in most
standard cases strictly monotonic increasing in the bounds (a,b) and
has therefore a unique inverse. The cdf of a discrete distribution,
however, is a step function, hence the inverse cdf, i.e., the percent
point function, requires a different definition:
ppf(q) = min{x : cdf(x) = q, x integer}
Cdf的计算要求一些额外的关注。在连续分布的情况下，累积分布函数在大多数标准情况下是严格递增的，所以有唯一的逆。而cdf在离散分布，无论如何，是阶跃函数，所以cdf的逆，分位点函数，要求一个不同的定义：
ppf(q) = min{x : cdf(x) = q, x integer}
For further info, see the docs here.
为了更多信息可以看这里。
We can look at the hypergeometric distribution as an example
from scipy.stats import hypergeom
[M, n, N] = [20, 7, 12]
我们可以看这个超几何分布的例子
from scipy.stats import hypergeom
[M, n, N] = [20, 7, 12]
If we use the cdf at some integer points and then evaluate the ppf at
those cdf values, we get the initial integers back, for example
如果我们使用在一些整数点使用cdf，它们的cdf值再作用ppf会回到开始的值。
x = np.arange(4)*2
x
array([0, 2, 4, 6])
prb = hypergeom.cdf(x, M, n, N)
prb
array([ 0.0001031991744066, 0.0521155830753351, 0.6083591331269301,
0.9897832817337386])
hypergeom.ppf(prb, M, n, N)
array([ 0., 2., 4., 6.])
If we use values that are not at the kinks of the cdf step function, we get the next higher integer back:
如果我们使用的值不是cdf的函数值，则我们得到一个更高的值。
hypergeom.ppf(prb + 1e-8, M, n, N)
array([ 1., 3., 5., 7.])
hypergeom.ppf(prb - 1e-8, M, n, N)
array([ 0., 2., 4., 6.])
python使用hist画频率直方图时，怎样修改填
示例代码：
#概率分布直方图
#高斯分布
#均值为0
mean = 0
#标准差为1，反应数据集中还是分散的值
sigma = 1
x=mean+sigma*np.random.randn(10000)
fig,(ax0,ax1) = plt.subplots(nrows=2,figsize=(9,6))
#第二个参数是柱子宽一些还是窄一些，越大越窄越密
ax0.hist(x,40,normed=1,histtype='bar',facecolor='yellowgreen',alpha=0.75)
##pdf概率分布图，一万个数落在某个区间内的数有多少个
ax0.set_title('pdf')
ax1.hist(x,20,normed=1,histtype='bar',facecolor='pink',alpha=0.75,cumulative=True,rwidth=0.8)
#cdf累计概率函数，cumulative累计。比如需要统计小于5的数的概率
ax1.set_title("cdf")
fig.subplots_adjust(hspace=0.4)
plt.show()
运行结果为：
怎样用python的matplotlib模块画累积分布图
下面的程序绘制随机变量X的累积分布函数和数组p的累加结果
pl.plot(t, X.cdf(t))
pl.plot(t2, np.add.accumulate(p)*(t2[1]-t2[0]))
图像处理的Python问题，怎么解决
imtools.py里面也要有numpy 的引用才对
def histeq(im,nbr_bins=256):
"""对一幅灰度图像进行直方图均衡化"""
#计算图像的直方图
imhist,bins = histogram(im.flatten(),nbr_bins,normed=True)
cdf = imhist.cumsum() #累计分布函数
cdf = 255 * cdf / cdf[-1] #归一化
#使用累计分布函数的线性插值，计算新的像素
im2 = interp(im.flatten(),bins[:-1],cdf)
return im2.reshape(im.shape),cdf
以上代码我定义在imtools.py文件里并且放在了python2.7里
然后我在num.py里引用他
Python code?
1
2
3
4
5
6
7
8
9
10
from PIL import Image
from pylab import *
from numpy import *
import imtools
im= array(Image.open('E:\\daima\\pydaima\\shijue\\tupian1\\gang2.jpg').convert('L'))
im2,cdf =imtools.histeq(im)
出现以下错误：
Traceback (most recent call last):
File "pyshell#56", line 1, in module
a=imtools.histeq(im)
File "E:\daima\pydaima\shijue\imtools.py", line 32, in histeq
NameError: global name 'histogram' is not defined
利用数据和专家模型实现对企业的信用评级
当一家企业申请贷款时，我们会对他进行一个评级，那么企业的评级是怎样产生的呢？首先我们来看一下什么是企业的评级，如下表所示：
从表中可以看出，企业的评级通常形式是AAA、AA、AA-等形式。那么是哪些因素决定了企业的评级呢？本文将简单阐述一个企业评级产生的过程。
首先，我们需要理解评级的含义，例如上文的AAA、AA+、AA-代表了什么?直观的理解肯定是AAA最好，但更为专业的解释是什么呢？我们引用标准普尔来作更为专业的解释：
从上图我们可以看到，评级更多的是对偿还债务能力的一种度量，通俗的来说，就是这个企业是否会产生违约。转化为数学描述，即该企业违约的概率。这样评级其实就转化为了一个二分类问题，我们就可以使用Logistic回归、GBDT、随机森林等机器学习方法来处理。
因此，我们Y值很容易可以判断为“是否违约”，我们X特征变量则可以表示为该“企业的相关数据”。
第一步我们需要确定我们的Y值，即是否违约。违约的定义可以参考《巴塞尔新资本协议》规定，同时需要考虑违约的时效性问题。比如用户在2018年违约了，但是在2019年偿还了，那么他只能算在2018年的违约样本中，如果用户在一年之内有多次违约，我们取最早的违约日期为违约日。因此，评级一般也是半年或者一年发起一次，遇到特殊情况，客户经理也可以手动发起评级。
一般企业的相关数据最重要的当然是财务报表，它主要包括资产负债表、损益表、现金流量表。财务报表的一般数据处理包括指标计算、量纲统一、勾稽关系检查、会计恒等式检测、新旧准则转换计算。财务报表口径又分为合并、本部、汇总，时间周期可以有月度、年度，在数据匹配时也应制定相应的优先级，在保证数据准确性的同时，也尽可能的增加匹配率。
财务报表匹配上有一个意思的问题是时间匹配。例如用户在2019年12月31日违约，那财报我们应该选择哪个时间点的呢？太近太远都会失去建模的意义，通常业内的做法是好用户匹配12个月以前24个月以内能匹配到的最近财报，违约用户则匹配违约时间前6个月以前18个月以内的能匹配到的最近财报。如下图所示：
为了更准确地对企业进行评级，除了财报之外，我们通常会通过购买等渠道，获取企业的征信、工商、税法数据。征信数据除了企业的征信，我们还会获取企业法人、实际控制人的征信数据。
工商、税法该类数据的主要问题是匹配度一般比较差，非常稀疏，通常不太适合应用在建模中，但我们可以在反欺诈、黑名单的第一道关卡中进行使用。
数据准备完成后，就是数据建模过程。通常企业评级的数据建模主要有数据清洗、单变量分析、相关性分析、变量转换、穷举模型、模型确定这几个步骤，本文不再详细描述建模过程，感兴趣的同学可以自行查找资料，建模是比较常规的过程。其中比较有意思的是变量转换，无论是个人还是企业，我们一般都喜欢最后生成一个分数来描述这个特征或者这个主体，简单的方法可以根据四分位法来确定0-100分的划分区间，需要同时考虑该指标是越大越好还是越小越好。
需要注意的是企业评级一般在数据量较小的时候会加入人工专家判断。例如在变量选择和模型确定上，除了考虑模型选择指标还需要考虑模型的可解释性，模型可解释性在金融领域是非常重要的，也涉及一定的监管要求。
由于企业授信一般额度特别大，但违约样本数据不足，业内很少有纯依靠机器学习的方式来进行评级，也要参考很多专家的意见，当我们在进行其他类似项目评级的时候，这种评级策略也给了我们很大的参考。试想，如果现在给你7个专家，你需要他们给你选出最重要的特征属性，你应该如何来汇总和综合他们的意见呢？这里我们使用层次分析法（AHP），方法的相关概念大家可以自行查找资料。我们会提供一种调查问卷的形式，将候选变量给专家们，让他们进行排序，生成一个表格如下：
其中1-9代表重要性越来越低，通过专家们的选择，我们有了一张各专家判断的指标优先级，也可以看出指标三是非常受到专家们的认可的。那么我们如何综合这些专家的意见呢？首先我们要得到每一个专家每一个指标的权重，然后再求得所有专家的每一个指标的平均权重。
计算步骤如下：
假设a 1 、a 2 、a 3 、……、a 9 代表重要性排序后的9个指标，以最重要的a 1 为锚点，各个指标和a 1 的重要性比较关系可以用公式表示为：
通过上面的模型计算后，最终我们产出了一个企业的分数，到这个时候就离最终评级AAA、AA+等等只有一步之遥了。
比较粗放的方式，我们可以直接拉出所有企业分数排序，然后根据你最终要的评级等级数量（AAA、BB、CC等级数量）进行均分，依次从高到低给予一定评级。但这种均匀分割存在一定弊端，例如AAA级和AA+级的间隔就和BBB和BB+级的间隔相同，这显然是不符合常识的。理论上来说，特别好的客户是相对较少的，特别坏的客户也是相对较少的，大量客户应集中在中部。我们自然的想到我们需要的评级区间的分布，Beta分布正好合适。
确定使用Beta分布后，我们需要调整α，β参数，调整参数就需要很多经验、不断尝试以及业务要求和监管政策要求了。比如每个评级等级分布增幅幅度、每个评级等级分布的违约概率，最终产生的整体违约概率等等。通过使用python的scipy.stats包pdf函数，我们可以很轻松地画出这个曲线。然后使用scipy.stats包的cdf函数就可以得到均匀分布点在Beta分布上累计分布函数值。
以α = 5.1，β=5为例，Beta分布曲线为：
最终通过这个分布区间来分割已排序好的所有企业的分数，映射到相应评级，这样我们就产生了一个企业的评级。
当有新企业来申请授信时，提交相应的财报数据，银行对财报进行勾稽检查和其他质量检测，获取相应的征信、工商、税法，应用我们配置好的模型，我们就可以快速的得到一个对应评级。
本文概述了面对一般企业的评级方法，详细阐述了比较有趣的专家模型产生的过程。当然，企业的评级远远比本文描述的复杂，例如遇到金融机构、事业单位、学校、医院等数据量较少且违约样本几乎为零的情况，我们如何对他们进行评级。同时，企业的数据质量相较个人来说存在很多问题，在特征处理时将耗费大量时间。并且企业的授信额度相对较大，模型更多还是参考为主，背后还得依赖于优秀的授信审批部同事。
得到企业的评级只是第一步，后续我们可以根据相关准入政策确定准入评级，并根据评级给予企业不同的额度。同时，评级是有时效性的，需要定期重新评级，并对企业进行实时监控，当突发情况发生时，及时调整评级。
最后，当我们得到了企业评级，预估了企业违约概率，但不同的债项发生违约时，有的可能全部损失，有的可能部分损失，这个应该如何量化呢？将在后续的文章中和大家分享。
