質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

88.58%

重回帰分析のマハラノビス距離でエラー

解決済

回答 1

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 966

_hh

score 56

Pythonのjupyter notebookで、UCI machine learning repositoryのあるデータ('flare.data2.txt')を重回帰分析しております。自作関数によりStepAICで変数選択を行ったのち、信頼区間・予測区間を求めようとすると、エラーが出ます。同じ手順で、"abalone.data.txt"というデータを分析した際には全く問題が無いのですが、なぜかわからずにおります。
エラーが出ているのは、マハラノビスの距離を求めるところで、分散共分散行列の逆行列がNANになってしまっている(下記の"mahara is..."以下の出力)ためなのですが、なぜそのようになるのかが分かりません。
引数のベクトルの値は、ちゃんと入っているようですし、"abalone.data.txt"分析時のデバッグでも、ほぼ同様にベクトルの値が出力され、行列式の値もほぼ近いのですが、不思議なことに、その場合はちゃんと分散共分散行列の逆行列が求まっています。。。変数vecの値の入り方が違うのでしょうか?

【環境】Window10 64bit, chrome

def step_aic(model, exog, endog, **kwargs):
    """
    This select the best exogenous variables with AIC
    Both exog and endog values can be either str or list.
    (Endog list is for the Binomial family.)

    Note: This adopt only "forward" selection

    Args:
        model: model from statsmodels.formula.api
        exog (str or list): exogenous variables
        endog (str or list): endogenous variables
        kwargs: extra keyword argments for model (e.g., data, family)

    Returns:
        model: a model that seems to have the smallest AIC
    """

    # exog, endogは強制的にリスト形式に変換しておく
    exog = np.r_[[exog]].flatten()
    endog = np.r_[[endog]].flatten()
    remaining = set(exog)
    selected = []  # 採用が確定された要因

    # 定数項のみのAICを計算
    formula_head = ' + '.join(endog) + ' ~ '
    formula = formula_head + '1'
    aic = model(formula=formula, **kwargs).fit().aic
    print('AIC: {}, formula: {}'.format(round(aic, 3), formula))

    current_score, best_new_score = np.ones(2) * aic

    # 全要因を採択するか,どの要因を追加してもAICが上がらなければ終了
    while remaining and current_score == best_new_score:
        scores_with_candidates = []
        for candidate in remaining:

            # 残っている要因を1つずつ追加したときのAICを計算
            formula_tail = ' + '.join(selected + [candidate])
            formula = formula_head + formula_tail
            aic = model(formula=formula, **kwargs).fit().aic
            print('AIC: {}, formula: {}'.format(round(aic, 3), formula))

            scores_with_candidates.append((aic, candidate))

        # 最もAICが小さかった要因をbest_candidateとする
        scores_with_candidates.sort()
        scores_with_candidates.reverse()
        best_new_score, best_candidate = scores_with_candidates.pop()

        # 候補要因追加でAICが下がったならば,それを確定要因として追加する
        if best_new_score < current_score:
            remaining.remove(best_candidate)
            selected.append(best_candidate)
            current_score = best_new_score

    formula = formula_head + ' + '.join(selected)
    print('The best formula: {}'.format(formula))
    return model(formula, **kwargs).fit()


import pandas as pd
import numpy as np
import scipy as sp
import scipy.stats as stats
import statsmodels.api as sm
import statsmodels.formula.api as smf

d = pd.read_csv('flare.data2.txt',header=None,skiprows=1,delim_whitespace=True,names=('x1','x2','x3','x4','x5','x6','x7','x8','x9','x10','y1','y2','y3'))
d.head()
df2 = pd.get_dummies(d)
model = step_aic(smf.ols,['x4','x5','x6','x7','x8','x9','x10','x1_B','x1_C','x1_D','x1_E','x1_F','x1_H','x2_A','x2_H','x2_K','x2_R','x2_S','x2_X','x3_C','x3_I','x3_O','x3_X'],['y1'],data=df2)



import pandas as pd
import numpy as np
import scipy as sp
import scipy.stats as stats

y=df2['y1']
X=df2[['x1_E', 'x3_I', 'x2_K', 'x4', 'x1_F', 'x1_D', 'x9', 'x6', 'x1_C']]

import statsmodels.api as sm
mod = sm.OLS(y, sm.add_constant(X))
res = mod.fit()
print(res.summary())


Sigma = sp.asmatrix(df2[['x1_E', 'x3_I', 'x2_K', 'x4', 'x1_F', 'x1_D', 'x9', 'x6', 'x1_C']].cov().as_matrix())
Sigma
pp = sp.linalg.det(Sigma)
print("det is",pp)


def Mahala2(vec_x, vec_mean, mat):
    length = mat.shape[0]
    vec = sp.asmatrix((vec_x - vec_mean).values.reshape(length, 1))
    print("vec_x is",vec_x)
    print("vec_mean is",vec_mean)
    print("vec is", vec)
    inv = sp.linalg.inv(mat)
    pp = sp.linalg.det(mat)
    print("det is",pp)
    mahala2 = vec.T.dot(inv.dot(vec))
    print("mahara is",mahala2)
    return mahala2[0, 0]


a1E = df2[['x1_E']].mean() * 1.01
a3I = df2[['x3_I']].mean() * 1.01

a2K = df2[['x2_K']].mean() * 1.01
a4 = df2[['x4']].mean() * 1.01

a1F = df2[['x1_F']].mean() * 1.01
a1D = df2[['x1_D']].mean() * 1.01

a9 = df2[['x9']].mean() * 1.01
a6 = df2[['x6']].mean() * 1.01

a1C = df2[['x1_C']].mean() * 1.01

import matplotlib.pyplot as plt
%matplotlib inline

import scipy as sp

hl_l = sp.linspace(0, 4, 100)
hat_y = []
# x6以外変数固定
(x1_E, x3_I, x2_K, x4, x1_F, x1_D, x9, x1_C) = (a1E,a3I,a2K,a4,a1F,a1D,a9,a1C)

print((x1_E, x3_I, x2_K, x4, x1_F, x1_D, x9, x1_C))
for hl in hl_l:
    X = sp.array([1,x1_E, x3_I, x2_K, x4, x1_F, x1_D, x9, hl, x1_C])
    hat_y.append(X.dot(res.params))
plt.plot(hl_l, hat_y)
n = 1066
phi_e = n - res.df_model - 1
t_0025 = stats.t.isf(q=0.05/2, df=phi_e)
# マハラノビス距離の二乗
#vec_mean = myrmesia_remNan_int[["HL", "ML", "WL", "PrW", "PpL", "GW"]].mean() # 平均ベクトル
vec_mean = df2[['x1_E', 'x3_I', 'x2_K', 'x4', 'x1_F', 'x1_D', 'x9', 'x6', 'x1_C']].mean()
print(vec_mean)
D2 = []
for hl in hl_l:
    D2_0 = Mahala2([x1_E, x3_I, x2_K, x4, x1_F, x1_D, x9, hl, x1_C], vec_mean, Sigma)
    D2.append(D2_0)
    print(D2_0)
D2 = sp.array(D2)


interval095 = t_0025 * sp.sqrt((1/n + D2 / (n-1)) * res.scale)
print(interval095)
plt.plot(hl_l, hat_y - interval095, "--g", label="95% interval conf.")
plt.plot(hl_l, hat_y + interval095, "--g")

pred_interval095 = t_0025 * sp.sqrt((1 + 1/n + D2 / (n-1)) * res.scale)
plt.plot(hl_l, hat_y - pred_interval095, "--r", label="95% prediction conf.")
plt.plot(hl_l, hat_y + pred_interval095, "--r")


# sample x6=2 のとき
hl=2
#X = sp.array([1, hl, ml, wl, prw, ppl, gw])
X = sp.array([1, x1_E, x3_I, x2_K, x4, x1_F, x1_D, x9, hl, x1_C])
plt.scatter(hl, X.dot(res.params), color="k", )

plt.grid()

【出力】

(x1_E    0.090009
dtype: float64, x3_I    0.211285
dtype: float64, x2_K    0.043583
dtype: float64, x4    1.165385
dtype: float64, x1_F    0.040741
dtype: float64, x1_D    0.226445
dtype: float64, x9    1.035582
dtype: float64, x1_C    0.199916
dtype: float64)
x1_E    0.089118
x3_I    0.209193
x2_K    0.043152
x4      1.153846
x1_F    0.040338
x1_D    0.224203
x9      1.025328
x6      1.059099
x1_C    0.197936
dtype: float64
vec_x is [x1_E    0.090009
dtype: float64, x3_I    0.211285
dtype: float64, x2_K    0.043583
dtype: float64, x4    1.165385
dtype: float64, x1_F    0.040741
dtype: float64, x1_D    0.226445
dtype: float64, x9    1.035582
dtype: float64, 0.0, x1_C    0.199916
dtype: float64]
vec_mean is x1_E    0.089118
x3_I    0.209193
x2_K    0.043152
x4      1.153846
x1_F    0.040338
x1_D    0.224203
x9      1.025328
x6      1.059099
x1_C    0.197936
dtype: float64
vec is [[x1_E    0.000891
dtype: float64]
[x3_I    0.002092
dtype: float64]
[x2_K    0.000432
dtype: float64]
[x4    0.011538
dtype: float64]
[x1_F    0.000403
dtype: float64]
[x1_D    0.002242
dtype: float64]
[x9    0.010253
dtype: float64]
[-1.0590994371482176]
[x1_C    0.001979
dtype: float64]]
det is 3.3851611336655366e-11
mahara is [[x1_C   NaN
x1_D   NaN
x1_E   NaN
x1_F   NaN
x2_K   NaN
x3_I   NaN
x4     NaN
x9     NaN
dtype: float64]]

【デバッグ】
`python コード

ちなみに、出力の差分としましては、うまくいっているデータの時は、
関数Mahala2の中でデバッグ表示しているベクトルについて、

vec is [[-0.40788125]
[ 0.00139516]
[ 0.00359367]
[ 0.00180594]
[ 0.00321283]]

となっているのですが、NGの時は以下の様になっております。

vec is [[x1_E    0.000891
dtype: float64]
[x3_I    0.002092
dtype: float64]
[x2_K    0.000432
dtype: float64]
[x4    0.011538
dtype: float64]
[x1_F    0.000403
dtype: float64]
[x1_D    0.002242
dtype: float64]
[x9    0.010253
dtype: float64]
[-1.0186953967441772]
[x1_C    0.001979
dtype: float64]]`

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

質問への追記・修正、ベストアンサー選択の依頼

  • y_waiwai

    2018/11/05 15:32

    このままではコードが読めません。質問を編集して、<code>ボタンを押して、’’’の枠の中にコードを貼り付けてください

    キャンセル

  • _hh

    2018/11/05 17:11

    ご連絡ありがとうございます。追加情報も加え修正させて頂きました。さっぱり原因が分かりません。

    キャンセル

  • ikedas

    2018/11/07 20:31

    Mahala2の中で、invの値はどうなっていますか。

    キャンセル

回答 1

check解決した方法

0

色々と有難うございます。def Mahala2(vec_x, vec_mean, mat):
以降で、以下の様にしたところ、解決致しました。

    vec_x = np.array(vec_x, dtype='float64')

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 88.58%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る