如何選擇咱們的損失函數 - 鏈接統計推斷和問題所在領域的橋樑

1. 從機率預測到決策函數之間的鴻溝

在統計學領域，統計學家經常不考慮能得到多少，而是考慮損失了多少，咱們認爲損失的負數就是得到。如何衡量損失是一件值得深刻思考的事情。

0x1：例子：只考慮預測結果不考慮不一樣預測結果背後的實際損失帶來的問題

考慮以下例子：

假設一個氣象學家正在預測颶風襲擊他所在城市的機率是多少。他估計颶風不襲擊該城市的可能性爲 99%-100%，且這一律率的可信度爲 95%。氣象學家很是滿意他的計算精度，並建議你們不用疏散。不幸的是颶風襲擊了城市，城市被淹沒了。

這個典型的例子說明了單純使用精度去度量結果時存在的缺陷，很顯然，這是咱們不少開發朋友包括筆者本身平時在項目中無時不刻遇到的狀況。

咱們很幸苦收集了大量打標樣本，設計了一套特徵工程以及算法模型後，又通過了好久的訓練獲得了一個預測模型，這時候，預測模型能夠獲得幾千到幾萬不等的預測結果，其中置信度機率從0.8 ~ 0.99不等，這個時候問題來了，咱們要拿哪些結果進行輸出，例如輸出給下游業務方或者輸出給客戶？

是 >= 0.99？那麼是否意味着漏報了不少數據？

是 >= 0.9？那麼是否意味着產生不少誤報？

退一步說， >= 0.99 就必定是100%沒誤報的嗎？咱們設定的這個所謂的 D（>= 0.99）的決策函數，100%不會出問題嗎？

形成這個問題的根本緣由是什麼呢？

筆者認爲根本緣由是：咱們對損失函數的定義還不夠完善，如本篇blog的tilte所說，損失函數的本質是」鏈接統計推斷和問題域的橋樑「，因此損失函數就不能單純只考慮統計推斷的機率結果，還要考慮問題域的實際狀況，即不一樣的決策結果帶來的實際現實世界的損失，損失函數必須是承上啓下的。

0x2：大體的正確比精確的錯誤更好

從貝葉斯推斷的世界觀來看，咱們不能太過於強調精確度的度量，儘管它每每是一個有吸引力和客觀的度量，但會忽視了執行這項統計推斷的初衷，那就是：推斷的結果。

此外，咱們但願有一種強調「決策收益的重要性」的方法，而不只僅是估計的精確度。

0x3：以尋找最優結果爲目標的損失評價函數

在機器學習項目中，咱們更在乎的是獲得最佳的預測結果，而不只僅是在全部可能的參數下達到最佳精度。
尋找貝葉斯行動等價於尋找一些參數，這些參數優化的不是參數精度，而是任意某種表現。不過咱們須要先定義表現（損失函數、AUC、ROC、準確率 / 召回率）

Relevant Link:

https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers

 

2. 損失函數

在開始討論這章以前，筆者但願先拋出一句話：損失函數是鏈接統計推斷和問題所在領域的橋樑。咱們整篇文章會圍繞這個觀點展開討論。

損失函數不只僅侷限於咱們在機器學習和深度學習書籍裏看到的最小二乘損失、交叉熵損失、0-1損失等等。實際上，損失函數是能夠根據實際狀況，在具體項目中進行自定義的定義，它的核心做用是將現實世界裏的問題經過數學方式數字化表達，讓損失這個抽象的邏輯概念能夠數值化並經過優化算法進行運算和優化。

0x1：損失函數的定義

咱們先來學習統計學和決策理論中的損失函數。損失函數是一個關於真實參數以及對該參數的估計的函數：

損失函數的重要性在於，他們可以衡量咱們的估計的好壞：損失越大，那麼損失函數獲得的估計就越差。

一個簡單而廣泛的例子是平方偏差損失函數，這是一種典型的，與偏差的平方成正比的損失函數。在線性迴歸、無誤差統計量計算、機器學習的許多領域都有普遍的應

0x2：常見的損失函數 

1. 0-1損失函數

0-1損失很是直觀和簡單，非對即錯，是一種階躍函數，在實際項目不多使用到。

2. 平方損失函數 - 非線性損失評價函數

平方損失偏差損失的缺點在於它過於強調大的異常值。這是由於隨着估計值的偏離，損失是平方增長的，而不是線性增長的。

對3個單位偏離的懲罰要小於對5個單位誤差的懲罰；可是對3個單位偏離的懲罰卻不比對1個單位偏離的懲罰大不少，雖然在一些場景下，這兩種狀況下偏離的差值是相同的。

均方損失的這種不線性，意味着較大的偏差會致使糟糕的結果。在實際項目中，是否可使用須要數據科學家本身根據經驗去權衡。

