python進行敏感性分析(SALib庫)

来源:https://www.cnblogs.com/liuliumei/archive/2022/10/06/16758669.html
-Advertisement-
Play Games

什麼是敏感性分析 敏感性分析(sensitivity analysis)是指從定量分析的角度研究有關因素髮生某種變化對某一個或一組關鍵指標影響程度的一種不確定分析技術。每個輸入的靈敏度用某個數值表示即敏感性指數(sensitivity index) 敏感性指數包括以下幾種: 一階指數:度量單個模型輸 ...


什麼是敏感性分析

 敏感性分析(sensitivity analysis)是指從定量分析的角度研究有關因素髮生某種變化對某一個或一組關鍵指標影響程度的一種不確定分析技術。每個輸入的靈敏度用某個數值表示即敏感性指數(sensitivity index) 敏感性指數包括以下幾種:
  • 一階指數:度量單個模型輸入對輸出方差的貢獻
  • 二階指數:度量兩個模型輸入的相互作用對輸出方差的貢獻
  • 總階指數:度量模型輸入對輸出方差的貢獻,包括一階及更高階
 

什麼是SALib

 

SALib是一個基於python進行敏感性分析的開源庫,SALib提供一個解耦的工作流,意味著它不直接與數學或計算模型交互,SALib 負責使用其中一個採樣函數(sample functions)生成模型輸入,並使用其中一個分析函數(analyze functions)計算模型輸出的靈敏度指數。使用 SALib 進行敏感性分析遵循四個步驟:

  • 確定模型輸入(參數)及採樣範圍
  • 運行採樣函數生成模型輸入
  • 使用生成的輸入評估模型,保存模型輸出
  • 基於模型輸出運行分析函數計算敏感性指數SALib提供了幾種靈敏度分析函數,如Sobol,Morris和FAST。有許多因素決定了哪種方法適用於特定應用。但是無論選擇哪種方法,都只需要用到兩種函數:sample,analyze

案例1

對Ishigami function進行Sobol敏感性分析,因為Ishigami函數表現出很強的非線性和非單調性,所以常用來測試不確定性和敏感性分析方法

 

 

 

1.導入庫

 

SALib的採樣和分析存儲在不同的模塊中,例如導入saltelli採樣函數和sobol分析函數,使用Ishigami作為測試函數,numpy用於存儲模型輸入和輸出

from SALib.sample import saltelli
from SALib.analyze import sobol
from SALib.test_functions import Ishigami
import numpy as np

2.定義模型輸入

problem = {
    'num_vars':3,
    'names':['x1','x2','x3'],
    'bounds':[[-3.14159265359, 3.14159265359],
              [-3.14159265359, 3.14159265359],
              [-3.14159265359, 3.14159265359]]
}

 

3.生成樣本

 

使用saltelli生成樣本



param_values = saltelli.sample(problem,1024)

 

param_values是一個numpy矩陣,其大小為(8192, 3),saltelli會生成N*(2D+2)個樣本,其中N=1024(傳入參數),D=3(模型輸入數量)。參數calc_second_order=False表示不包括二階指數,採樣數變為N*(D+2)


 

param_values

 

 

array([[-3.13238877, -0.77619428, -0.32827189],
       [-0.08283496, -0.77619428, -0.32827189],
       [-3.13238877,  0.3589515 , -0.32827189],
       ...,
       [-0.93572828,  0.80073797,  0.99095159],
       [-0.93572828,  0.81914574,  2.70901007],
       [-0.93572828,  0.81914574,  0.99095159]])
param_values.shape

 

(8192, 3)

4.運行模型

 

SALib不直接參与數學或計算模型的評估,如果模型是用python書寫,可以直接迴圈遍歷每個樣本輸入和評估模型



 

Y = np.zeros([param_values.shape[0]])
for i,X in enumerate(param_values):
    Y[i] = evaluate_model(X)

如果模型不是python書寫,可以保存模型的輸入輸出



np.savetext('param_values.txt',param_values)
Y = np.loadtxt('outputs.txt',float)

本例中使用Ishigami函數評估樣本數據


 

Y = Ishigami.evaluate(param_values)
Y
array([ 3.426362  ,  3.3527401 ,  0.85463176, ...,  2.72470174,
       -1.40463805,  2.85339365])

5.分析

 

