FIFA18のデータからサンフレッチェの選手の価格を予想してみた

データ分析のコンペティションで有名なKaggleはさまざまなデータが公開されています。例えば EA SPORTSからFIFA18のデータが公開されています。そこで今回はこのデータを使用して何か分析をしてみようと思います。（私はFIFA18をしないですが。。。）

https://www.kaggle.com/thec03u5/fifa-18-demo-player-dataset

FIFA18のデータは海外の有名選手やJリーグの選手含めて登録されています。また、データの項目として選手の能力や価格があります。そこでサンフレッチェの選手が他の選手の能力と比較してどのぐらいの価格が妥当か予測してみました。分析するデータはFIFA18に登録されているデータを元にしてますので、結果が悪くてもEA SPORTSに言ってくださいｗ。また、FIFA18といいつつも選手のデータは2017年のデータですので今年度のデータはありません。

FIFA18に登録されているサンフレッチェ広島の価格について

フィールドプレーヤーはごらんのとおりです。単位はEUROです。DF陣が低いですね。そして、皆川が思ったより高いのが

	name	eur_value
4689	Felipe Silva	2500000
4877	T. Aoyama	1500000
5053	Y. Kashiwa	1700000
5397	Patric	1800000
5436	K. Shibasaki	1400000
5965	Anderson Lopes	1400000
6269	D. Niwa	700000
6477	M. Mikić	350000
7871	M. Kudo	900000
8072	K. Chiba	450000
8961	H. Mizumoto	425000
10110	K. Morisaki	70000
10651	T. Miyayoshi	725000
10877	N. Burns	550000
11074	S. Inagaki	575000
12756	K. Mukuhara	300000
13090	Y. Chajima	375000
13429	Y. Minagawa	375000
14240	T. Morishima	375000
14859	T. Marutani	190000
15740	S. Sasaki	110000
15769	Y. Nogami	130000
16312	S. Takahashi	140000
16706	O. Iyoha	180000
17938	T. Matsumoto	60000

GKはこちらの通りです。

	name	eur_value
10100	T. Hayashi	220000
14571	R. Hironaga	160000
15826	H. Nakabayashi	80000

データの解析

まず、データを読み込みます。気を付けたい点は文字コードをUTF-8で読み込むことです。そうしないとクロアチア人のデータが読み込めません。

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from scipy.stats import norm
from sklearn.preprocessing import StandardScaler
from scipy import stats
%matplotlib inline
base = pd.read_csv('complete.csv', encoding='utf_8')

サンフレとそれ以外のデータに分けます。また、フィールドプレーヤーとGKに分けます。

base2 = base.copy()
sanf_base = base2[base2['club'].isin(['Sanfrecce Hiroshima'])]
base3 = base2[(base2['club'] != 'Sanfrecce Hiroshima')]
sanf_f = sanf_base[sanf_base['gk'].isnull()]
sanf_g = sanf_base[sanf_base['cf'].isnull()]
base_f = base3[base3['gk'].isnull()]
base_g = base3[base3['cf'].isnull()]

サンフレの選手はこんな感じです。

	ID	name	full_name	club	special	age	league
4689	192982	Felipe Silva	Felipe de Oliveira Silva	Sanfrecce Hiroshima	1659	27	Japanese J1 League
4877	222988	T. Aoyama	Toshihiro Aoyama	Sanfrecce Hiroshima	1879	31	Japanese J1 League
5053	232512	Y. Kashiwa	Yoshifumi Kashiwa	Sanfrecce Hiroshima	1714	29	Japanese J1 League
5397	202401	Patric	Anderson Patric Aguiar Oliveira	Sanfrecce Hiroshima	1670	29	Japanese J1 League
5436	232619	K. Shibasaki	Kosei Shibasaki	Sanfrecce Hiroshima	1812	32	Japanese J1 League
5965	229949	Anderson Lopes	Anderson José Lopes de Souza	Sanfrecce Hiroshima	1722	23	Japanese J1 League
6269	232598	D. Niwa	Daiki Niwa	Sanfrecce Hiroshima	1459	31	Japanese J1 League
6477	141265	M. Mikić	Mihael Mikić	Sanfrecce Hiroshima	1792	37	Japanese J1 League
7871	231972	M. Kudo	Masato Kudo	Sanfrecce Hiroshima	1603	27	Japanese J1 League
8072	162388	K. Chiba	Kazuhiko Chiba	Sanfrecce Hiroshima	1541	32	Japanese J1 League
8961	197422	H. Mizumoto	Hiroki Mizumoto	Sanfrecce Hiroshima	1568	31	Japanese J1 League
10110	232509	K. Morisaki	Kazuyuki Morisaki	Sanfrecce Hiroshima	1557	36	Japanese J1 League
10651	232896	T. Miyayoshi	Takumi Miyayoshi	Sanfrecce Hiroshima	1522	24	Japanese J1 League
10877	181491	N. Burns	Nathan Burns	Sanfrecce Hiroshima	1777	29	Japanese J1 League
11074	232991	S. Inagaki	Sho Inagaki	Sanfrecce Hiroshima	1747	25	Japanese J1 League
12756	237788	K. Mukuhara	Kenta Mukuhara	Sanfrecce Hiroshima	1478	27	Japanese J1 League
13090	232513	Y. Chajima	Yusuke Chajima	Sanfrecce Hiroshima	1649	25	Japanese J1 League
13429	232614	Y. Minagawa	Yusuke Minagawa	Sanfrecce Hiroshima	1400	25	Japanese J1 League
14240	232898	T. Morishima	Tsukasa Morishima	Sanfrecce Hiroshima	1453	20	Japanese J1 League
14859	232613	T. Marutani	Takuya Marutani	Sanfrecce Hiroshima	1509	28	Japanese J1 League
15740	232511	S. Sasaki	Sho Sasaki	Sanfrecce Hiroshima	1402	27	Japanese J1 League
15769	235089	Y. Nogami	Yuki Nogami	Sanfrecce Hiroshima	1341	26	Japanese J1 League
16312	232621	S. Takahashi	Soya Takahashi	Sanfrecce Hiroshima	1458	21	Japanese J1 League
16706	236761	O. Iyoha	Osamu Henry Iyoha	Sanfrecce Hiroshima	1306	19	Japanese J1 League
17938	237535	T. Matsumoto	Taishi Matsumoto	Sanfrecce Hiroshima	1345	18	Japanese J1 League