若是咱們的問題是迴歸分析的場景，咱們是但願擬合模型能儘可能地歸納數據的分佈，離得越遠就意味着迴歸擬合能力越差，由於懲罰也越大，所以可能均方損失的這種非線性特性反而是合理的。

所以要選擇損失函數的時候，須要咱們對問題場景要很是的瞭解，知道選擇不一樣的損失函數對最終得到的結果意味着什麼。

3. 對數損失函數（交叉熵損失函數）

交叉熵損失函數基於信息論中的熵值概念，將損失當作是當前系統不肯定性的度量，預測值和目標值差距越大，意味着不肯定性越大，熵值也就越大，損失也就越大。

4. 絕對值損失函數 - 強調偏離度的線性損失評價函數

一種更穩健的損失函數是偏差的絕對值函數，公式以下：

能夠看到，絕對值損失函數不考慮預測值和目標值是」左偏「仍是」右偏「，它只是線性地衡量」預測值和目標值的距離「。

絕對值損失函數是機器學習和穩定統計中常用的損失函數。

筆者插入：能夠很容易地看到，絕對損失函數沒有對偏離的正負進行區分。but！這個重要嗎？截止目前爲止，這個問題筆者很難做出回答，我本身也在摸索中，從筆者目前爲止爲數很少的項目經驗來看，兩種狀況都會遇到，你須要本身靈活判斷。

假如你須要的預測的一臺機器的惡意網絡發包行爲的可疑程度，咱們對網絡established行爲進行建模，進行迴歸分析，迴歸模型擬合的值實際上是一個閾值，低於或者高於這個迴歸值的行爲模式都被認爲是有必定可疑程度的。可是很顯然，從安全業務場景裏去定義這個問題時，低於這個閾值被認爲是可疑程度相對較低的，高於這個閾值的被認爲是可疑程度相對較高的，這個時候，絕對值損失函數就不能100%客觀地表徵我對這個業務場景的理解了，我可能須要將絕對值損失函數進行改造，改成分段函數，並在負分段和正分段分別賦予不一樣的懲罰因子。

咱們在下個小節會討論的分段損失函數，能夠緩解這個問題。

0x3：面向具體業務場景的自定義損失函數

從筆者的理解上來看，目前數據科學家和數據工程師們選擇損失函數主要基於以下兩方面：

1. 數學計算上的便利性，因爲計算機在數學上的便利性，咱們能夠自由地設計本身的損失函數；
2. 應用上的穩健性，即損失函數能夠客觀的對損失進行度量。損失函數確實是客觀的，它們都是估計值和真值之差的函數，無論該偏差是正仍是負，和估計的收益對象無關；
3. 損失函數是」可優化的「，這涉及到凸優化理論，理論上說，只有凸函數能夠經過計算進行優化，可是值得注意的是，如今對非凸函數的優化也有大量的研究，有興趣的朋友能夠google相關資料；

咱們文章開頭說過，損失函數的定義須要承上啓下，既要可以對機率問題進行數值化，也要考慮實際問題域的狀況。這就要求咱們對常規損失函數進行定製。

再次考慮以上颶風的例子，統計學家也能夠這樣預測：

颶風襲擊的機率在 0%到 1%之間，可是若是發生颶風的損失是巨大的。99%的機率沒有颶風，沒有颶風帶來的收益有限。基於這種損失評價體系，統計學家最後給出的建議可能大相徑庭。

很顯然，這種改動意味着咱們把關注的重心從更精確的機率估計轉到機率估計帶來的現實結果上來，即損失函數要更加關注決策結果和業務收益。

1. 分段區間非對稱平方偏差損失函數

設想一種場景，預測值低於目標值形成的損失，小於預測值高於目標值。咱們能夠用一種非對稱的均方偏差損失函數來描述這種關係：

上式代表，這種損失函數對預測值小於目標值的狀況，懲罰更嚴重。在這種損失函數的驅動下，模型會盡可能避免讓預測值小於目標值。

這種損失函數形式更適合於估計下個月的網絡流量，由於爲了不服務器資源分配不足狀況，在對服務器資源分配進行估計時要多分配一些。

從這個例子能夠看到，損失函數並不必定都是簡單的說去衡量模型的預測值和目標值之間的數值差別，而是更加貼近咱們的業務目標。損失函數衡量的是咱們的模型預測結果對咱們的業務目標的貢獻程度。

2. 強調參數估計結果更接近 0 或者 1 的損失函數

該損失函數更傾向於估計靠近0或者1位置的參數 𝜃。或者理解爲噹噹目標值爲靠近0或者1時，對預測值的一點點偏離都會產生比較大的懲罰。

由於若是真值𝜃靠近 0 或者 1 時，損失函數將變的很是大，除非估計值𝜃̂也 接近 0 或者 1。

基於這種損失函數進行參數優化，模型獲得的參數是傾向於靠近 0 或者 1 的。

這種損失函數更適合政治評論家使用，由於他們更傾向於給出 「是/否」的答案（即目標值爲0或者1）。

