0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

回帰係数がバラつく原因と対策:ラグ変数 × PCA × Bootstrap

Posted at

user:

目的変数Yに対して、説明変数Xをラグ変数にして、ラグ変数を用いた回帰をすると、ラグの増加に伴う回帰係数の推移が、滑らかでなくガタガタになることがよくあります。滑らかにする工夫はありませんか?

assistant:

はい、ラグ変数を用いた回帰(たとえばARXやVARなど)で、ラグに伴う回帰係数がガタガタ(不規則)に変化する問題はよくあります。この現象は、ラグが多くなることで推定の分散が大きくなることや、多重共線性、過学習、サンプルサイズの不足などが影響しています。

以下に、回帰係数を滑らかにするための実践的な工夫・方法をいくつか紹介します:


✅ 1. 正則化(リッジやLASSO)を使う

  • ラグが多くなると説明変数が増え、係数推定が不安定になります。
  • 特に リッジ回帰(L2正則化) は、係数の変動を抑える効果があり、滑らかにするのに有効です。
  • Elastic Net(L1とL2の混合)も選択肢です。
from sklearn.linear_model import Ridge
model = Ridge(alpha=1.0)
model.fit(X_lagged, y)

✅ 2. スプライン・フィルタベースの平滑化(ラグに構造を課す)

  • ラグ係数に滑らかさの構造を事前に課す方法もあります。
  • 例:係数がラグの関数(たとえば多項式やBスプライン)で表現されるようにする。
# ラグを z=1,2,... とし、lagごとに多項式展開
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
poly = PolynomialFeatures(degree=3, include_bias=False)
Z = poly.fit_transform(np.arange(1, p+1).reshape(-1, 1))
# 回帰モデルにZを使って滑らかな係数を回帰

✅ 3. ベイズ的アプローチ(滑らかな事前分布)

  • ベイズ回帰で「係数が滑らかに変化する」という事前分布(例:ガウス過程やランダムウォーク)を課すと滑らかになります。
  • 例:状態空間モデル(動的線形モデル)で、係数が時間(またはラグ)に応じて変化する。

✅ 4. ラグのグルーピング

  • ラグを個別に使うのではなく、ラグウィンドウ(例:3~5期前の平均)として使うとノイズが減ります。
df['x_lag_avg'] = df[['x_lag3', 'x_lag4', 'x_lag5']].mean(axis=1)

✅ 5. 主成分回帰(PCR)やPLS

  • ラグ変数の次元が高くなる場合、主成分回帰で低次元の成分だけ使うと安定化する。
  • 滑らかさを保証するものではないですが、変動を抑える効果があります。

✅ 6. 係数をプロットする際の平滑化

  • 係数のプロットがガタガタに見えるだけであれば、**可視化の平滑化(例:LOESSや移動平均)**も検討できます。
import pandas as pd
smoothed_coeffs = pd.Series(coefficients).rolling(window=3, center=True).mean()

🔁 まとめ

方法 滑らかさの付与方法 メリット 注意点
リッジ回帰 ペナルティ項 簡単・安定化 過剰に平滑になる可能性
スプライン回帰 ラグ構造を明示 解釈しやすい 実装少し複雑
状態空間モデル(ベイズ) 動的事前分布 柔軟性高い モデル設定が複雑
ラグの平均化 事前集約 ノイズ除去 解像度低下の恐れ
主成分回帰やPLS 次元削減 安定化 解釈が難しい

必要に応じて、具体的な実装例(sklearnやPyMC、statsmodelsなど)もお手伝いできます。どのアプローチを試してみたいですか?

user:

ありがとうございます。「ラグ変数の次元が高くなる場合、主成分回帰で低次元の成分だけ使うと安定化する」について質問です。低次元とは、どのくらいの低次元にすれば良いか、見積もる方法はありますか?

assistant:

とても良い質問です。「主成分回帰(PCR)でどのくらいの次元を残すべきか?」というのは、ラグ変数が多くなる場合の実践上の重要な論点です。以下のような観点で次元数(主成分数)を見積もる方法があります。


✅ 方法1:累積寄与率(Cumulative Explained Variance)に基づく選択

