ポケモンと線形回帰


##  はじめに

このプログラムは、ポケットモンスター（ポケモン）シリーズのゲーム作品に関するデータ（発売年・売上本数・リメイクの有無など）を対象として、以下のような総合的なデータサイエンス手法を実践的に学ぶことを目的としています


##  参考リンクまとめ

https://ja.wikipedia.org/wiki/%E3%83%9D%E3%82%B1%E3%83%83%E3%83%88%E3%83%A2%E3%83%B3%E3%82%B9%E3%82%BF%E3%83%BC_(%E3%82%B2%E3%83%BC%E3%83%A0%E3%82%B7%E3%83%AA%E3%83%BC%E3%82%BA)



##  Pythonコード


```python
# --- ライブラリのインポート / Import Libraries ---
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import boxcox, shapiro, spearmanr, kendalltau, zscore
from sklearn.linear_model import Lasso, Ridge, BayesianRidge
from sklearn.svm import SVR
from sklearn.neighbors import KNeighborsRegressor
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
import statsmodels.api as sm
import statsmodels.formula.api as smf
import xgboost as xgb

# --- データ準備 / Prepare Pokémon Sales Data ---
df = pd.DataFrame({
    'Title': ['Red/Green', 'Blue', 'Pikachu', 'Gold/Silver', 'Crystal', 'Ruby/Sapphire', 'FireRed/LeafGreen',
              'Emerald', 'Diamond/Pearl', 'Platinum', 'HeartGold/SoulSilver', 'Black/White',
              'Black2/White2', 'X/Y', 'Omega Ruby/Alpha Sapphire', 'Sun/Moon', 'Ultra Sun/Ultra Moon',
              "Let's Go!", 'Sword/Shield', 'BD/SP', 'Legends Arceus', 'Scarlet/Violet'],
    'Year': [1996, 1996, 1998, 1999, 2000, 2002, 2004, 2004, 2006, 2008, 2009, 2010, 2012, 2013, 2014,
             2016, 2017, 2018, 2019, 2021, 2022, 2022],
    'Sales': [46.02, 2.01, 14.64, 29.49, 6.39, 16.14, 12.00, 6.80, 17.67, 7.60, 12.72, 15.64, 8.52,
              16.76, 14.63, 16.33, 9.23, 15.07, 26.60, 14.92, 14.83, 26.38],
    'Remake': [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0,
               0, 1, 0, 0, 1, 0, 1, 0, 0]
})

# --- 外れ値除去（Zスコアベース）/ Outlier Removal ---
df['z_score'] = zscore(df['Sales'])               # zスコアを計算 / Compute z-score
df = df[df['z_score'].abs() < 2.5].reset_index(drop=True)  # z > 2.5 を除去 / Remove outliers
df.drop(columns='z_score', inplace=True)          # z_score列を削除 / Drop helper column

# --- データ標準化 / Standardize Features ---
scaler = StandardScaler()
X = scaler.fit_transform(df[['Year', 'Remake']])
y = df['Sales'].values

# --- Box-Cox変換 / Box-Cox Transformation ---
y_bc, lmbda = boxcox(y + 1)  # avoid zero
stat, p_value = shapiro(y_bc)
print(f'Shapiro-Wilk test (Box-Cox Sales): stat={stat:.4f}, p={p_value:.4f}')

# --- ANOVA（リメイク効果）/ ANOVA for Remake Effect ---
df['Remake'] = df['Remake'].astype('category')
model = smf.ols('Sales ~ Remake', data=df).fit()
anova_result = sm.stats.anova_lm(model)
print(anova_result)

# --- 非線形相関 / Spearman & Kendall ---
corr_spearman, _ = spearmanr(df['Year'], df['Sales'])
corr_kendall, _ = kendalltau(df['Year'], df['Sales'])
print(f"Spearman: {corr_spearman:.3f}, Kendall: {corr_kendall:.3f}")

# --- 線形回帰（YearとRemake）/ Linear Regression (OLS) ---
X_lr = sm.add_constant(df[['Year', 'Remake']])  # 定数項追加 / Add constant
lr_model = sm.OLS(df['Sales'], X_lr).fit()
print(lr_model.summary())

# --- 線形回帰プロット / Regression Plot ---
plt.figure(figsize=(8, 6))
sns.regplot(data=df, x='Year', y='Sales', scatter_kws={'s': 50})
plt.title('Linear Regression: Year vs Sales')
plt.xlabel('Year')
plt.ylabel('Sales (in millions)')
plt.grid(True)
plt.show()

# --- モデル比較 / Regression Model Comparison ---
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)

models = {
    'Lasso': Lasso(alpha=0.1),
    'Ridge': Ridge(alpha=1.0),
    'BayesianRidge': BayesianRidge(),
    'SVR': SVR(),
    'KNN': KNeighborsRegressor(),
    'XGBoost': xgb.XGBRegressor(objective='reg:squarederror')
}

model_results = {}
for name, model in models.items():
    model.fit(X_train, y_train)
    preds = model.predict(X_test)
    rmse = np.sqrt(mean_squared_error(y_test, preds))
    r2 = r2_score(y_test, preds)
    model_results[name] = {'RMSE': rmse, 'R2': r2}
    print(f"{name}: RMSE={rmse:.2f}, R2={r2:.3f}")

# --- PCA可視化 / PCA Visualization ---
pca = PCA(n_components=2)
pca_result = pca.fit_transform(X)
df['PCA1'] = pca_result[:, 0]
df['PCA2'] = pca_result[:, 1]

plt.figure(figsize=(8, 6))
sns.scatterplot(data=df, x='PCA1', y='PCA2', hue='Remake', palette='Set1')
plt.title('PCA Visualization of Pokémon Sales')
plt.xlabel('PCA Component 1')
plt.ylabel('PCA Component 2')
plt.grid(True)
plt.show()

# --- t-SNE可視化 / t-SNE Visualization ---
tsne_result = TSNE(n_components=2, perplexity=5, random_state=42).fit_transform(X)
df['tSNE1'] = tsne_result[:, 0]
df['tSNE2'] = tsne_result[:, 1]

plt.figure(figsize=(8, 6))
sns.scatterplot(data=df, x='tSNE1', y='tSNE2', hue='Remake', palette='Set2')
plt.title('t-SNE Visualization of Pokémon Sales')
plt.xlabel('t-SNE Dimension 1')
plt.ylabel('t-SNE Dimension 2')
plt.grid(True)
plt.show()
# --- ライブラリ / Libraries ---
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import boxcox, shapiro, spearmanr, kendalltau, zscore
from sklearn.linear_model import LinearRegression, Lasso, Ridge, BayesianRidge
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.svm import SVR
from sklearn.neighbors import KNeighborsRegressor
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
from sklearn.model_selection import train_test_split, cross_val_score, KFold
from sklearn.metrics import mean_squared_error, r2_score
import statsmodels.api as sm
import statsmodels.formula.api as smf
import xgboost as xgb

# --- データ定義 / Pokémon Sales Data ---
df = pd.DataFrame({
    'Title': ['Red/Green', 'Blue', 'Pikachu', 'Gold/Silver', 'Crystal', 'Ruby/Sapphire', 'FireRed/LeafGreen',
              'Emerald', 'Diamond/Pearl', 'Platinum', 'HeartGold/SoulSilver', 'Black/White',
              'Black2/White2', 'X/Y', 'Omega Ruby/Alpha Sapphire', 'Sun/Moon', 'Ultra Sun/Ultra Moon',
              "Let's Go!", 'Sword/Shield', 'BD/SP', 'Legends Arceus', 'Scarlet/Violet'],
    'Year': [1996, 1996, 1998, 1999, 2000, 2002, 2004, 2004, 2006, 2008, 2009, 2010, 2012, 2013, 2014,
             2016, 2017, 2018, 2019, 2021, 2022, 2022],
    'Sales': [46.02, 2.01, 14.64, 29.49, 6.39, 16.14, 12.00, 6.80, 17.67, 7.60, 12.72, 15.64, 8.52,
              16.76, 14.63, 16.33, 9.23, 15.07, 26.60, 14.92, 14.83, 26.38],
    'Remake': [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0,
               0, 1, 0, 0, 1, 0, 1, 0, 0]
})

# --- 外れ値除去 / Outlier Removal ---
df['z'] = zscore(df['Sales'])
df = df[df['z'].abs() < 2.5].drop(columns='z').reset_index(drop=True)

# --- 標準化 / Standardization ---
scaler = StandardScaler()
X = scaler.fit_transform(df[['Year', 'Remake']])
y = df['Sales'].values

# --- Shapiro-Wilk + Q-Qプロット ---
y_bc, _ = boxcox(y + 1)
stat, p = shapiro(y_bc)
print(f'Shapiro-Wilk: stat={stat:.4f}, p={p:.4f}')
sm.qqplot(y_bc, line='s')
plt.title('Q-Q Plot of Box-Cox Transformed Sales')
plt.grid(True)
plt.show()

# --- リメイク別平均 + 誤差棒付きプロット ---
group_stats = df.groupby('Remake')['Sales'].agg(['mean', 'sem']).reset_index()
plt.errorbar(group_stats['Remake'].astype(str), group_stats['mean'], yerr=group_stats['sem'], fmt='o', capsize=5)
plt.title('Sales Mean by Remake (w/ SEM)')
plt.xlabel('Remake')
plt.ylabel('Sales')
plt.grid(True)
plt.show()

# --- 分布比較（Boxplot / Violin） ---
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
sns.boxplot(data=df, x='Remake', y='Sales')
plt.title('Boxplot: Sales by Remake')
plt.subplot(1, 2, 2)
sns.violinplot(data=df, x='Remake', y='Sales')
plt.title('Violinplot: Sales by Remake')
plt.tight_layout()
plt.show()

# --- 線形回帰 + 残差プロット ---
X_lr = sm.add_constant(df[['Year', 'Remake']])
lr_model = sm.OLS(df['Sales'], X_lr).fit()
print(lr_model.summary())

sns.residplot(x=lr_model.fittedvalues, y=lr_model.resid, lowess=True, line_kws={'color': 'red'})
plt.axhline(0, linestyle='--', color='gray')
plt.title('Residual Plot')
plt.xlabel('Fitted Values')
plt.ylabel('Residuals')
plt.grid(True)
plt.show()

# --- 多項式回帰 ---
poly = PolynomialFeatures(degree=2)
X_poly = poly.fit_transform(df[['Year']])
lr_poly = LinearRegression().fit(X_poly, df['Sales'])
years_sorted = np.linspace(df['Year'].min(), df['Year'].max(), 100)
X_pred = poly.transform(years_sorted.reshape(-1, 1))
sales_pred = lr_poly.predict(X_pred)
plt.scatter(df['Year'], df['Sales'], label='Data')
plt.plot(years_sorted, sales_pred, color='red', label='Poly Fit')
plt.title('Polynomial Regression')
plt.xlabel('Year')
plt.ylabel('Sales')
plt.legend()
plt.grid(True)
plt.show()

# --- モデル交差検証スコア（K-Fold） ---
model_cv = Ridge()
scores = cross_val_score(model_cv, X, y, cv=KFold(n_splits=5), scoring='r2')
print(f"Ridge CV R² scores: {scores}")
print(f"Mean CV R²: {scores.mean():.3f}")

# --- 年代別トレンドプロット ---
plt.figure(figsize=(8, 5))
sns.lineplot(data=df, x='Year', y='Sales', marker='o')
sns.regplot(data=df, x='Year', y='Sales', scatter=False, color='red', label='Trend')
plt.title('Sales Over Time with Trend Line')
plt.xlabel('Year')
plt.ylabel('Sales')
plt.grid(True)
plt.legend()
plt.show()

# --- 売上ヒートマップ（相関行列） ---
corr = df[['Year', 'Sales', 'Remake']].corr()
sns.heatmap(corr, annot=True, cmap='coolwarm', fmt='.2f')
plt.title('Correlation Heatmap')
plt.show()

```