3. 傾向於偏差小的估計的損失函數

該損失函數的取值範圍爲 0 到 1，這種損失函數代表其使用者不太關心偏差大的估計結果，即預測值和目標值偏離程度並不會顯著致使損失增長，這和 0-1 損失函數比較相似，可是在真值附件的懲罰程度又沒有 0-1 損失函數那麼強烈。0-1損失的懲罰太階躍了。

4. 帶有對某個結果傾向性的損失函數

天氣預報人員使用的損失函數經常是帶有明顯的傾向性的。氣象預報員有很強的動力盡可能準確地預報下雨的機率，也一樣有動力去錯誤地暗示可能有雨。爲何會出現這種傾向性？

人們更喜歡有所準備，即便可能不會下雨，也不喜歡下雨形成的措手不及。出於這個緣由，預報員傾向於人爲地增長降雨機率和誇大這一估計，由於這能帶來更大「收益」。

這個和咱們在開頭討論的洪水預測有相似的意思，正是中國那句老話：寧肯信其有，不可信其無。

 

3. 貝葉斯推斷理論體系下的損失函數

0x1：樣本集數據驅動下的經驗損失函數

須要明白的是，咱們以前在討論損失函數時，都假設咱們知道了真值（目標值），可是實際狀況下，損失函數理論下的參數的真實值基本上是不可能已知的，不然，若是咱們已經知道真實值，也就沒有必要再去估計了。討論和真值之間的損失只是爲咱們提供一個理論的上界，做爲一個benchmark參考系。

所以無論是在傳統機器學習/深度學習中，仍是貝葉斯推斷中，損失函數都不是針對真是參數值的，而是針對訓練樣本集的，即經驗損失函數。

在這一點上，貝葉斯統計推斷和機器學習是同樣的，本質上都是訓練樣本數據驅動的。

0x2：貝葉斯後驗分佈估計與貝葉斯點估計 - 從機率分佈和最大後驗估計的角度看後驗分佈結果

在貝葉斯推斷中，把未知參數當作是與先驗和後驗分佈的隨機變量。貝葉斯推斷不會給出一個具體的估計值，而是該出一個估計值的後驗分佈。

就後驗分佈而言，從後驗分佈獲得的樣本均可能是真實的參數值（真實的參數值可能就是分佈中的某個點）。

而貝葉斯點估計是什麼呢？咱們能夠理解爲在貝葉斯後驗估計的機率分佈空間中取一個點做爲輸出結果。這麼講有點抽象，舉一個例子說明。

一維高斯函數：

 

上圖中的曲線能夠理解爲對均值和方差【u，σ】的貝葉斯後驗估計，它是一個後驗機率分佈。

而咱們根據最大後驗似然估計的原則，選擇【u，σ】做爲這個分佈的最終輸出，獲得的【u，σ】就是一個貝葉斯估計點。

0x3：貝葉斯估計的指望損失 - 貝葉斯估計損失的理論值

當咱們知道了未知參數的後驗分佈，就能夠基於樣本集，計算某種參數估計相關的損失函數了。

可是問題是後驗分佈的損失函數仍是一個分佈，咱們沒法作出決策，咱們須要一個肯定的結果或者實值，例如在不少項目中，當咱們經過ML獲得一個預測機率值的時候，這時咱們面臨一種不肯定性時，咱們每每會經過決策函數（例如一個經驗閾值）來將不肯定性提煉成一個單獨的肯定性決策動做（0/1）。

咱們須要提煉咱們的後驗分佈爲單一的值（或向量（在多變量的狀況下））。若是後驗分佈的值取得合理，能夠避免頻率學派沒法提供不肯定性的缺陷，同時這種結果的信息更豐富。

回到貝葉斯估計獲得的後驗估計話題，比起貝葉斯後驗估計，咱們更關心的是在某種估計下的指望損失，即一個肯定的實值。

這麼作的緣由有不少，首先指望損失有利於解釋貝葉斯點估計。現實世界中的系統和機器沒法把先驗分佈做爲輸入參數。實際上最終關心的是估計結果，後驗分佈只是計算估計結果的中間步驟， 若是直接給出後驗分佈是沒有太大做用的。

與之相似，咱們須要使後驗分佈變得陡峭，理想狀況是變成一個點。若是後驗分佈的值取得合理，能夠避免頻率學派沒法提供不肯定性的缺陷，這種結果的信息更豐富。

從貝葉斯後驗中選擇點， 就是貝葉斯點估計，也即求指望損失。

假設 𝑃(𝜃|𝑋) 是在觀測數據 X 的條件下 𝜃 的後驗機率，𝜃 的指望損失  爲

從後驗分佈獲得 𝑁 個樣本 𝜃𝑖 , 𝑖 = 1, ⋯ , 𝑁， 損失函數爲𝐿，能夠經過大數定理，利用 𝜃̂ 的估計值近似計算出指望損失