● 基本アイデア:

  • 主成分分析(PCA)の各成分が説明する**分散の割合(寄与率)**を見る。
  • 通常、80~95%の累積寄与率を超えるように次元数を選ぶ。

● 実装例(scikit-learn):

from sklearn.decomposition import PCA

# X_lagged: ラグ変数(標準化しておくと良い)
pca = PCA()
pca.fit(X_lagged)

# 累積寄与率
cumulative = pca.explained_variance_ratio_.cumsum()

# どの成分数で95%を超えるか
import numpy as np
n_components = np.argmax(cumulative >= 0.95) + 1
print(f"95%の分散を説明する主成分数: {n_components}")

✅ 方法2:交差検証(Cross-validation)で予測性能が最も高い次元数を選ぶ

● 基本アイデア:

  • 主成分数を変えながらPCRを実行し、交差検証で予測精度(MSEなど)を評価
  • MSEが最小になる次元数を採用。

● 実装例:

from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline

mse_list = []
for n in range(1, X_lagged.shape[1] + 1):
    model = make_pipeline(PCA(n_components=n), LinearRegression())
    scores = cross_val_score(model, X_lagged, y, scoring='neg_mean_squared_error', cv=5)
    mse_list.append(-scores.mean())

# MSEが最小の次元数
optimal_n = np.argmin(mse_list) + 1
print(f"最適な主成分数: {optimal_n}")

✅ 方法3:スクリープロット(Scree Plot)で「ひざ」位置を見る

  • 主成分の固有値(=説明される分散)を折れ線グラフにし、**急減から緩やかになる「ひざ」**の位置を見る。
import matplotlib.pyplot as plt

plt.plot(np.arange(1, len(pca.explained_variance_)+1), pca.explained_variance_, marker='o')
plt.xlabel("Principal Component")
plt.ylabel("Explained Variance")
plt.title("Scree Plot")
plt.grid(True)
plt.show()
  • 「エルボー法」に似た方法で、人間の判断も入りますが視覚的にわかりやすいです。

✅ 方法4:Information Criterion(AICやBIC)を使う

  • より統計的に厳密にやるなら、次元数に応じてAIC/BICを計算し最小のものを選ぶ
  • statsmodelsやPyMCでベイズ的主成分回帰を使う場合によく使われます。

🔁 まとめ:何次元を選ぶべきか?

方法 判断基準 メリット 注意点
累積寄与率 80%~95%程度 解釈が簡単 説明力に偏る可能性
交差検証 CVによる最小MSE 汎化性能に強い 計算コスト大きめ
スクリープロット 「ひざ」位置 視覚的に直感的 主観的判断が入る
AIC/BIC 情報量規準で最小化 モデル選択の理論的基盤 実装がやや面倒

ご希望であれば、上記の方法のうちどれかをあなたのデータに適用するコードを書きますので、データ構造や使用中のライブラリを教えてください。

user:

ありがとうございます。では、次の要件を満たす Python 関数を書いてください。

時刻 t における目的変数 Y_t の値を、それ以前の目的変数のラグ Y_{t-1}, Y_{t-2}, … , Y_{t-max_lag} および、ひとつの説明変数のラグ X_{t-1}, X_{t-2}, … , X_{t-max_lag} から主成分回帰する。主成分分析により写像された主成分に対する回帰係数は、 写像する前の元の変数 Y_{t-1}, Y_{t-2}, … , Y_{t-max_lag}, X_{t-1}, X_{t-2}, … , X_{t-max_lag} に分配して、元の変数に対する偏回帰係数を得る。用いる主成分の数は、累積寄与率ベースで選択する。主成分回帰でブートストラップすることで、95%信頼区間も得る。

入力
・ df : データ
・ target_var : df中の目的変数を表す行の名前
・explanatory_var : df中の説明変数を表す行の名前
・max_lag : ラグの最大値、デフォルトは 24
・contribution_rate : 主成分の数を決める時に用いる、累積寄与率の閾値。デフォルトは 0.8