ちなみに他のデータはこんな感じです。なぜメッシよりR.ロナウドが上なのか分りません。

	ID	name	full_name	club	special	age
0	20801	Cristiano Ronaldo	C. Ronaldo dos Santos Aveiro	Real Madrid CF	2228	32	Spanish Primera División
1	158023	L. Messi	Lionel Messi	FC Barcelona	2158	30	Spanish Primera División
2	190871	Neymar	Neymar da Silva Santos Jr.	Paris Saint-Germain	2100	25	French Ligue 1
3	176580	L. Suárez	Luis Suárez	FC Barcelona	2291	30	Spanish Primera División
5	188545	R. Lewandowski	Robert Lewandowski	FC Bayern Munich	2146	28	German Bundesliga

選手の価格帯のを調べます。明らかに正規分布にはなってないですね。。こちらはフィールドプレーヤーですがGKも同じような分布になりました。

sns.distplot(base_f['eur_value']);

次に欠損値について調べます。フィールドプレーヤーで調べたのでGKが100%となっています。eur_release_clauseは今回使用しないので列ごと削除します。league、clubのデータが0ですが、これは未所属であり価格が0であり、データ分析に有用とはならないので削除します。また価格が0の選手もすべて削除します。

base_f_na = (base_f.isnull().sum() / len(base_f)) * 100
base_f_na = base_f_na.drop(base_f_na[base_f_na == 0].index).sort_values(ascending=False)[:30]
missing_data = pd.DataFrame({'Missing Ratio' :base_f_na})
missing_data.head(20)

	Missing Ratio
gk	100.000000
eur_release_clause	8.590419
league	1.348131
club	1.348131

分析用にデータを分けます。

base_f =base_f.drop(['eur_release_clause','gk'], axis=1)
base_f = base_f.dropna(how='any')
base_f = base_f.drop(base_f[(base_f['eur_value']==0)].index)
base_g = base_g.drop(base_g[(base_g['eur_value']==0)].index)                          
base_g =base_g.drop(['eur_release_clause','rs','rw','rf','ram','rcm','rm','rdm','rcb','rb','rwb','st','lw','cf','cam','cm','lm','cdm','cb','lb','lwb','ls','lf','lam','lcm','ldm','lcb'], axis=1)
base_g = base_g.dropna(how='any')
train_f_ID =base_f.loc[:,['ID','name','full_name']]
train_g_ID =base_g.loc[:,['ID','name','full_name']]
test_f_ID =sanf_f.loc[:,['ID','name','full_name']]
test_g_ID =sanf_g.loc[:,['ID','name','full_name']]
y_f_train = base_f['eur_value']
y_g_train = base_g['eur_value']
y_f_test = sanf_f['eur_value']
y_g_test = sanf_g['eur_value']
X_f_train =base_f.drop(['ID','name','full_name','eur_value','eur_wage'], axis=1)
X_g_train =base_g.drop(['ID','name','full_name','eur_value','eur_wage'], axis=1)
X_f_test =sanf_f.drop(['ID','name','full_name','eur_value','eur_wage','eur_release_clause','gk'], axis=1)
X_g_test =sanf_g.drop(['ID','name','full_name','eur_value','eur_wage','eur_release_clause','rs','rw','rf','ram','rcm','rm','rdm','rcb','rb','rwb','st','lw','cf','cam','cm','lm','cdm','cb','lb','lwb','ls','lf','lam','lcm','ldm','lcb'], axis=1)

文字列、ブール型になっている値を LabelEncoder使って実数に変換します。

from sklearn.preprocessing import LabelEncoder

for i in range(X_f_train.shape[1]):
    if X_f_train.iloc[:,i].dtypes == object:
        lbl = LabelEncoder()
        lbl.fit(list(X_f_train.iloc[:,i].values) + list(X_f_test.iloc[:,i].values))
        X_f_train.iloc[:,i] = lbl.transform(list(X_f_train.iloc[:,i].values))
        X_f_test.iloc[:,i] = lbl.transform(list(X_f_test.iloc[:,i].values))

for i in range(X_f_train.shape[1]):
    if X_f_train.iloc[:,i].dtypes == bool:
        lbl = LabelEncoder()
        lbl.fit(list(X_f_train.iloc[:,i].values) + list(X_f_test.iloc[:,i].values))
        X_f_train.iloc[:,i] = lbl.transform(list(X_f_train.iloc[:,i].values))
        X_f_test.iloc[:,i] = lbl.transform(list(X_f_test.iloc[:,i].values))

for i in range(X_g_train.shape[1]):
    if X_g_train.iloc[:,i].dtypes == bool:
        lbl = LabelEncoder()
        lbl.fit(list(X_g_train.iloc[:,i].values) + list(X_g_test.iloc[:,i].values))
        X_g_train.iloc[:,i] = lbl.transform(list(X_g_train.iloc[:,i].values))
        X_g_test.iloc[:,i] = lbl.transform(list(X_g_test.iloc[:,i].values))

for i in range(X_g_train.shape[1]):
    if X_g_train.iloc[:,i].dtypes == object:
        lbl = LabelEncoder()
        lbl.fit(list(X_g_train.iloc[:,i].values) + list(X_g_test.iloc[:,i].values))
        X_g_train.iloc[:,i] = lbl.transform(list(X_g_train.iloc[:,i].values))
        X_g_test.iloc[:,i] = lbl.transform(list(X_g_test.iloc[:,i].values))

価格をlogを使用して擬似的に正規分布にします。

y_f_train = np.log(y_f_train)
y_g_train = np.log(y_g_train)

重要な項目について累乗します。