筆者插入：這裏要注意，指望損失計算的對象貝葉斯估計的理論值，在實際工程中，咱們是沒有上帝視角，咱們只有訓練樣本，所以基於訓練樣本獲得的貝葉斯後驗分佈以及貝葉斯估計結果，自己就包含了必定的偏置（bias），咱們將基於樣本獲得的貝葉斯估計結果成爲經驗損失。

0x2：貝葉斯點估計相比於頻率派估計的優點

注意，與 MAP 相比，計算損失函數的指望值利用了後驗分佈中的更多信息，由於在 MAP 中僅僅是尋找後驗分佈的最大值，而忽略了後驗分佈的形狀。

MAP 中忽略掉的信息可能使你暴露在長尾部分的風險當中，像颶風這種可能性很小可是存在的風險卻很巨大，MAP 將致使估計結果無視參數的未知性。

相似地，頻率學派的目標只是最小化化偏差，不考慮與偏差結果相關的損失。頻率派方法幾乎能夠保證永遠不會絕對地準確（咱們都是預測模型幾乎不多能獲得100%的預測結果）。

貝葉斯點估計方法是經過提早計劃解決這個問題：若是你的估計將是錯誤的，那還不如以正確的姿態犯錯 --- 模糊地正確性賽過精確的錯誤。

 

4. 貝葉斯統計實例：優化「價格競猜」遊戲的展品出價

這個小節咱們使用貝葉斯推斷中的核心思想：使用基於面向最終結果收益的損失函數，進行策略的優化。

0x1：業務問題場景

咱們的目標是在一個價格競猜的遊戲中取得最大收益，「價格競猜」的遊戲規則以下：

1. 比賽雙方爭奪競猜站臺上的商品的價格。
2. 每位參賽者看到的商品都不同，價格都是獨一無二的，因此並不存在價格各項干擾的狀況。
3. 觀看後，每位參賽者被要求給出對於本身那套獎品的投標價格。
4. 若是投標價格超過實際價格，投標者被取消獲獎資格。
5. 若是投標價格低於真正的價格，且差距在250$之內，投標者得到兩套獎品。

遊戲的難度在於如何平衡價格的不肯定性，保持你的出價足夠低，以便不超過實際價格，同時又不能過低，須要儘可能接近真實價格。

0x2：對問題進行數學建模

很明顯，這個問題是一個迴歸預測問題，解決這個問題的方法有不少種，咱們可使用機器學習迴歸分析，也能夠構建一個深度學習網絡進行有監督學習預測，同時，咱們也可使用本章介紹的主題：貝葉斯機率，將咱們的先驗知識、訓練樣本、後驗估計結果所有機率化，在機率分佈的框架內對問題進行建模分析。

1. 肯定隨機變量

無論是機器學習/深度學習/機率統計方法，在建模的時候都暗含了一個要求，即確認隨機變量。

因此接下來的問題是，在此次編程建模中，隨機變量是什麼？

1. 基於歷史價格對對待預測價格的先驗估計：是輸入x；
2. 後驗估計分佈：是預測值（y）；
3. 對其的信心分佈參數：是函數模型的權重參數；
4. 在訓練中觀察到兩個獎項更新後的真實價格：是Y目標值；

2. 肯定先驗分佈

先驗分佈的選擇，是一個很是複雜的話題，在PAC可學習理論中，先驗假設是對假設搜索空間的一種約束，先驗越強，這種約束就越強。先驗分佈的好處是使搜索結果更快且更靠近先驗約束，缺點就是若是先驗引入的很差，可能會致使最後的搜索結果欠擬合。

這裏咱們假設咱們記錄了以前的「價格競猜」比賽，咱們將歷史的競猜結果做爲真實價格的先驗分佈。

咱們假設它遵循一個正態分佈：真實價格 ~ Normal（𝜇p, 𝜎p）

經過以前的記錄，獲得的先驗假設𝜇p = 35000，𝜎p = 7500。

筆者思考_1：注意，這裏咱們的先驗不是一個肯定的實值，而是一個機率分佈，在貝葉斯機率編程體系中，讀者朋友必定要深入去理解這種」一切皆有可能「的思想，不管是先驗知識，仍是後驗估計，全部的東西都是一個機率分佈，世界不是肯定的，而是機率的。

筆者思考_2：這個問題是一個憑空預測價格的場景，並非提供商品的某些特徵而後基於這些特徵來預測價格，因此這裏其實暗含了一個前提，全部待預測價格的商品，價格都是接近的，收斂在某個價格區間中，因此不管是先驗分佈仍是信念分佈，訓練過程本質上都是在儘可能地擬合（靠近）那個價格區間，若是沒有這個前提，這個問題是無解的，或者至少不能用簡單的一個信念函數來解，請讀者朋友仔細體察這點。

3. 決策函數的選擇 - 函數模型F() - 信念函數