在得到模型的輸出後可以計算敏感性指數。本例中使用sobol.analyze,會計算一階,二階和總階指數

Si = sobol.analyze(problem,Y,print_to_console=True)
          ST   ST_conf
x1  0.555860  0.080045
x2  0.441898  0.034177
x3  0.244675  0.025569
          S1   S1_conf
x1  0.316832  0.068707
x2  0.443763  0.046636
x3  0.012203  0.064176
                S2   S2_con
(x1, x2)  0.009254  0.093058
(x1, x3)  0.238172  0.111655
(x2, x3) -0.004888  0.066105

 

Si是一個字典,關鍵詞有"S1", “S2”, “ST”, “S1_conf”, “S2_conf”, and “ST_conf”。_conf存儲相應的置信區間,置信水平在95%。可以使用print_to_console=True 列印所有的指數,或者直接取鍵值。



Si
{'S1': array([0.31683154, 0.44376306, 0.01220312]),
 'S1_conf': array([0.06314249, 0.05230396, 0.05764901]),
 'ST': array([0.55586009, 0.44189807, 0.24467539]),
 'ST_conf': array([0.08582851, 0.04184123, 0.02424759]),
 'S2': array([[        nan,  0.00925429,  0.23817211],
        [        nan,         nan, -0.0048877 ],
        [        nan,         nan,         nan]]),
 'S2_conf': array([[       nan, 0.08325501, 0.10813299],
        [       nan,        nan, 0.06117807],
        [       nan,        nan,        nan]])}
print(Si['S1'])
[0.31683154 0.44376306 0.01220312]

可以看出x1和x2表現出了一階靈敏性,但是x3沒有一階效應

如果總階指數基本上比一階指數大,則可能發生了高階交互作用,可以查看二階指數

print('x1-x2:',Si['S2'][0,1])
print('x1-x3:',Si['S2'][0,2])
print('x2-x3:',Si['S2'][1,2])
x1-x2: 0.00925429303490799
x1-x3: 0.2381721095685646
x2-x3: -0.004887704633467273

 


x1和x3之間有較強的交互,有時也會出現計算誤差,如x2-x3指數為負,隨著樣本的增加,這些誤差會縮小。
輸出也可以變成Pandas DataFrame從而進行其它分析

total_si,first_si,second_si = Si.to_df()
second_si

 

 

 

6.繪圖

 

為了方便起見,SALib提供了基本的繪圖功能




Si.plot()

 

 

 

案例2

 當你想要分析的模型依賴於不是敏感性分析的參數(如位置和時間),則可以單獨對每個時間/位置進行分析。以拋物線函數為例

 

 

參數a,b將接受敏感性分析,但是x不會
首先導入需要的庫

import numpy as np
import matplotlib.pyplot as plt

from SALib.sample import saltelli
from SALib.analyze import sobol

定義拋物線函數

def parabola(x,a,b):
    return a + b*x**2

字典描述只包含a,b的問題

problem = {
    'num_vars':2,
    'names':['a','b'],
    'bounds':[[0,1]]*2
}

採樣,評估,分析。此例中選舉100個x的取值,針對需要進行敏感性分析的a,b的每個樣本(共384個)計算y值,因此y的大小為(384, 100)

param_values = saltelli.sample(problem,2**6)
print(param_values.shape)
(384, 2)
x = np.linspace(-1,1,100)
y = np.array([parabola(x,*params) for params in param_values])
print(x.shape)
print(y.shape)
(100,)
(384, 100)
print(y)
[[0.421875   0.40593913 0.39032847 ... 0.39032847 0.40593913 0.421875  ]
 [1.21875    1.20281413 1.18720347 ... 1.18720347 1.20281413 1.21875   ]
 [0.859375   0.82594091 0.79318915 ... 0.79318915 0.82594091 0.859375  ]
 ...
 [0.640625   0.61531508 0.59052169 ... 0.59052169 0.61531508 0.640625  ]
 [1.21875    1.19219021 1.16617246 ... 1.16617246 1.19219021 1.21875   ]
 [1.1875     1.16219008 1.13739669 ... 1.13739669 1.16219008 1.1875    ]]

 

此例中的敏感性指數是一個長度為100的列表,每個元素是一個如上例中的字典

