python畫cdf（python畫圖案）

本文目錄一覽：

• 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分布曲線為：

最終通過這個分布區間來分割已排序好的所有企業的分數，映射到相應評級，這樣我們就產生了一個企業的評級。

當有新企業來申請授信時，提交相應的財報數據，銀行對財報進行勾稽檢查和其他質量檢測，獲取相應的徵信、工商、稅法，應用我們配置好的模型，我們就可以快速的得到一個對應評級。

本文概述了面對一般企業的評級方法，詳細闡述了比較有趣的專家模型產生的過程。當然，企業的評級遠遠比本文描述的複雜，例如遇到金融機構、事業單位、學校、醫院等數據量較少且違約樣本幾乎為零的情況，我們如何對他們進行評級。同時，企業的數據質量相較個人來說存在很多問題，在特徵處理時將耗費大量時間。並且企業的授信額度相對較大，模型更多還是參考為主，背後還得依賴於優秀的授信審批部同事。

得到企業的評級只是第一步，後續我們可以根據相關准入政策確定準入評級，並根據評級給予企業不同的額度。同時，評級是有時效性的，需要定期重新評級，並對企業進行實時監控，當突發情況發生時，及時調整評級。

最後，當我們得到了企業評級，預估了企業違約概率，但不同的債項發生違約時，有的可能全部損失，有的可能部分損失，這個應該如何量化呢？將在後續的文章中和大家分享。