咱們把先驗分佈理解爲機器學習中的輸入x，決策函數至關於F(x)的形式，決策函數負責將先驗分佈轉化爲後驗機率估計。同時它也是決定咱們應該怎麼玩遊戲的策略判斷依據。

咱們如今有一個關於價格的大概想法（即先驗分佈估計），但這種猜想可能與真實價格顯著不一樣（在舞臺上壓力會成倍增長，你能夠看到爲何有些出價這麼瘋狂）。因此咱們須要定義一個信念函數，這個信念函數決定咱們有多大可能性會相信本身的先驗判斷。

咱們假設咱們關於獎品價格的信念也是正態分佈：

Prize i~ Normal（𝜇𝑖, 𝜎𝑖），i = 1，2。這裏 i = 1，2 表示每一個參賽者本身的那套獎品只有兩個獎品（但這能夠擴展到任何數量）。

該套獎品的價格，能夠由 Prize1 + Prize2 + w 決定，其中 w 是一個偏差項。

訓練的結果是獲得一個最匹配訓練樣本（歷史價格分佈）的模型參數，即信念函數，這個信念函數是一個連續正態函數，根據不一樣的價格給出不一樣的信念，基本上來講，越偏離信念均值中心的價格，信念就越低。

4. 肯定後驗估計函數 - 對真實價格估計的模型

首先明確一點，貝葉斯估計是面向機率分佈的一種算法，算法獲得的後驗估計不是一個肯定的實值，而是一個分佈。

咱們選擇高斯分佈做爲後驗估計的模型，理由是高斯分佈很是適合進行決策。高斯分佈的均值中心自然就包含了最佳後驗機率估計的特性。

同時由於獎品有2個，全部最終的結果是兩個分佈累加的結果，即後驗估計模型函數F(x)是一個複合函數，F(x) = F(x1) + F(x2)。 

0x3：基於PyMC預測獎品價格實例

1. 先驗價格分佈

讓咱們取一些具體的值，假設套裝有兩個獎品：一趟奇妙的加拿大多倫多之旅，一個可愛的新吹雪機。

咱們對這些獎品的真實價格有一些猜想（先驗假設），可是咱們也很是不肯定。按照上面所討論的，咱們能夠用正態分佈來表達這種不肯定性。

即吹雪機 ~ Normal（3000，500），多倫多之旅 ~ Normal（12000，3000）

這個先驗假設包含了一個事實，咱們相信前往多倫多旅行的真實價格爲 12000 美圓，但有68.2% 的機率價格會降低1個標準差，也就是說，有 68.2% 的機率行程價格在 [9000,15000]區間中。

# -*- coding: utf-8 -*-

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt


norm_pdf = stats.norm.pdf

plt.subplot(311)
x = np.linspace(0, 60000, 200)
# 歷史價格先驗分佈
sp1 = plt.fill_between(x, 0, norm_pdf(x, 35000, 7500),
                       color="#348ABD", lw=3, alpha=0.6,
                       label="historical total prices")
p1 = plt.Rectangle((0, 0), 1, 1, fc=sp1.get_facecolor()[0])
plt.legend([p1], [sp1.get_label()])

plt.subplot(312)
x = np.linspace(0, 10000, 200)
# 吹雪機價格先驗分佈
sp2 = plt.fill_between(x, 0, norm_pdf(x, 3000, 500),
                       color="#A60628", lw=3, alpha=0.6,
                       label="snowblower price guess")

p2 = plt.Rectangle((0, 0), 1, 1, fc=sp2.get_facecolor()[0])
plt.legend([p2], [sp2.get_label()])

plt.subplot(313)
x = np.linspace(0, 25000, 200)
# 多倫多之旅價格先驗分佈
sp3 = plt.fill_between(x, 0, norm_pdf(x, 12000, 3000),
                       color="#7A68A6", lw=3, alpha=0.6,
                       label="Trip price guess")
plt.autoscale(tight=True)
p3 = plt.Rectangle((0, 0), 1, 1, fc=sp3.get_facecolor()[0])
plt.legend([p3], [sp3.get_label()])

plt.show()

2. PYMC訓練及預測後驗機率估計

# -*- coding: utf-8 -*-

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm

norm_pdf = stats.norm.pdf

plt.subplot(311)
x = np.linspace(0, 60000, 200)
# 歷史價格先驗分佈
sp1 = plt.fill_between(x, 0, norm_pdf(x, 35000, 7500),
                       color="#348ABD", lw=3, alpha=0.6,
                       label="historical total prices")
p1 = plt.Rectangle((0, 0), 1, 1, fc=sp1.get_facecolor()[0])
plt.legend([p1], [sp1.get_label()])

plt.subplot(312)
x = np.linspace(0, 10000, 200)
# 吹雪機價格先驗分佈
sp2 = plt.fill_between(x, 0, norm_pdf(x, 3000, 500),
                       color="#A60628", lw=3, alpha=0.6,
                       label="snowblower price guess")