sobol_indices = [sobol.analyze(problem,Y) for Y in y.T]
sobol_indices[0]
{'S1': array([0.49526584, 0.49526584]),
 'S1_conf': array([0.17061222, 0.21853518]),
 'ST': array([0.49745084, 0.49672251]),
 'ST_conf': array([0.15376801, 0.16325137]),
 'S2': array([[       nan, 0.00436999],
        [       nan,        nan]]),
 'S2_conf': array([[       nan, 0.39998034],
        [       nan,        nan]])}
len(sobol_indices)
100

 

接下來單獨分析每個x對應的指數,提取一階指數繪圖

# 提取100個a,b一階指數
S1s = np.array([s['S1'] for s in sobol_indices])

fig = plt.figure(figsize=(10,6),constrained_layout = True)
gs = fig.add_gridspec(2,2)

ax0 = fig.add_subplot(gs[:,0])
ax1 = fig.add_subplot(gs[0,1])
ax2 = fig.add_subplot(gs[1,1])

for i,ax in enumerate([ax1,ax2]):
    ax.plot(x,S1s[:,i],
           label=r'S1$_\mathregular{{{}}}$'.format(problem["names"][i]),
           color = 'black')
    ax.set_xlabel('x')
    ax.set_ylabel('First-order Sobol index')
    ax.set_ylim(0,1.04)
    ax.yaxis.set_label_position("right")
    ax.yaxis.tick_right()
    ax.legend(loc='upper right')
    
ax0.plot(x,np.mean(y,axis=0),label="Mean", color='black')
prediction_interval = 95
ax0.fill_between(x,
                 np.percentile(y, 50 - prediction_interval/2., axis=0),
                 np.percentile(y, 50 + prediction_interval/2., axis=0),
                 alpha=0.5, color='black',
                 label=f"{prediction_interval} % prediction interval")
ax0.set_xlabel("x")
ax0.set_ylabel("y")
ax0.legend(title=r"$y=a+b\cdot x^2$",
           loc='upper center')._legend_box.align = "left"

plt.show()

 

 

左圖為每個x對應不同a,b取值下y的均值以及95%置信區間,右圖為參數a,b的一階指數。由圖可知,在x=0時,y完全由參數a決定,參數b由於x而消失,x的絕對值越大,參數b對變化貢獻越大,參數a相應越小




您的分享是我們最大的動力!

-Advertisement-
Play Games
更多相關文章
  • 今天我們來聊一個老生常談的話題,跨域!又是跨域,煩不煩 ?網上跨域的文章那麼多,跨的我眼睛都疲勞了,不看了不看了
  • 一篇文章帶你掌握主流服務層框架——SpringMVC 在之前的文章中我們已經學習了Spring的基本內容,SpringMVC隸屬於Spring的一部分內容 但由於SpringMVC完全針對於服務層使用,所以我們在介紹時常常把SpringMVC單獨當作一個大章節來學習 溫馨提醒:在學習SpringMV ...
  • 介紹了布隆過濾器的原理,結合分析guava框架如何實現JVM層面的布隆過濾器,參照guava編寫Redis實現的分散式布隆過濾器 ...
  • 本文詳細介紹Kmeans聚類演算法的原理和程式實現。首先介紹利用該演算法的原理及理解,詳細介紹基於MATLAB設計一個自定義的Kmeans函數過程,然後利用該函數對UCI的數據集進行聚類以測試聚類結果。代碼見文末介紹,後續章節將介紹的主要部分有:Kmeans演算法的原理與理解、基本原理、演算法流程、編程實現... ...
  • 什麼是API API (Application Programming Interface) :應用程式編程介面 java中的API 指的就是 JDK 中提供的各種功能的 Java類,這些類將底層的實現封裝了起來,我們不需要關心這些類是如何實現的,只需要學習這些類如何使用即可,我們可以通過幫助文檔來 ...
  • admin後臺管理 django給您提供了一個可視化圖形界面,來方便您來對資料庫里的表進行增刪改查的管理 但是!使用admin後臺管理你自己註冊的模型表時,需要自行進行先註冊該表! 在應用下的admin.py里進行註冊: from django.contrib import admin from a ...
  • Google翻譯大概是目前機器翻譯中翻譯最為準確的了,本文分別使用了兩種可行的方式編寫了簡易的谷歌翻譯軟體。將詳細介紹調用谷歌翻譯API和自行定義谷歌翻譯介面的方式,最後講解如何通過pyqt5實現軟體UI界面並附上全部程式文件。要點如下:調用Google翻譯API方式、自行定義谷歌翻譯介面方式、谷歌... ...
  • 這篇博文在早前本人寫的介紹拼圖游戲的基礎上推出帶有GUI用戶界面的增強版,這裡將通過上、中、下三篇博文詳細介紹利用MATLAB GUI設計的拼圖游戲完整實現過程,每篇都會附上相應代碼及解釋。中篇主要講解拼圖游戲中游戲難度選擇(拼圖階數設置)與拼圖塊數字提示功能的詳細實現過程。中篇的要點如下:拼圖游戲... ...
一周排行
    -Advertisement-
    Play Games
  • 1.部署歷史 猿友們好,作為初來實習的我,已經遭受社會的“毒打”,所以請容許我在下麵環節適當吐槽,3Q! 傳統部署 ​ 回顧以往在伺服器部署webapi項目(非獨立發佈),dotnet環境、守護進程兩個逃都逃不掉,正常情況下還得來個nginx代理。不僅僅這仨,可能牽扯到yum或npm。node等都要 ...
  • 隨著技術的進步,跨平臺開發已經成為了標配,在此大背景下,ASP.NET Core也應運而生。本文主要基於ASP.NET Core+Element+Sql Server開發一個校園圖書管理系統為例,簡述基於MVC三層架構開發的常見知識點,前一篇文章,已經簡單介紹瞭如何搭建開發框架,和登錄功能實現,本篇... ...
  • 這道題只要會自定義cmp恰當地進行排序,其他部分沒有什麼大問題。 上代碼: 1 #include<bits/stdc++.h> 2 using namespace std; 3 int n,s,h1,h2,cnt; 4 struct apple{ 5 int height,ns;//height為蘋 ...
  • 這篇文章主要描述RPC的路由策略,包括為什麼需要請求隔離,為什麼不在註冊中心中實現請求隔離以及不同粒度的路由策略。 ...
  • 簡介: 中介者模式,屬於行為型的設計模式。用一個中介對象來封裝一系列的對象交互。中介者是各對象不需要顯式地相互引用,從而使其耦合鬆散,而且可以獨立地改變他們之間的交互。 適用場景: 如果平行對象間的依賴複雜,可以使用中介者解耦。 優點: 符合迪米特法則,減少成員間的依賴。 缺點: 不適用於系統出現對 ...
  • 【前置內容】Spring 學習筆記全系列傳送門: Spring學習筆記 - 第一章 - IoC(控制反轉)、IoC容器、Bean的實例化與生命周期、DI(依賴註入) Spring學習筆記 - 第二章 - 註解開發、配置管理第三方Bean、註解管理第三方Bean、Spring 整合 MyBatis 和 ...
  • 簡介: 享元模式,屬於結構型的設計模式。運用共用技術有效地支持大量細粒度的對象。 適用場景: 具有相同抽象但是細節不同的場景中。 優點: 把公共的部分分離為抽象,細節依賴於抽象,符合依賴倒轉原則。 缺點: 增加複雜性。 代碼: //用戶類 class User { private $name; fu ...
  • 這次設計一個通用的多位元組SPI介面模塊,特點如下: 可以設置為1-128位元組的SPI通信模塊 可以修改CPOL、CPHA來進行不同的通信模式 可以設置輸出的時鐘 狀態轉移圖和思路與多位元組串口發送模塊一樣,這裡就不給出了,具體可看該隨筆。 一、模塊代碼 1、需要的模塊 通用8位SPI介面模塊 `tim ...
  • AOP-03 7.AOP-切入表達式 7.1切入表達式的具體使用 1.切入表達式的作用: 通過表達式的方式定義一個或多個具體的連接點。 2.語法細節: (1)切入表達式的語法格式: execution([許可權修飾符] [返回值類型] [簡單類名/全類名] [方法名]([參數列表]) 若目標類、介面與 ...
  • 測試一、虛繼承與繼承的區別 1.1 單個繼承,不帶虛函數 1>class B size(8): 1> + 1> 0 | + (base class A) 1> 0 | | _ia //4B 1> | + 1> 4 | _ib //4B 有兩個int類型數據成員,占8B,基類邏輯存在前面 1.2、單個 ...