手順
・説明変数・目的変数をZスコア化する。
・説明変数・目的変数から、max_lag 個のラグ変数を作る。
・説明変数のラグ変数・目的変数のラグ変数を主成分分析し、主成分を得る。主成分の数は累積寄与率ベースで決める。
・目的変数に対し、主成分を説明変数とするブートストラップ主成分回帰を行い、偏回帰係数と95%信頼区間を得る。
・主成分に対する偏回帰係数から、ラグ変数に対する偏回帰係数を得る。

出力
・ラグ変数に対する偏回帰係数と、95%信頼区間を格納した pandas dataframe

assistant:

以下の関数 bootstrap_pcr_with_confidence_intervals は、指定されたデータフレームに対して以下を行います:

  • 目的変数と説明変数のラグを構築
  • 主成分分析(PCA)で次元圧縮(累積寄与率ベースで主成分数を選択)
  • 主成分回帰(PCR)をブートストラップ付きで実行
  • 元のラグ変数に対する偏回帰係数とその95%信頼区間を出力

✅ 使用例:

result = bootstrap_pcr_with_confidence_intervals(
    df=my_df,
    target_var='Y',
    explanatory_var='X',
    max_lag=24,
    contribution_rate=0.8
)
print(result)

📋 出力:

variable coef_mean ci_lower ci_upper
Y_lag1 ... ... ...
X_lag1 ... ... ...
... ... ... ...

必要に応じて、グラフ表示や、特定のラグだけ抽出するフィルタリングなども追加できます。ご希望があれば続けて実装します!

user:

ありがとうございます。resultを入力すると、ラグの増加に伴って偏回帰係数がどのように推移するかを縦棒グラフで出力する関数を書いてください。95%信頼区間はエラーバーとして表示してください。説明変数のラグに対する図と、目的変数のラグに対する図は、別の図にしてください。

assistant:

# 再定義(コード実行環境がリセットされたため)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def plot_lag_coefficients(result_df, target_var, explanatory_var):
    # 目的変数ラグと説明変数ラグで分ける
    target_lags = result_df[result_df['variable'].str.startswith(f'{target_var}_lag')]
    explanatory_lags = result_df[result_df['variable'].str.startswith(f'{explanatory_var}_lag')]

    def plot_bar_with_error(data, title):
        lags = data['variable'].str.extract(r'_(?:lag)(\d+)').astype(int)[0]
        sorted_idx = np.argsort(lags)
        data = data.iloc[sorted_idx]
        lags_sorted = lags.iloc[sorted_idx]

        plt.figure(figsize=(10, 5))
        plt.bar(lags_sorted, data['coef_mean'], yerr=[data['coef_mean'] - data['ci_lower'], 
                                                       data['ci_upper'] - data['coef_mean']], 
                capsize=4, color='skyblue', edgecolor='black')
        plt.xlabel("Lag")
        plt.ylabel("Regression Coefficient")
        plt.title(title)
        plt.grid(True)
        plt.tight_layout()
        plt.show()

    # 目的変数ラグのプロット
    plot_bar_with_error(target_lags, f"Regression Coefficients for Lags of {target_var}")

    # 説明変数ラグのプロット
    plot_bar_with_error(explanatory_lags, f"Regression Coefficients for Lags of {explanatory_var}")

関数 plot_lag_coefficients(result_df, target_var, explanatory_var) を定義しました。この関数は次のように使えます:


✅ 使い方

plot_lag_coefficients(result, target_var='Y', explanatory_var='X')

🔍 機能の概要

  • result は、前ステップで作成した関数 bootstrap_pcr_with_confidence_intervals の出力。
  • 目的変数 Y と説明変数 X のラグに対して、それぞれ 別々の縦棒グラフ を表示。
  • 縦棒の高さが 偏回帰係数の平均エラーバーが95%信頼区間
  • ラグ番号順にソートして表示。

必要であれば、グラフの保存や色・フォント調整、横棒バージョンなども対応できます。

user:

ありがとうございます。

上記のコードをベースに、次のように書き換えてください。
・偏回帰係数の平均値が正の値の時と負の値の時で色を変える(デフォルトでは、正の値の時に steelblue 、負の値の時に tomato)
・PowerPointにそのまま貼り付けても図や文字が見えるように整える

assistant:

以下の改良された関数は、PowerPoint向けに視認性を高め、正負の偏回帰係数に応じた色分けも行います。

def plot_lag_coefficients(result_df, target_var, explanatory_var):
    import matplotlib.pyplot as plt
    import numpy as np

    # フォントと図のスタイルを設定(PowerPoint用に見やすく)
    plt.rcParams.update({
        'font.size': 14,
        'axes.labelsize': 16,
        'axes.titlesize': 18,
        'xtick.labelsize': 12,
        'ytick.labelsize': 12,
        'legend.fontsize': 12
    })

    def plot_bar_with_error(data, title):
        # ラグ番号抽出とソート
        lags = data['variable'].str.extract(r'_(?:lag)(\d+)').astype(int)[0]
        sorted_idx = np.argsort(lags)
        data = data.iloc[sorted_idx]
        lags_sorted = lags.iloc[sorted_idx].values

        # 色の設定(正:steelblue、負:tomato)
        colors = ['steelblue' if val >= 0 else 'tomato' for val in data['coef_mean']]

        # 図の描画
        fig, ax = plt.subplots(figsize=(12, 6))
        bars = ax.bar(
            lags_sorted,
            data['coef_mean'],
            yerr=[data['coef_mean'] - data['ci_lower'],
                  data['ci_upper'] - data['coef_mean']],
            capsize=5,
            color=colors,
            edgecolor='black'
        )

        ax.set_xlabel("Lag")
        ax.set_ylabel("Regression Coefficient")
        ax.set_title(title)
        ax.axhline(0, color='black', linewidth=1)
        ax.grid(True, linestyle='--', alpha=0.6)

        # 軸の調整(PowerPoint用に余白を多めに)
        ax.set_xlim([0, max(lags_sorted) + 1])
        ax.set_ylim([
            min(data['ci_lower'].min(), -0.1),
            max(data['ci_upper'].max(), 0.1)
        ])

        plt.tight_layout()
        plt.show()

    # 対象の変数に応じて別々に描画
    target_lags = result_df[result_df['variable'].str.startswith(f'{target_var}_lag')]
    explanatory_lags = result_df[result_df['variable'].str.startswith(f'{explanatory_var}_lag')]

    plot_bar_with_error(target_lags, f"Regression Coefficients for Lags of {target_var}")
    plot_bar_with_error(explanatory_lags, f"Regression Coefficients for Lags of {explanatory_var}")

✅ 特徴

  • 正の係数:steelblue負の係数:tomato(自動判別)
  • フォントサイズ大きめ、エラーバー付、網掛けグリッド、図の余白拡張
  • PowerPointへコピペしても読みやすい設計

✅ 使い方

plot_lag_coefficients(result, target_var='Y', explanatory_var='X')

必要に応じて、グラフの保存(PNG, SVG)や、ラベルの日本語対応なども可能です。ご希望があれば追加します。

user:

ありがとうございます。

先ほどの bootstrap_pcr_with_confidence_intervals 関数について、目的変数を1つ、説明変数を1つ入力する前提で設計していただきました。それを拡張して、次の2つの使い方ができるように書き換えていただけませんか?

・explanatory_var is None のとき:target_var のラグ変数だけを用いて主成分回帰する
・explanatory_var is not None のとき:これまでのbootstrap_pcr_with_confidence_intervals関数どおり、target_var のラグ変数とexplanatory_varのラグ変数を用いて主成分回帰する

assistant:

関数 bootstrap_pcr_with_confidence_intervals を次の2パターンに対応するように拡張しました:

# 再定義(環境リセットのため)

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression

def bootstrap_pcr_with_confidence_intervals(df, target_var, explanatory_var=None, max_lag=24, contribution_rate=0.8, n_bootstrap=1000, random_state=42):
    np.random.seed(random_state)

    # 1. 対象変数の抽出と標準化
    columns = [target_var] if explanatory_var is None else [target_var, explanatory_var]
    df_std = df[columns].copy()
    df_std = (df_std - df_std.mean()) / df_std.std()

    # 2. ラグ変数の作成
    lagged_data = []
    for lag in range(1, max_lag + 1):
        lagged_columns = {}
        lagged_columns[f"{target_var}_lag{lag}"] = df_std[target_var].shift(lag)
        if explanatory_var is not None:
            lagged_columns[f"{explanatory_var}_lag{lag}"] = df_std[explanatory_var].shift(lag)
        lagged_data.append(pd.DataFrame(lagged_columns))
    
    df_lagged = pd.concat([df_std] + lagged_data, axis=1).dropna()

    Y = df_lagged[target_var].values
    feature_cols = [col for col in df_lagged.columns if col not in [target_var] + ([explanatory_var] if explanatory_var else [])]
    X_lagged = df_lagged[feature_cols].values

    # 3. PCAと次元選択
    pca = PCA()
    X_pca = pca.fit_transform(X_lagged)
    cumulative_variance = np.cumsum(pca.explained_variance_ratio_)
    n_components = np.argmax(cumulative_variance >= contribution_rate) + 1
    X_pca_reduced = X_pca[:, :n_components]
    components_matrix = pca.components_[:n_components, :]

    # 4. ブートストラップPCR
    n_features = X_lagged.shape[1]
    coefs = np.zeros((n_bootstrap, n_features))

    for i in range(n_bootstrap):
        indices = np.random.choice(len(Y), len(Y), replace=True)
        X_sample = X_pca_reduced[indices]
        y_sample = Y[indices]
        model = LinearRegression().fit(X_sample, y_sample)
        coef_pca = model.coef_
        coef_original = np.dot(coef_pca, components_matrix)
        coefs[i, :] = coef_original

    # 5. 結果の統計量
    coef_mean = np.mean(coefs, axis=0)
    coef_lower = np.percentile(coefs, 2.5, axis=0)
    coef_upper = np.percentile(coefs, 97.5, axis=0)

    result_df = pd.DataFrame({
        "variable": feature_cols,
        "coef_mean": coef_mean,
        "ci_lower": coef_lower,
        "ci_upper": coef_upper
    })

    return result_df

✅ 対応パターン

  1. explanatory_var is None の場合
     → 目的変数 target_var のラグのみを使って主成分回帰

  2. explanatory_var is not None の場合
     → target_varexplanatory_var のラグを両方使って主成分回帰(従来の動作)


✅ 使用例

# 自己ラグだけで回帰
result1 = bootstrap_pcr_with_confidence_intervals(df, target_var='Y', explanatory_var=None)

# 自己ラグ+説明変数のラグで回帰
result2 = bootstrap_pcr_with_confidence_intervals(df, target_var='Y', explanatory_var='X')

どちらのパターンでも同様に plot_lag_coefficients(result_df, ...) を使って可視化できます。

他に多変量拡張(説明変数が複数)などの要望があれば、それにも対応可能です。

user:

ありがとうございました。今回はこの辺りで大丈夫です。ここまでの議論をブログに記したいので、タイトル案をいくつか考えていただけますか?

assistant:

もちろんです!今回の内容は、ラグ変数の主成分回帰による係数の安定化とブートストラップ信頼区間の可視化という高度で実践的なテーマでした。それを踏まえたブログタイトル案を、目的や読者層別にご提案します:


🔬 技術的に正確&専門家向け

  1. 「ラグ変数の主成分回帰で係数を安定化:ブートストラップで信頼区間も可視化」
  2. 「自己・他己ラグの次元圧縮と感度分析:主成分回帰 × ブートストラップ入門」
  3. 「ラグ系列のガタつきを抑える:PCA+PCRによる信頼区間付き回帰分析」

🧑‍🏫 学習者・実務家向け(少し親しみやすく)

  1. 「ラグ変数が多すぎて不安定?主成分回帰+ブートストラップで滑らかに」
  2. 「時系列回帰がガタガタする?主成分で圧縮して回帰係数を安定させる方法」
  3. 「回帰係数がバラつく原因と対策:ラグ変数 × PCA × Bootstrap」

✍️ 読み物・ブログらしい柔らかめの表現

  1. 「回帰係数が暴れる…!そんなときの主成分回帰+ブートストラップ」
  2. 「ラグのせいでモデルが不安定?主成分で整理してから回帰しよう」
  3. 「PCAで滑らかに!ラグ変数だらけの回帰を整理するテクニック」

ご希望に応じて、SEO対策寄り、シリーズ化、サブタイトル付きなどもご提案できます。気に入ったスタイルがあれば、それに寄せたリライト案も可能です。

user:

「回帰係数がバラつく原因と対策:ラグ変数 × PCA × Bootstrap」で行きたいと思います。最後に、ここまでの内容の詳しいまとめを書いてください。

assistant:

素晴らしいタイトルのご選択です!
以下に、ブログ記事「回帰係数がバラつく原因と対策:ラグ変数 × PCA × Bootstrap」に相応しい、**詳しいまとめ(記事の締めや導入にも使える構成)**を記します。


🧠 回帰係数がバラつく原因と対策:ラグ変数 × PCA × Bootstrap — 詳しいまとめ

🔍 なぜラグ変数の回帰係数はガタガタになるのか?

時系列データの解析では、目的変数や説明変数の**過去の値(ラグ)**を使って予測モデルを構築することが一般的です。しかし、ラグを多く含めると、以下のような問題に直面します:

  • 多重共線性:ラグ変数同士が似た情報を持ち、回帰係数が不安定になる。
  • 高次元化による過学習:モデルがデータに過剰に適合し、係数がバラつく。
  • データ不足:サンプル数に比して変数が多すぎると、推定の信頼性が低下。

これにより、ラグごとの回帰係数が滑らかに変化せず、ガタガタとばらついた推移を示すことがよくあります。


🛠️ 対策1:PCAでラグ変数を次元圧縮する

この問題への有効な対策の一つが、主成分分析(PCA)によるラグ変数の次元圧縮です。

  • 目的:ラグ変数群の情報を保持しつつ次元を縮小し、回帰を安定化。
  • 方法:ラグ変数をPCAで主成分に変換し、主要な成分(累積寄与率80%など)だけを使う。
  • 結果:不要なノイズを除き、滑らかで意味のある係数推定が可能になる。

🛠️ 対策2:ブートストラップで信頼区間を算出する

係数推定の不確実性を定量的に把握するには、ブートストラップ法が有効です。

  • データを繰り返し再サンプリングして主成分回帰を実行。
  • 得られた係数の分布から、**95%信頼区間(CI)**を導出。
  • 結果の可視化により、係数が有意か・不安定かを明確に把握できる。

🔁 処理の流れ(Pythonで実装)

今回の分析では以下のステップをPythonで実装しました:

  1. ラグ変数の生成(目的変数・説明変数)
  2. Zスコア標準化
  3. PCAによる主成分抽出と主成分回帰
  4. ブートストラップで回帰を繰り返し、信頼区間を算出
  5. 主成分空間から元のラグ変数空間へ係数を逆写像
  6. 係数とCIを棒グラフで可視化(正負で色分け)

📈 成果物:滑らかな係数推定と可視化

  • 主成分を用いた回帰により、ばらつきの少ない安定した係数を推定
  • 信頼区間の可視化により、どのラグが重要かが明確に
  • 正の影響と負の影響を色分け表示し、解釈性を向上

📝 まとめ:時系列回帰を“構造化”する技術

時系列回帰における「ガタガタ問題」は、単なる見た目の問題ではなく、モデルの信頼性・頑健性に直結する課題です。

  • PCAによる次元圧縮は、ラグ変数の構造的冗長性を整理する手法。
  • ブートストラップは、係数推定の不確実性を可視化する現実的な方法。
  • 両者を組み合わせることで、信頼性の高い感度分析が可能になります。

🔗 今後の展望

  • 説明変数を複数に拡張する(多変量主成分回帰)
  • 動的主成分回帰(時系列に特化したPCA)
  • ベイズ主成分回帰による事前情報の組み込み

記事の最後には、コード全文やダウンロードリンク、使用データの紹介を載せておくと読者の実践に繋がります。必要であればそのまとめもお手伝いできます!

user:

ありがとうございました!

assistant:

どういたしまして!😊
とても深くて実践的なテーマに取り組まれていて素晴らしいです。
ブログ執筆、応援しています!また何かあればいつでもお声かけくださいね。

🧠📈✍️「ラグ × PCA × Bootstrap」の次の一歩も、いつでもご一緒します。
お疲れさまでした!

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?