p2 = plt.Rectangle((0, 0), 1, 1, fc=sp2.get_facecolor()[0])
plt.legend([p2], [sp2.get_label()])

plt.subplot(313)
x = np.linspace(0, 25000, 200)
# 多倫多之旅價格先驗分佈
sp3 = plt.fill_between(x, 0, norm_pdf(x, 12000, 3000),
                       color="#7A68A6", lw=3, alpha=0.6,
                       label="Trip price guess")
plt.autoscale(tight=True)
p3 = plt.Rectangle((0, 0), 1, 1, fc=sp3.get_facecolor()[0])
plt.legend([p3], [sp3.get_label()])

# plt.show()

# 吹雪機和多倫多之旅的價格猜想先驗分佈
data_mu = [3e3, 12e3]
data_std = [5e2, 3e3]

# 對價格的信念分佈就是歷史價格的先驗分佈
mu_prior = 35e3
std_prior = 75e2

true_price = pm.Normal("true_price", mu_prior, 1.0 / std_prior ** 2)

prize_1 = pm.Normal("first_prize", data_mu[0], 1.0 / data_std[0] ** 2)
prize_2 = pm.Normal("second_prize", data_mu[1], 1.0 / data_std[1] ** 2)
price_estimate = prize_1 + prize_2


@pm.potential
def error(true_price=true_price, price_estimate=price_estimate):
        return pm.normal_like(true_price, price_estimate, 1 / (3e3) ** 2)


mcmc = pm.MCMC([true_price, prize_1, prize_2, price_estimate, error])
mcmc.sample(50000, 10000)

# train
price_trace = mcmc.trace("true_price")[:]


# predict
x = np.linspace(5000, 40000)
plt.plot(x, stats.norm.pdf(x, 35000, 7500), c="k", lw=2,
         label="prior dist. of suite price")

_hist = plt.hist(price_trace, bins=35, normed=True, histtype="stepfilled")
plt.title("Posterior of the true price estimate")
plt.vlines(mu_prior, 0, 1.1 * np.max(_hist[0]), label="prior's mean",
           linestyles="--")
plt.vlines(price_trace.mean(), 0, 1.1 * np.max(_hist[0]),
           label="posterior's mean", linestyles="-.")
plt.legend(loc="upper left")

plt.show()

從上圖能夠看到，基於吹雪機和多倫多旅行的價格先驗，整體價格從歷史先驗的35000降低到了20000（圖中藍色色塊）。

3. 結合實際業務場景的損失函數

從最大後驗的角度來看，咱們此時已經解決問題了，直接給出最終出價便可。可是在貝葉斯統計領域，咱們有更多的信息，咱們應該將損失歸入咱們的最終出價中。
咱們須要使用結合業務場景的損失函數來找到最佳出價，即對於咱們的損失函數來講是最優解。

def showcase_loss(guess, true_price, risk=80000):
    if true_price < guess:
        return risk
    elif abs(true_price - guess) <= 250:
        return -2*np.abs(true_price)
    else:
        return np.abs(true_price - guess - 250)

其中 riks 是一個參數，代表若是你的猜想高於真正的價格的損失程度。
這裏選擇了一個值: 80000，風險較低意味着你更能忍受出價高於真實價格的狀況。
若是咱們出價低於真實價格，而且差別小於250，咱們將得到兩套獎品；不然，當咱們出價比真實價格低，咱們要儘量接近，所以其損失是猜想和真實價格之間距離的遞增函數。

接下來的問題是，這個 risk 參數如何選擇呢？

對於每一個可能的出價，咱們計算與該出價相關聯的指望損失，經過改變 risk 參數來看它如何影響咱們的最終損失。

# -*- coding: utf-8 -*-

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm

norm_pdf = stats.norm.pdf

# plt.show()

# 吹雪機和多倫多之旅的價格猜想先驗分佈
data_mu = [3e3, 12e3]
data_std = [5e2, 3e3]

# 對價格的信念分佈就是歷史價格的先驗分佈
mu_prior = 35e3
std_prior = 75e2

true_price = pm.Normal("true_price", mu_prior, 1.0 / std_prior ** 2)

prize_1 = pm.Normal("first_prize", data_mu[0], 1.0 / data_std[0] ** 2)
prize_2 = pm.Normal("second_prize", data_mu[1], 1.0 / data_std[1] ** 2)
price_estimate = prize_1 + prize_2


@pm.potential
def error(true_price=true_price, price_estimate=price_estimate):
        return pm.normal_like(true_price, price_estimate, 1 / (3e3) ** 2)

mcmc = pm.MCMC([true_price, prize_1, prize_2, price_estimate, error])
mcmc.sample(50000, 10000)

# train
price_trace = mcmc.trace("true_price")[:]