X_f_train['overall2'] = X_f_train['overall'] ** 2
X_f_train['potential2'] = X_f_train['potential'] ** 2
X_f_test['overall2'] = X_f_test['overall'] ** 2
X_f_test['potential2'] = X_f_test['potential'] ** 2
X_g_train['gk2'] = X_g_train['gk'] ** 2
X_g_train['overall2'] = X_g_train['overall'] ** 2
X_g_train['potential2'] = X_g_train['potential'] ** 2
X_g_test['gk2'] = X_g_test['gk'] ** 2
X_g_test['overall2'] = X_g_test['overall'] ** 2
X_g_test['potential2'] = X_g_test['potential'] ** 2

あとはXGBoostとLightGBMを使用して価格を予測します。こちらはフィールドプレーヤーですがGKもおんなじような感じです。

reg = xgb.XGBRegressor(max_depth=2,n_estimators=300,min_child_weight=1,gamma=0.1,subsample=0.2,colsample_bytree=0.4,learning_rate=0.1,reg_alpha=0.0001)
reg.fit(X_f_train, y_f_train)

import lightgbm as lgb
lgb_train = lgb.Dataset(X_f_train, y_f_train)
lgb_eval  = lgb.Dataset(X_f_test, y_f_test, reference=lgb_train)
params = {
    'task': 'train',
    'boosting_type': 'gbdt',
    'objective': 'regression',
    'metric': {'l2'},
    'num_leaves': 128,
    'learning_rate': 0.01,
    'num_iterations':100,
    'feature_fraction': 0.38,
    'bagging_fraction': 0.68,
    'bagging_freq': 5,
    'verbose': 0
}


gbm = lgb.train(params,
                lgb_train,
                num_boost_round=1500,
                valid_sets=lgb_eval,
                early_stopping_rounds=15)

 
X_f_train2 = pd.DataFrame( {'XGB': reg.predict(X_f_train),
     'LGB': gbm.predict(np.array(X_f_train), num_iteration=gbm.best_iteration),
    })

from sklearn import linear_model

reg_all = linear_model.LinearRegression()
reg_all.fit(X_f_train2, y_f_train)

X_f_test2 = pd.DataFrame( {'XGB': reg.predict(X_f_test),
     'LGB': gbm.predict(np.array(X_f_test), num_iteration=gbm.best_iteration),
    })

predictions = np.exp(reg_all.predict(X_f_test2))
 
submission = pd.DataFrame({ 'ID': test_f_ID['ID'],
                           'name' : test_f_ID['name'],
                           'full_name' : test_f_ID['full_name'],
                            'eur_value<200b>' : predictions})
submission

結果はこのようになりました。皆川が高いのが違和感あるなあ。

	ID	eur_value?	full_name	name
0	192982	2376986	Felipe de Oliveira Silva	Felipe Silva
1	222988	1324306	Toshihiro Aoyama	T. Aoyama
2	232512	1634586	Yoshifumi Kashiwa	Y. Kashiwa
3	202401	1693254	Anderson Patric Aguiar Oliveira	Patric
4	232619	1366923	Kosei Shibasaki	K. Shibasaki
5	229949	1548074	Anderson Jos? Lopes de Souza	Anderson Lopes
6	232598	685995	Daiki Niwa	D. Niwa
7	141265	329592	Mihael Miki?	M. Miki?
8	231972	842412	Masato Kudo	M. Kudo
9	162388	463991	Kazuhiko Chiba	K. Chiba
10	197422	439531	Hiroki Mizumoto	H. Mizumoto
11	232509	99654	Kazuyuki Morisaki	K. Morisaki
12	232896	694564	Takumi Miyayoshi	T. Miyayoshi
13	181491	537812	Nathan Burns	N. Burns
14	232991	641852	Sho Inagaki	S. Inagaki
15	237788	311987	Kenta Mukuhara	K. Mukuhara
16	232513	356286	Yusuke Chajima	Y. Chajima
17	232614	374843	Yusuke Minagawa	Y. Minagawa
18	232898	367593	Tsukasa Morishima	T. Morishima
19	232613	182939	Takuya Marutani	T. Marutani
20	232511	108015	Sho Sasaki	S. Sasaki
21	235089	131221	Yuki Nogami	Y. Nogami
22	232621	152696	Soya Takahashi	S. Takahashi
23	236761	198621	Osamu Henry Iyoha	O. Iyoha
24	237535	65290	Taishi Matsumoto	T. Matsumoto

GKはこのような結果となりました。廣永のほうが高いのが意外ですね。

	ID	eur_value​	full_name	name
10100	232505	233405	Takuto Hayashi	T. Hayashi
14571	232507	169154	Ryotaro Hironaga	R. Hironaga
15826	237426	78575	Hirotsugu Nakabayashi	H. Nakabayashi

今期のサンフレッチェ広島の観客動員数を機械学習で予測した 結果編

今期もJリーグが開幕しました。

私が応援するサンフレッチェもタイの英雄ティーラシンのゴールを守り切って勝つことができました。 まあ課題はありますが、このまま一つずつ勝っていけばいいのではないでしょうか？

では今日の観客動員がどうだったか結果発表です。こちらで今日の観客動員を17645人を予測しましたが 実際どうだったでしょうか？ masser199.hatenablog.com

結果は17,026人でした。619人差でしたが、なかなかいい線行ってたのではないでしょうか？

KaggleのHouse Pricesに挑戦 その１

今回kaggleの初心者向き課題の一つであるHousingPriceに挑戦しました。とはいいつつも項目数が70以上もあり、参加している人もレベルがかなり高いので前回挑戦したタイタニックに比べればはるかに難しいです。

HousingPriceとは

米国アイオワ州のエイムズという都市の物件価格を予測する問題となります。データは、”築年数”、”設備”、”広さ”、”エリア”、”ガレージに入る車の数”など79個の項目数および1460戸の学習データが与えられます。そのデータをもとに、残りの1459戸の家の価格を予測します。

はじめに

今回課題を取り組むにあたって以下のkernel、ブログを参考にしました。

Stacked Regressions : Top 4% on LeaderBoard | Kaggle

こちらのリンクでは欠損値の埋め方等を参考にしました。

KaggleのHousePrices問題を決定木系のアンサンブルで解く - mirandora.commirandora.com

