pythonでPearsonの相関係数の信頼区間を算出する

pythonで相関係数の信頼区間をだしたい時ってありますよね。
統計系はRかpythonかくらいの流れなのに、（少なくとも僕が調べた限りは・・・）pythonに相関係数の信頼区間を求めるよく使われているfunctionはscipyなどでは定義されていなさそうです。

ですので、自前実装してみました。

Code

def corr_CI(a, b, alpha = 0.95):
    r = stats.pearsonr(a, b)[0]
    n = len(a)
    if n <= 3:
        AssertionError("Not enough amount data")
    z= 0.5*np.log((1+r)/(1-r))
    za = stats.norm.ppf(0.5 + 0.5 * alpha) 
    zl = z - za * math.sqrt(1/(n-3))
    zu = z + za * math.sqrt(1/(n-3))
    rhol = (math.exp(2 * zl) - 1 )/ (math.exp(2 * zl) +1 )
    rhou = (math.exp(2 * zu) - 1 )/ (math.exp(2 * zu) +1 )
    return rhol, rhou

利用例

例えば、iris datasetで利用してみるなら

import pandas as pd
from sklearn import datasets
iris = datasets.load_iris()
df = pd.DataFrame(iris.data, columns=iris.feature_names)
sl = df["sepal length (cm)"]
sw = df["sepal width (cm)"]
corr_CI(sl, sw)

とすると信頼区間の下限、上限が出力される

参考資料

https://bellcurve.jp/statistics/course/9591.html
というかほぼコピペ・・・