def showdown_loss(guess, true_price, risk=80000):
        loss = np.zeros_like(true_price)
        ix = true_price < guess
        loss[~ix] = np.abs(guess - true_price[~ix])
        close_mask = [abs(true_price - guess) <= 250]
        loss[close_mask] = -2 * true_price[close_mask]
        loss[ix] = risk
        return loss


guesses = np.linspace(5000, 50000, 70)
risks = np.linspace(30000, 150000, 6)
expected_loss = lambda guess, risk: \
    showdown_loss(guess, price_trace, risk).mean()

for _p in risks:
    results = [expected_loss(_g, _p) for _g in guesses]
    plt.plot(guesses, results, label="%d" % _p)

plt.title("Expected loss of different guesses, \nvarious risk-levels of \
overestimating")
plt.legend(loc="upper left", title="Risk parameter")
plt.xlabel("price bid")
plt.ylabel("expected loss")
plt.xlim(5000, 30000)

plt.show()

 

最大限度地減小咱們的損失是咱們的目標，這對應於上面圖中的每條曲線的最小值點。更正式地說，咱們但願經過尋找如下公式的解來減小咱們的指望損失

指望損失的最小值被稱爲貝葉斯行動。咱們能夠基於SciPy的fmin api進行智能搜索，尋找任一單變量或多變量函數的極值（不必定是全局極值）

# -*- coding: utf-8 -*-

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.optimize as sop

norm_pdf = stats.norm.pdf

# plt.show()

# 吹雪機和多倫多之旅的價格猜想先驗分佈
data_mu = [3e3, 12e3]
data_std = [5e2, 3e3]

# 對價格的信念分佈就是歷史價格的先驗分佈
mu_prior = 35e3
std_prior = 75e2

true_price = pm.Normal("true_price", mu_prior, 1.0 / std_prior ** 2)

prize_1 = pm.Normal("first_prize", data_mu[0], 1.0 / data_std[0] ** 2)
prize_2 = pm.Normal("second_prize", data_mu[1], 1.0 / data_std[1] ** 2)
price_estimate = prize_1 + prize_2


@pm.potential
def error(true_price=true_price, price_estimate=price_estimate):
        return pm.normal_like(true_price, price_estimate, 1 / (3e3) ** 2)

mcmc = pm.MCMC([true_price, prize_1, prize_2, price_estimate, error])
mcmc.sample(50000, 10000)

# train
price_trace = mcmc.trace("true_price")[:]


ax = plt.subplot(111)


def showdown_loss(guess, true_price, risk=80000):
        loss = np.zeros_like(true_price)
        ix = true_price < guess
        loss[~ix] = np.abs(guess - true_price[~ix])
        close_mask = [abs(true_price - guess) <= 250]
        loss[close_mask] = -2 * true_price[close_mask]
        loss[ix] = risk
        return loss

guesses = np.linspace(5000, 50000, 70)
risks = np.linspace(30000, 150000, 6)
expected_loss = lambda guess, risk: \
    showdown_loss(guess, price_trace, risk).mean()


for _p in risks:
    _color = next(ax._get_lines.prop_cycler)
    _min_results = sop.fmin(expected_loss, 15000, args=(_p,),disp = False)
    _results = [expected_loss(_g, _p) for _g in guesses]
    plt.plot(guesses, _results , color = _color['color'])
    plt.scatter(_min_results, 0, s = 60, \
                color= _color['color'], label = "%d"%_p)
    plt.vlines(_min_results, 0, 120000, color = _color['color'], linestyles="--")
    print("minimum at risk %d: %.2f" % (_p, _min_results))

plt.title("Expected loss & Bayes actions of different guesses, \n \
various risk-levels of overestimating")
plt.legend(loc="upper left", scatterpoints=1, title="Bayes action at risk:")
plt.xlabel("price guess")
plt.ylabel("expected loss")
plt.xlim(7000, 30000)
plt.ylim(-1000, 80000)

plt.show()

當咱們下降風險閾值（不那麼擔心出價太高），咱們能夠提升咱們的出價，想要更接近真實價格。

這個例子中，由於維度較低，咱們能夠肉眼找到極值點，可是在高維函數中，肉眼是沒法找到的極值的。
這裏列舉幾個損失函數的貝葉斯行動能夠用顯式公式表達的例子：

1. 使用均方損失: 貝葉斯行動是後驗分佈的均值；
2. 當後驗分佈的中位數將絕對指望損失函數最小化時，用樣本的中位數來接近是很是準確的；
3. MAP估計是某個損失函數收縮到0-1損失的解；

Relevant Link:

https://www.zhihu.com/question/21134457
http://www.xuyankun.cn/2017/05/13/bayes/
http://mindhacks.cn/2008/09/21/the-magical-bayesian-method/

 

5. 貝葉斯行動預測的另外一個例子 - 金融預測

0x1：問題場景