機械学習3ヶ月の初心者がKaggleHousingPriceに挑戦した(Top29%) 前半 - Qiita

Kaggleの練習問題（Regression）を解いてKagglerになる - Qiita

こちらのサイトではどのアルゴリズムを採用するか、また特徴量の考え方、抽出について非常に参考になりました。

価格の分析

まずそもそもアイオワのエイムズに関して知識もないし、価格も全く分からないので、まずは価格帯についてどうなのか調べました。

train = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')
train['SalePrice'].describe()

count      1460.000000
mean     180921.195890
std       79442.502883
min       34900.000000
25%      129975.000000
50%      163000.000000
75%      214000.000000
max      755000.000000
Name: SalePrice, dtype: float64

sns.distplot(train['SalePrice']);

正規分布になってないので、logを使用して一時的に正規分布になるように変換します。

train_ID = train['Id']
test_ID = test['Id']
y_train = train['SalePrice']
X_train = train.drop(['Id','SalePrice'], axis=1)
X_test = test.drop('Id', axis=1)
ntrain = X_train.shape[0]
ntest = X_test.shape[0]
Xmat = pd.concat((X_train,X_test))
y_train = np.log(y_train)
sns.distplot(y_train);

これで正規分布になりました。もちろん最後には戻す必要があります。

欠損値について

次にデータの中に欠損値があるかどうか調べます。

Xmat_na = (Xmat.isnull().sum() / len(Xmat)) * 100
Xmat_na = Xmat_na.drop(Xmat_na[Xmat_na == 0].index).sort_values(ascending=False)[:30]
Xmat_na_missing_data = pd.DataFrame({'Missing Ratio' :Xmat_na})
Xmat_na_missing_data

	Missing Ratio
PoolQC	99.657417
MiscFeature	96.402878
Alley	93.216855
Fence	80.438506
FireplaceQu	48.646797
LotFrontage	16.649538
GarageFinish	5.447071
GarageYrBlt	5.447071
GarageQual	5.447071
GarageCond	5.447071
GarageType	5.378554
BsmtExposure	2.809181
BsmtCond	2.809181
BsmtQual	2.774923
BsmtFinType2	2.740665
BsmtFinType1	2.706406
MasVnrType	0.822199
MasVnrArea	0.787941
MSZoning	0.137033
BsmtFullBath	0.068517
BsmtHalfBath	0.068517
Utilities	0.068517
Functional	0.068517
Exterior2nd	0.034258
Exterior1st	0.034258
SaleType	0.034258
BsmtFinSF1	0.034258
BsmtFinSF2	0.034258
BsmtUnfSF	0.034258
Electrical	0.034258

ここから欠損値がある項目から数値0、文字の場合はNONEを入れます。一部のみ処理を記載します。

for col in ('BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF','TotalBsmtSF', 'BsmtFullBath', 'BsmtHalfBath'):
    Xmat[col] = Xmat[col].fillna(0)
for col in ('BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2'):
    Xmat[col] = Xmat[col].fillna('None')

また、Utilitiesの項目はほぼAllPubであり特徴量の解析に寄与しないため削除します。

Xmat = Xmat.drop(['Utilities'], axis=1)

そしてラベルコンバーターよりすべての値を実数にします。

from sklearn.preprocessing import LabelEncoder
cols = ('FireplaceQu', 'BsmtQual', 'BsmtCond', 'GarageQual', 'GarageCond', 
        'ExterQual', 'ExterCond','HeatingQC', 'PoolQC', 'KitchenQual', 'BsmtFinType1', 
        'BsmtFinType2', 'Functional', 'Fence', 'BsmtExposure', 'GarageFinish', 'LandSlope',
        'LotShape', 'PavedDrive', 'Street', 'Alley', 'CentralAir', 'MSSubClass', 'OverallCond', 
        'YrSold', 'MoSold', 'MSZoning', 'LandContour', 'LotConfig', 'Neighborhood', 'Condition1',
       'Condition2', 'BldgType', 'HouseStyle', 'RoofStyle', 'RoofMatl',
       'Exterior1st', 'Exterior2nd', 'MasVnrType', 'Foundation', 'Heating',
       'Electrical', 'GarageType', 'MiscFeature', 'SaleType', 'SaleCondition')

for c in cols:
    lbl = LabelEncoder() 
    lbl.fit(list(Xmat[c].values)) 
    Xmat[c] = lbl.transform(list(Xmat[c].values))

データはトレーニングデータ、テストデータ含めて加工しましたが、ここでトレーニングデータ、テストデータそれぞれに戻します。

X_train_1 = Xmat[:ntrain]
X_test_1 = Xmat[ntrain:]

今期のサンフレッチェ広島の観客動員数を機械学習で予測した その２（XGBoostのチューニング含む）

今週末にもJリーグが開幕するので、今週発表された週間天気を基にXGBoostのアルゴリズムを使って再度今期の観客を予測しました。今回はKaggleのHouse PricesのコンペをやっていくうちにXGBoostの使い方を覚えたので、XGBoostのチューニングの仕方も併せて乗せてみます。

今週末の天気

今週末は晴れで気温が11度の予想です。

XGBoostとは

Kaggleのようなデータ分析コンペで実績を残しているアルゴリズムです。 具体的にどのような動きをするかというのはこういうところが参考になります。 今でこそ目的関数は分かりますがMachine Learningを取る前だとそもそも何かわからなかったかも。

パッケージユーザーのための機械学習(12)：Xgboost (eXtreme Gradient Boosting) - 六本木で働くデータサイエンティストのブログ XGBoostやパラメータチューニングの仕方に関する調査 | かものはしの分析ブログ

XGBoostのチューニング

チューニングの元ネタはこちらを参考にしました。

Analytics_Vidhya/XGBoost models.ipynb at master · aarshayj/Analytics_Vidhya · GitHub

実際にチューニングした結果はこちらとなります。 まずはライブラリの読み込みと、テーブルの読み込みです。 テーブルはこちらで作成したものを使用します。

今期のサンフレッチェ広島の観客動員数を機械学習で予測した - Masser’s Blog

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from scipy.stats import norm
from sklearn.preprocessing import StandardScaler
from scipy import stats
%matplotlib inline

train = pd.read_csv('sanftrain.csv')
test = pd.read_csv('sanftest.csv')
train_ID = train['Key']
test_ID = test['Key']
y_train = train['Audience']
X_train = train.drop(['Key','Audience'], axis=1)
X_test = test.drop('Key', axis=1)

まずは、max_depthと、n_estimatorsをチューニングします。

import xgboost as xgb
from sklearn.model_selection import GridSearchCV

xgb_model = xgb.XGBRegressor(max_depth=5,min_child_weight=1,gamma=0,subsample=0.2,colsample_bytree=0.2,learning_rate=0.1)
reg = GridSearchCV(xgb_model,
                   {'max_depth': [2,3,4,5,6],
                    'n_estimators': [50,75,100,125,150,175,200,225,250,275,300]}, verbose=1)
reg.fit(X_train, y_train)

結果はこちらで分かります。これからmax_depthは4, n_estimatorsは300だと分かります。

reg.grid_scores_,reg.best_params_, reg.best_score_

   ([mean: 0.30059, std: 0.01853, params: {'max_depth': 2, 'n_estimators': 50},
      mean: 0.39571, std: 0.00777, params: {'max_depth': 2, 'n_estimators': 75},
      mean: 0.44242, std: 0.01296, params: {'max_depth': 2, 'n_estimators': 100},
      mean: 0.47352, std: 0.02664, params: {'max_depth': 2, 'n_estimators': 125},
      mean: 0.49547, std: 0.01577, params: {'max_depth': 2, 'n_estimators': 150},
      mean: 0.54260, std: 0.02425, params: {'max_depth': 2, 'n_estimators': 175},
      mean: 0.54957, std: 0.03926, params: {'max_depth': 2, 'n_estimators': 200},
      mean: 0.52969, std: 0.04004, params: {'max_depth': 2, 'n_estimators': 225},
      mean: 0.54029, std: 0.05254, params: {'max_depth': 2, 'n_estimators': 250},
      mean: 0.54485, std: 0.04447, params: {'max_depth': 2, 'n_estimators': 275},
      mean: 0.55701, std: 0.04541, params: {'max_depth': 2, 'n_estimators': 300},
      mean: 0.29411, std: 0.02525, params: {'max_depth': 3, 'n_estimators': 50},
      mean: 0.39210, std: 0.00699, params: {'max_depth': 3, 'n_estimators': 75},
      mean: 0.44110, std: 0.01439, params: {'max_depth': 3, 'n_estimators': 100},
      mean: 0.47381, std: 0.02580, params: {'max_depth': 3, 'n_estimators': 125},
      mean: 0.49691, std: 0.01668, params: {'max_depth': 3, 'n_estimators': 150},
      mean: 0.54583, std: 0.02631, params: {'max_depth': 3, 'n_estimators': 175},
      mean: 0.55439, std: 0.04097, params: {'max_depth': 3, 'n_estimators': 200},
      mean: 0.53532, std: 0.04101, params: {'max_depth': 3, 'n_estimators': 225},
      mean: 0.54449, std: 0.05164, params: {'max_depth': 3, 'n_estimators': 250},
      mean: 0.54941, std: 0.04269, params: {'max_depth': 3, 'n_estimators': 275},
      mean: 0.56011, std: 0.04418, params: {'max_depth': 3, 'n_estimators': 300},
      mean: 0.29457, std: 0.02471, params: {'max_depth': 4, 'n_estimators': 50},
      mean: 0.39260, std: 0.00701, params: {'max_depth': 4, 'n_estimators': 75},
      mean: 0.44148, std: 0.01402, params: {'max_depth': 4, 'n_estimators': 100},
      mean: 0.47413, std: 0.02567, params: {'max_depth': 4, 'n_estimators': 125},
      mean: 0.49719, std: 0.01653, params: {'max_depth': 4, 'n_estimators': 150},
      mean: 0.54507, std: 0.02510, params: {'max_depth': 4, 'n_estimators': 175},
      mean: 0.55360, std: 0.03993, params: {'max_depth': 4, 'n_estimators': 200},
      mean: 0.53499, std: 0.04068, params: {'max_depth': 4, 'n_estimators': 225},
      mean: 0.54440, std: 0.05068, params: {'max_depth': 4, 'n_estimators': 250},
      mean: 0.54975, std: 0.04196, params: {'max_depth': 4, 'n_estimators': 275},
      mean: 0.56034, std: 0.04323, params: {'max_depth': 4, 'n_estimators': 300},
      mean: 0.29457, std: 0.02471, params: {'max_depth': 5, 'n_estimators': 50},
      mean: 0.39260, std: 0.00701, params: {'max_depth': 5, 'n_estimators': 75},
      mean: 0.44148, std: 0.01402, params: {'max_depth': 5, 'n_estimators': 100},
      mean: 0.47413, std: 0.02567, params: {'max_depth': 5, 'n_estimators': 125},
      mean: 0.49719, std: 0.01653, params: {'max_depth': 5, 'n_estimators': 150},
      mean: 0.54507, std: 0.02510, params: {'max_depth': 5, 'n_estimators': 175},
      mean: 0.55360, std: 0.03993, params: {'max_depth': 5, 'n_estimators': 200},
      mean: 0.53499, std: 0.04068, params: {'max_depth': 5, 'n_estimators': 225},
      mean: 0.54440, std: 0.05068, params: {'max_depth': 5, 'n_estimators': 250},
      mean: 0.54975, std: 0.04196, params: {'max_depth': 5, 'n_estimators': 275},
      mean: 0.56034, std: 0.04323, params: {'max_depth': 5, 'n_estimators': 300},
      mean: 0.29457, std: 0.02471, params: {'max_depth': 6, 'n_estimators': 50},
      mean: 0.39260, std: 0.00701, params: {'max_depth': 6, 'n_estimators': 75},
      mean: 0.44148, std: 0.01402, params: {'max_depth': 6, 'n_estimators': 100},
      mean: 0.47413, std: 0.02567, params: {'max_depth': 6, 'n_estimators': 125},
      mean: 0.49719, std: 0.01653, params: {'max_depth': 6, 'n_estimators': 150},
      mean: 0.54507, std: 0.02510, params: {'max_depth': 6, 'n_estimators': 175},
      mean: 0.55360, std: 0.03993, params: {'max_depth': 6, 'n_estimators': 200},
      mean: 0.53499, std: 0.04068, params: {'max_depth': 6, 'n_estimators': 225},
      mean: 0.54440, std: 0.05068, params: {'max_depth': 6, 'n_estimators': 250},
      mean: 0.54975, std: 0.04196, params: {'max_depth': 6, 'n_estimators': 275},
      mean: 0.56034, std: 0.04323, params: {'max_depth': 6, 'n_estimators': 300}],
     {'max_depth': 4, 'n_estimators': 300},
     0.5603390632246239)

次に、gammaをチューニングします。

import xgboost as xgb
from sklearn.model_selection import GridSearchCV

xgb_model = xgb.XGBRegressor(max_depth=4,n_estimators=300,min_child_weight=1,gamma=0,subsample=0.2,colsample_bytree=0.2,learning_rate=0.1)
reg = GridSearchCV(xgb_model,
                   {'gamma':[i/10.0 for i in range(0,5)]}, verbose=1)
reg.fit(X_train, y_train)

結果はこちらで分かります。これからgammaは0.0だと分かります。

reg.grid_scores_,reg.best_params_, reg.best_score_

    ([mean: 0.56034, std: 0.04323, params: {'gamma': 0.0},
      mean: 0.56034, std: 0.04323, params: {'gamma': 0.1},
      mean: 0.56034, std: 0.04323, params: {'gamma': 0.2},
      mean: 0.56034, std: 0.04323, params: {'gamma': 0.3},
      mean: 0.56034, std: 0.04323, params: {'gamma': 0.4}],
     {'gamma': 0.0},
     0.5603390632246239)

次に、subsampleとcolsample_bytreeをチューニングします。

import xgboost as xgb
from sklearn.model_selection import GridSearchCV

xgb_model = xgb.XGBRegressor(max_depth=2,n_estimators=100,min_child_weight=1,gamma=0,subsample=0.2,colsample_bytree=0.2,learning_rate=0.1)
reg = GridSearchCV(xgb_model,
                    {'subsample':[i/10.0 for i in range(4,10)],
    'colsample_bytree':[i/10.0 for i in range(4,10)]}, verbose=1)
reg.fit(X_train, y_train)

結果はこちらで分かります。これからsubsampleとcolsample_bytreeはどちらも0.7だと分かります。

reg.grid_scores_,reg.best_params_, reg.best_score_

    ([mean: 0.55839, std: 0.05922, params: {'colsample_bytree': 0.4, 'subsample': 0.4},
      mean: 0.55703, std: 0.08538, params: {'colsample_bytree': 0.4, 'subsample': 0.5},
      mean: 0.55295, std: 0.09469, params: {'colsample_bytree': 0.4, 'subsample': 0.6},
      mean: 0.54332, std: 0.07073, params: {'colsample_bytree': 0.4, 'subsample': 0.7},
      mean: 0.54041, std: 0.06807, params: {'colsample_bytree': 0.4, 'subsample': 0.8},
      mean: 0.55073, std: 0.06947, params: {'colsample_bytree': 0.4, 'subsample': 0.9},
      mean: 0.56410, std: 0.06286, params: {'colsample_bytree': 0.5, 'subsample': 0.4},
      mean: 0.54650, std: 0.08238, params: {'colsample_bytree': 0.5, 'subsample': 0.5},
      mean: 0.55583, std: 0.09343, params: {'colsample_bytree': 0.5, 'subsample': 0.6},
      mean: 0.55248, std: 0.09086, params: {'colsample_bytree': 0.5, 'subsample': 0.7},
      mean: 0.54473, std: 0.08220, params: {'colsample_bytree': 0.5, 'subsample': 0.8},
      mean: 0.55531, std: 0.08295, params: {'colsample_bytree': 0.5, 'subsample': 0.9},
      mean: 0.56410, std: 0.06286, params: {'colsample_bytree': 0.6, 'subsample': 0.4},
      mean: 0.54650, std: 0.08238, params: {'colsample_bytree': 0.6, 'subsample': 0.5},
      mean: 0.55583, std: 0.09343, params: {'colsample_bytree': 0.6, 'subsample': 0.6},
      mean: 0.55248, std: 0.09086, params: {'colsample_bytree': 0.6, 'subsample': 0.7},
      mean: 0.54473, std: 0.08220, params: {'colsample_bytree': 0.6, 'subsample': 0.8},
      mean: 0.55531, std: 0.08295, params: {'colsample_bytree': 0.6, 'subsample': 0.9},
      mean: 0.55025, std: 0.07755, params: {'colsample_bytree': 0.7, 'subsample': 0.4},
      mean: 0.54751, std: 0.09015, params: {'colsample_bytree': 0.7, 'subsample': 0.5},
      mean: 0.55105, std: 0.09073, params: {'colsample_bytree': 0.7, 'subsample': 0.6},
      mean: 0.56454, std: 0.09483, params: {'colsample_bytree': 0.7, 'subsample': 0.7},
      mean: 0.54205, std: 0.10489, params: {'colsample_bytree': 0.7, 'subsample': 0.8},
      mean: 0.53911, std: 0.09839, params: {'colsample_bytree': 0.7, 'subsample': 0.9},
      mean: 0.55933, std: 0.07821, params: {'colsample_bytree': 0.8, 'subsample': 0.4},
      mean: 0.54618, std: 0.11937, params: {'colsample_bytree': 0.8, 'subsample': 0.5},
      mean: 0.55731, std: 0.11517, params: {'colsample_bytree': 0.8, 'subsample': 0.6},
      mean: 0.55182, std: 0.11408, params: {'colsample_bytree': 0.8, 'subsample': 0.7},
      mean: 0.55087, std: 0.09362, params: {'colsample_bytree': 0.8, 'subsample': 0.8},
      mean: 0.53608, std: 0.10538, params: {'colsample_bytree': 0.8, 'subsample': 0.9},
      mean: 0.55311, std: 0.09649, params: {'colsample_bytree': 0.9, 'subsample': 0.4},
      mean: 0.53675, std: 0.12155, params: {'colsample_bytree': 0.9, 'subsample': 0.5},
      mean: 0.56001, std: 0.10962, params: {'colsample_bytree': 0.9, 'subsample': 0.6},
      mean: 0.55277, std: 0.10686, params: {'colsample_bytree': 0.9, 'subsample': 0.7},
      mean: 0.55744, std: 0.09592, params: {'colsample_bytree': 0.9, 'subsample': 0.8},
      mean: 0.52927, std: 0.09682, params: {'colsample_bytree': 0.9, 'subsample': 0.9}],
     {'colsample_bytree': 0.7, 'subsample': 0.7},
     0.56453621270808296)

最後に、reg_alphaをチューニングします。

import xgboost as xgb
from sklearn.model_selection import GridSearchCV

xgb_model = xgb.XGBRegressor(max_depth=2,n_estimators=100,min_child_weight=1,gamma=0,subsample=0.7,colsample_bytree=0.7,learning_rate=0.1)
reg = GridSearchCV(xgb_model,
                    {'reg_alpha':[1e-5, 1e-2, 0.1, 1, 100]}, verbose=1)
reg.fit(X_train, y_train)

結果はreg_alphaは1となりました。

reg.grid_scores_,reg.best_params_, reg.best_score_

    ([mean: 0.56454, std: 0.09483, params: {'reg_alpha': 1e-05},
      mean: 0.56454, std: 0.09483, params: {'reg_alpha': 0.01},
      mean: 0.56454, std: 0.09483, params: {'reg_alpha': 0.1},
      mean: 0.56522, std: 0.09479, params: {'reg_alpha': 1},
      mean: 0.56443, std: 0.10188, params: {'reg_alpha': 100}],
     {'reg_alpha': 1},
     0.56521892471311341)

これらのチューニングを元にして予測します。

reg = xgb.XGBRegressor(max_depth=4,n_estimators=300,min_child_weight=1,gamma=0,subsample=0.7,colsample_bytree=0.7,learning_rate=0.1,reg_alpha=1)
reg.fit(X_train, y_train)

結果は17645人となりました。ちょっと多いかな？実際の観客がどれくらいになるか楽しみです。

Coursera Machine Learning を学習中

ここ最近機械学習を学ぶ上で評判がいいCoursera Machine Learningを学習しています。

Machine Learning | Coursera

スタンフォード大学のAndrew Ngが講義をしていますが、機械学習の権威だとか。けどこういうのが無料で学習できるのがすごくいいですね。すごく便利な世の中になりました。あと有料ですが、きちんと認定書を得ることができるとのことです。現在は通勤時間を利用して学習しています。講義の中でテストやプログラミング演習もあるとのことです。

前提知識

高校理系の数学の知識が必要となります。また、有志の方が講義の動画について日本語の字幕を追加してくださってますが、課題は英語で出されるので、英語がわからないとかなり厳しいです。

感想

Week1まで完了しました。まだプログラミングの課題はなくテストのみですが、英語で正誤問題を回答するのは単語の意味を理解しないといけないので結構厳しいですね。あとこのサイトが非常に便利です。

数学の苦手なバイオの学生がCourseraの機械学習コースを修了して気づいたこと - Qiita

廃線予定の三江線に乗ってみた

年末に帰省のついでに2018年3月31日に廃線予定の三江線に乗ってみました。

まず広島から芸備線で三次に行きました。芸備線に乗るのは高校以来です。今はバスのほうが便利だしねえ。三次で昼食を取った後、14：11発の汽車で江津方面に向かいました。やはりほとんど乗客は乗り鉄でした。中には車で並走してる人もいました。地元の人らしき人は１，２名でした。

自然豊かな景色といかにもローカル駅というしょぼい駅舎がなんとも情緒を掻き立てられます。秋はすごく混んでたとのことでしたが、紅葉時期に行けばきれいなんだろうなあと思いました。

一旦口羽駅で降り、次の電車が3時間待ち！だったので、約7km歩いて天空の城で有名な宇都井駅に向かいました。宇都井駅は結構来てる人も多く、取材も来てました。

そして18：09発の電車に乗って川戸駅に向かいました。乗り鉄しかいなかったものの、座れないほどではなかったです。

当日は川戸駅前の美川旅館に泊まりました。21時到着だったのですが、晩飯は電車の到着時間に合わせてくれて、日本酒もおいしく評判のいい宿ってのが納得でした。

そして、翌日は川戸駅7：42発の汽車に乗ってに江津に到着してほぼほぼ三江線制覇しました。その後はまだまで向かい、バスで広島に帰りました。行きは電車で9時間かかりましたが、バスは2時間ちょっとで帰りました。

けど乗ってみて分かったのですが廃線は納得です。逆に今までよく廃線しなかったなと思いました。汽車のスピードが25km/hしかでなく、時間がかかりすぎるので地域連絡線になりえず、しかも地元の人も使う人があまり見かけませんでした。せいぜい川戸の人が江津行くのに使うぐらいで、例えば石見川本の人は江津ではなく車で大田に向かうんですよね。しかも外部の人が利用するにも紅葉時期ぐらいしか需要を見込めずろくな観光地がありません。それに加えて地元の人に路線を残そうというやる気が感じられませんでした。

戦前は木材運ぶのに需要があったみたいですが今はトラックでもできるのでとっくに役目を終えてる路線ですね。しかも地元の方が反対してたのは乗る人が困る理由ではなく、JR西日本が払う固定資産税のみが目的だったという話もあります。この路線の建設にあたっても地元で相当もめていて、建設が終わった1975年はすでにモータリゼーションの世の中になってました。

ameblo.jp

地域のインフラを生かすも殺すもその地元次第ということが痛切によくわかりました。地方の悲惨な状況を考えるいい機会となりますので是非乗っていただければと思います。

今期のサンフレッチェ広島の観客動員数を機械学習で予測した

私は広島出身ということもあり、サンフレッチェ広島のサポーターでありますが、これを機械学習で生かせないか考えてみました。そこで2009年から2017年までの観客データやその日の天候、その他の要素から2018年の観客動員数を機械学習で予測してみました。ただ、Python学習してまだ2日であり、あくまで機械学習の勉強のために作成したので、結果としてはしょぼい結果となってしまいますがご了承のほどよろしくお願いします。

観客動員数の分布

2009年から2017年までの観客動員数の分布をグラフで表すとこんな感じとなりました。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
df = pd.read_csv("sanfdata1.csv", encoding="SHIFT_JIS")
y_Audi =df['Audience']
ax = sns.distplot(y_Audi)
plt.show()

上のソースコードで下記のグラフが表示できます。

予測のための要素

観客動員に影響する要素として以下の項目を考えました。

* 節数

ホーム開幕戦（1節、2節）、ホーム最終戦（33節、34節）、その他に分けて平均を算出しました。算出したソースコードは次のようになります。

from pandas import DataFrame
grouped = df.groupby('Setu')
grouped.mean().Audience

結果の平均はこのようになりました。ホーム開幕戦、ホーム最終戦は多い傾向があります。

対戦相手	動員数
ホーム開幕戦	19230.909091
ホーム最終戦	23149.000000
その他	14583.143939

* 対戦相手

対戦相手ごとの観客動員の平均を算出しました。算出したソースコードは次のようになります。

grouped = df.groupby('Team')
grouped.mean().Audience

結果の平均はこのようになりました。やはり浦和が一番多く、第二グループとして鹿島、F東京、横浜M、大阪の2チームでしょうか。面白いのは湘南が平均2位ですが、これは2013年、2015年と優勝が懸かったホーム最終戦と被ったためです。あと京都、福岡もホーム最終戦の影響で多いですね。

対戦相手	動員数
鹿島	15760
浦和	23658.777778
大宮	14724
千葉	11572
鹿島	13872.75
F東京	17636.25
川崎	13802
横浜M	15502.222222
湘南	21319.75
甲府	12737.833333
松本	7956
新潟	13785.444444
清水	14709.75
磐田	14162.857143
名古屋	14064.714286
京都	17748
G大阪	16970
C大阪	16996.428571
神戸	12691.625
徳島	11776
福岡	16413.5
鳥栖	15164.166667
大分	10829.5

* 順位

順位ごとの観客動員の平均を算出しました。算出したソースコードは次のようになります。

grouped = df.groupby('Place')
grouped.mean().Audience

結果は次のようになりました。1位が飛びぬけて多く、勝たないと客が来ない広島県民の気質が表れてますね。カープの観客動員数で集計するともっと顕著かもしれません。

順位	動員数
1位	21162.000000
2位ー9位	15093.833333
10位以下	12980.318182

* 天気

順位ごとの観客動員の平均を算出しました。算出したソースコードは次のようになります。

grouped = df.groupby('Weather')
grouped.mean().Audience

結果は次のようになりました。予想通り雨が降ると観客動員が減りますね。

天気	動員数
晴れ、曇り	15839.6875
雨	13224.0000

* 気温

気温ごとの平均を算出しました。算出したソースコードは次のようになります。

grouped = df.groupby('Heat')
grouped.mean().Audience

結果は次のようになりました。寒いときに観客動員が多いのはホーム開幕戦、ホーム最終戦の影響があるかもしれません。

気温	動員数
10度以下	17890.923077
11度ー15度	17294.454545
16度ー20度	14713.500000
21度ー25度	14258.352941
26度以上	15699.647059

* 時間帯

時間帯ごとの観客動員の平均を算出しました。算出したソースコードは次のようになります。

grouped = df.groupby('Time')
grouped.mean().Audience

結果は次のようになりました。夜開催が少ないのは平日開催の影響が多いでしょうね。

天気	動員数
1200-1600	16375.242857
1700以降	14616.939024

* 平日、休日

平日開催、休日開催ごとの観客動員の平均を算出しました。算出したソースコードは次のようになります。

grouped = df.groupby('Time')
grouped.mean().Audience

結果は次のようになりました。明らかに平日開催のほうが少ないですね。

平日、休日	動員数
平日	10813.210526
休日	16085.751880

* カープ開催あり、なし

同日にカープの試合があるかどうかで観客動員の平均を算出しました。ただ同日でも時間帯が違う場合は別開催としています。算出したソースコードは次のようになります。

grouped = df.groupby('Carp')
grouped.mean().Audience

結果は次のようになりました。カープ開催なしのほうが観客動員多いです。やはりこの辺はなんとかならないですかね。

カープ開催あり、なし	動員数
なし	15868.693069
あり	14551.333333

ランダムフォレストを用いた機械学習による予測

ランダムフォレストを使いまして第1節札幌戦の予測をしてみました。（本当はxgboostとか試したかったがあまりいい結果がでなかった）

import pandas as pd
from sklearn.ensemble import RandomForestClassifier
test_df = pd.read_csv("sanftest.csv", encoding="SHIFT_JIS")
train_df = pd.read_csv("sanftrain.csv", encoding="SHIFT_JIS")
train_X = train_df[['Setu','Team','Place','Weather','Heat','Time','Holiday','Carp']]
test_X = test_df[['Setu','Team','Place','Weather','Heat','Time','Holiday','Carp']]
train_Y = train_df['Audience']
clf = RandomForestClassifier(random_state=0)
clf = clf.fit(train_X, train_Y)
clf.predict(test_X)

結果はこうなりました。天気予報が出たらまた改めて予測してみたいですね。また、もっと良さそうなアルゴリズムがあったら試したいです。

対戦相手	予想動員数
第1節	14671

考察

今回は天候、時間帯、対戦相手等表に出ているデータのみで算出しましたが、表に出ていない年パスの数、歩いて行ける範囲の駐車場数、直前のチケットの売り上げ枚数の要素を加えるともっと正確な数が予測できるのではないかと思います。おそらくクラブはやっているはずで、これを基にしてシャトルバスの本数とか予測してるでしょう。あともう少しデータがほしいですね。