咱們須要設計一個模型，對某隻股票將來的價格進行預測，咱們的利潤和損失將直接依賴於基於預測價格的行動。
很顯然，若是預測將來價格是上行的(價格上漲)，而實際也上漲，那不管是精確預測仍是預測結果高於實際價格，至少是盈利的；可是若是股票價格實際下行了，預測倒是上行，那麼不管預測結果是靠近仍是偏離，都是虧損的，區別只是虧多虧少而已，固然，虧少確定比虧多要好。
同時還有一個約束條件，做爲金融機構來講，穩健的收益回報比大比例的盈利和虧損更被看重，因此即便是預測上行對了，也儘可能會避免過大的偏離實際值。

0x2：模型設計

1. 肯定損失評價函數

如何衡量模型的預測結果和將來實際價格之間差距的損失？假如咱們使用平方偏差損失函數，那麼對正負號是不加以區分的，預測值 -0.01 和 0.03 和實際價格 0.02 之間的懲罰相同: 

(0.01 - (-0.01))2 = (0.01 - 0.03)2 = 0.0004

若是咱們基於採用了這種損失策略的模型的預測下了堵住，那麼 0.03 的預測值會使你賺錢，而 -0.01 的預測值會使你賠錢，但損失函數沒法提示這一點。咱們須要一個更好的損失函數，即考慮了預測價格的正負號和真正的價值的損失函數。

# -*- coding: utf-8 -*-

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.optimize as sop


def stock_loss(true_return, yhat, alpha=100.):
    if true_return * yhat < 0:
        # opposite signs, not good
        return alpha * yhat ** 2 - np.sign(true_return) * yhat \
            + abs(true_return)
    else:
        return abs(true_return - yhat)


true_value = .05
pred = np.linspace(-.04, .12, 75)

plt.plot(pred, [stock_loss(true_value, _p) for _p in pred],
         label="Loss associated with\n prediction if true value = 0.05", lw=3)
plt.vlines(0, 0, .25, linestyles="--")

plt.xlabel("prediction")
plt.ylabel("loss")
plt.xlim(-0.04, .12)
plt.ylim(0, 0.25)

true_value = -.02
plt.plot(pred, [stock_loss(true_value, _p) for _p in pred], alpha=0.6,
         label="Loss associated with\n prediction if true value = -0.02", lw=3)
plt.legend()
plt.title("Stock returns loss if true value = 0.05, -0.02")

plt.show()

上圖中咱們能夠看到幾個關鍵的特徵：

1. 損失函數反應了用戶並不想猜想符號，猜錯符號的損失遠大於猜對符號；
2. 無論是否猜對符號，用戶都不想大幅度猜錯，即大幅度偏離真實值。這也體現了金融機構應對下行風險(預測方向是錯的，量級很大)和上線風險(預測方向是正確的，量級很大)的態度是類似的。二者都被視爲危險的行爲而不被孤立。所以，當咱們進一步遠離真實價格，咱們的損失會增長，但在正確方向上的極端損失相對錯誤方向上的較小。

2. 肯定函數模型F()

咱們將會在被認爲可以預測將來回報的交易信號上作一個迴歸。數據是人工虛構的，由於絕大多數時候金融數據都不是線性的。

下圖中，咱們沿着最小方差線畫出了數據分佈狀況：

# -*- coding: utf-8 -*-

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.optimize as sop


# Code to create artificial data
N = 100
X = 0.025 * np.random.randn(N)
Y = 0.5 * X + 0.01 * np.random.randn(N)

ls_coef_ = np.cov(X, Y)[0, 1] / np.var(X)
ls_intercept = Y.mean() - ls_coef_ * X.mean()

plt.scatter(X, Y, c="k")
plt.xlabel("trading signal")
plt.ylabel("returns")
plt.title("Empirical returns vs trading signal")
plt.plot(X, ls_coef_ * X + ls_intercept, label="Least-squares line")
plt.xlim(X.min(), X.max())
plt.ylim(Y.min(), Y.max())
plt.legend(loc="upper left")

plt.show()

線性迴歸的函數形式以下：

對一個特定的交易信號 x，回報的分佈有以下形式：

對於給定的損失函數，咱們但願找到：

下圖中咱們對比了在平方損失函數和咱們自定義損失函數條件下，不一樣的交易信號的貝葉斯行動。

從上圖中能夠看出：

1. 當交易信號接近爲0時，正負回報均可能出現，咱們最好的行爲是預測接近爲0，也就是說處於中立；
2. 只有當咱們都很是自信時，咱們才進場下注，當對不肯定性感到不安時，選擇不做爲；
3. 當信號變得愈來愈極端時，咱們對正負回報的預測愈來愈自信，自定義損失函數的模型將會收斂到最小二乘；

具備以上特色的模型被稱爲係數預測模型。稀疏模型不是想要想方設法地去擬合數據，稀疏模型試圖找打針對咱們定義的股票損失的最好預測，而最小二乘只試圖找到相對於平方損失下的數據的最佳擬合。

Relevant Link: