In [2]:
# from __future__ import exam_success
from __future__ import absolute_import
from __future__ import print_function
%matplotlib inline
import sklearn
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import random
import pandas as pd
import scipy.stats as stats
# Sk cheats
from sklearn.cross_validation import cross_val_score # cross val
from sklearn.ensemble import RandomForestRegressor
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.preprocessing import Imputer # get rid of nan
Predictions is made in the USA corn growing states (mainly Iowa, Illinois, Indiana) during the season with the highest rainfall (as illustrated by Iowa for the april to august months)
The Kaggle page indicate that the dataset have been shuffled, so working on a subset seems acceptable
The test set is not a extracted from the same data as the training set however, which make the evaluation trickier
In [3]:
%%time
filename = "data/reduced_train_100000.csv"
#filename = "data/reduced_train_100000.csv"
raw = pd.read_csv(filename)
raw = raw.set_index('Id')
In [4]:
raw['Expected'].describe()
Out[4]:
Per wikipedia, a value of more than 421 mm/h is considered "Extreme/large hail"
If we encounter the value 327.40 meter per hour, we should probably start building Noah's ark
Therefor, it seems reasonable to drop values too large, considered as outliers
In [5]:
# Considering that the gauge may concentrate the rainfall, we set the cap to 1000
# Comment this line to analyse the complete dataset
l = len(raw)
raw = raw[raw['Expected'] < 1000]
print("Dropped %d (%0.2f%%)"%(l-len(raw),(l-len(raw))/float(l)*100))
In [6]:
raw.head(5)
Out[6]:
In [7]:
raw.describe()
Out[7]:
We regroup the data by ID
In [8]:
# We select all features except for the minutes past,
# because we ignore the time repartition of the sequence for now
features_columns = list([u'Ref', u'Ref_5x5_10th',
u'Ref_5x5_50th', u'Ref_5x5_90th', u'RefComposite',
u'RefComposite_5x5_10th', u'RefComposite_5x5_50th',
u'RefComposite_5x5_90th', u'RhoHV', u'RhoHV_5x5_10th',
u'RhoHV_5x5_50th', u'RhoHV_5x5_90th', u'Zdr', u'Zdr_5x5_10th',
u'Zdr_5x5_50th', u'Zdr_5x5_90th', u'Kdp', u'Kdp_5x5_10th',
u'Kdp_5x5_50th', u'Kdp_5x5_90th'])
def getXy(raw):
selected_columns = list([ u'minutes_past',u'radardist_km', u'Ref', u'Ref_5x5_10th',
u'Ref_5x5_50th', u'Ref_5x5_90th', u'RefComposite',
u'RefComposite_5x5_10th', u'RefComposite_5x5_50th',
u'RefComposite_5x5_90th', u'RhoHV', u'RhoHV_5x5_10th',
u'RhoHV_5x5_50th', u'RhoHV_5x5_90th', u'Zdr', u'Zdr_5x5_10th',
u'Zdr_5x5_50th', u'Zdr_5x5_90th', u'Kdp', u'Kdp_5x5_10th',
u'Kdp_5x5_50th', u'Kdp_5x5_90th'])
data = raw[selected_columns]
docX, docY = [], []
for i in data.index.unique():
if isinstance(data.loc[i],pd.core.series.Series):
m = [data.loc[i].as_matrix()]
docX.append(m)
docY.append(float(raw.loc[i]["Expected"]))
else:
m = data.loc[i].as_matrix()
docX.append(m)
docY.append(float(raw.loc[i][:1]["Expected"]))
X , y = np.array(docX) , np.array(docY)
return X,y
In [9]:
#noAnyNan = raw.loc[raw[features_columns].dropna(how='any').index.unique()]
noAnyNan = raw.dropna()
In [10]:
noFullNan = raw.loc[raw[features_columns].dropna(how='all').index.unique()]
In [11]:
fullNan = raw.drop(raw[features_columns].dropna(how='all').index)
As a first try, we make predictions on the complete data, and return the 50th percentile and uncomplete and fully empty data
In [12]:
%%time
X,y=getXy(noAnyNan)
In [13]:
%%time
#XX = [np.array(t).mean(0) for t in X]
XX = [np.append(np.nanmean(np.array(t),0),(np.array(t)[1:] - np.array(t)[:-1]).sum(0) ) for t in X]
In [17]:
t = X[0]
In [30]:
t
Out[30]:
In [31]:
np.percentile(t,50,axis=0)
Out[31]:
In [32]:
XX=[]
for t in X:
#print(idx)
tmp = np.append(np.nanmean(np.array(t),0),(np.array(t)[1:] - np.array(t)[:-1]).sum(0) )
tmp = np.append(tmp,np.percentile(t,10,axis=0))
tmp = np.append(tmp,np.percentile(t,50,axis=0))
tmp = np.append(tmp,np.percentile(t,90,axis=0))
test = t
try:
taa=test[:,0]
except TypeError:
taa=[test[0][0]]
valid_time = np.zeros_like(taa)
valid_time[0] = taa[0]
for n in xrange(1,len(taa)):
valid_time[n] = taa[n] - taa[n-1]
valid_time[-1] = valid_time[-1] + 60 - np.sum(valid_time)
valid_time = valid_time / 60.0
sum=0
try:
column_ref=test[:,2]
except TypeError:
column_ref=[test[0][2]]
for dbz, hours in zip(column_ref, valid_time):
# See: https://en.wikipedia.org/wiki/DBZ_(meteorology)
if np.isfinite(dbz):
mmperhr = pow(pow(10, dbz/10)/200, 0.625)
sum = sum + mmperhr * hours
#XX.append(np.append(np.array(sum),tmp))
#XX.append(np.array([sum]))
XX.append(tmp)
In [16]:
XX[0]
Out[16]:
In [263]:
XX[2]
Out[263]:
In [33]:
def splitTrainTest(X, y, split=0.2):
tmp1, tmp2 = [], []
ps = int(len(X) * (1-split))
index_shuf = range(len(X))
random.shuffle(index_shuf)
for i in index_shuf:
tmp1.append(X[i])
tmp2.append(y[i])
return tmp1[:ps], tmp2[:ps], tmp1[ps:], tmp2[ps:]
In [34]:
X_train,y_train, X_test, y_test = splitTrainTest(XX,y)
In [225]:
from sklearn import svm
In [283]:
svr = svm.SVR(C=100000)
In [284]:
%%time
srv = svr.fit(X_train,y_train)
In [240]:
err = (svr.predict(X_test)-y_test)**2
err.sum()/len(err)
Out[240]:
In [285]:
err = (svr.predict(X_test)-y_test)**2
err.sum()/len(err)
Out[285]:
In [241]:
%%time
svr_score = cross_val_score(svr, XX, y, cv=5)
print("Score: %s\tMean: %.03f"%(svr_score,svr_score.mean()))
In [367]:
from sklearn.neighbors import KNeighborsRegressor
from sklearn import grid_search
In [378]:
knn = KNeighborsRegressor(n_neighbors=6,weights='distance')
In [ ]:
parameters = {'weights':('distance','uniform'),'algorithm':('auto', 'ball_tree', 'kd_tree', 'brute')}
grid_knn = grid_search.GridSearchCV(knn, parameters,n_jobs=2)
In [ ]:
%%time
grid_knn.fit(X_train,y_train)
In [ ]:
print(grid_knn.grid_scores_)
print("Best: ",grid_knn.best_params_)
In [375]:
knn = grid_knn.best_estimator_
In [379]:
knn.fit(X_train,y_train)
Out[379]:
In [380]:
err = (knn.predict(X_test)-y_test)**2
err.sum()/len(err)
Out[380]:
In [35]:
etreg = ExtraTreesRegressor(n_estimators=200, max_depth=None, min_samples_split=1, random_state=0)
In [69]:
parameters = {'n_estimators':range(100,200,20)}
grid_rf = grid_search.GridSearchCV(etreg, parameters,n_jobs=2)
In [70]:
%%time
grid_rf.fit(X_train,y_train)
Out[70]:
In [91]:
print(grid_rf.grid_scores_)
print("Best: ",grid_rf.best_params_)
In [86]:
grid_rf.best_params_
Out[86]:
In [289]:
es = etreg
#es = grid_rf.best_estimator_
In [290]:
%%time
es = es.fit(X_train,y_train)
In [291]:
err = (es.predict(X_test)-y_test)**2
err.sum()/len(err)
Out[291]:
In [296]:
from keras.models import Sequential
from keras.layers.core import Dense, Dropout, Activation
from keras.optimizers import SGD
in_dim = len(XX[0])
out_dim = 1
model = Sequential()
# Dense(64) is a fully-connected layer with 64 hidden units.
# in the first layer, you must specify the expected input data shape:
# here, 20-dimensional vectors.
model.add(Dense(128, input_shape=(in_dim,)))
model.add(Activation('tanh'))
model.add(Dropout(0.5))
model.add(Dense(1, init='uniform'))
model.add(Activation('linear'))
#sgd = SGD(lr=0.1, decay=1e-6, momentum=0.9, nesterov=True)
#model.compile(loss='mean_squared_error', optimizer=sgd)
rms = RMSprop()
model.compile(loss='mean_squared_error', optimizer=rms)
#model.fit(X_train, y_train, nb_epoch=20, batch_size=16)
#score = model.evaluate(X_test, y_test, batch_size=16)
In [309]:
prep = []
for i in y_train:
prep.append(min(i,20))
In [327]:
prep=np.array(prep)
mi,ma = prep.min(),prep.max()
fy = (prep-mi) / (ma-mi)
#my = fy.max()
#fy = fy/fy.max()
In [329]:
model.fit(np.array(X_train), fy, batch_size=10, nb_epoch=10, validation_split=0.1)
Out[329]:
In [335]:
pred = model.predict(np.array(X_test))*ma+mi
In [336]:
err = (pred-y_test)**2
err.sum()/len(err)
Out[336]:
In [357]:
r = random.randrange(len(X_train))
print("(Train) Prediction %0.4f, True: %0.4f"%(model.predict(np.array([X_train[r]]))[0][0]*ma+mi,y_train[r]))
r = random.randrange(len(X_test))
print("(Test) Prediction %0.4f, True: %0.4f"%(model.predict(np.array([X_test[r]]))[0][0]*ma+mi,y_test[r]))
In [38]:
def marshall_palmer(ref, minutes_past):
#print("Estimating rainfall from {0} observations".format(len(minutes_past)))
# how long is each observation valid?
valid_time = np.zeros_like(minutes_past)
valid_time[0] = minutes_past.iloc[0]
for n in xrange(1, len(minutes_past)):
valid_time[n] = minutes_past.iloc[n] - minutes_past.iloc[n-1]
valid_time[-1] = valid_time[-1] + 60 - np.sum(valid_time)
valid_time = valid_time / 60.0
# sum up rainrate * validtime
sum = 0
for dbz, hours in zip(ref, valid_time):
# See: https://en.wikipedia.org/wiki/DBZ_(meteorology)
if np.isfinite(dbz):
mmperhr = pow(pow(10, dbz/10)/200, 0.625)
sum = sum + mmperhr * hours
return sum
def simplesum(ref,hour):
hour.sum()
# each unique Id is an hour of data at some gauge
def myfunc(hour):
#rowid = hour['Id'].iloc[0]
# sort hour by minutes_past
hour = hour.sort('minutes_past', ascending=True)
est = marshall_palmer(hour['Ref'], hour['minutes_past'])
return est
In [140]:
info = raw.groupby(raw.index)
In [39]:
estimates = raw.groupby(raw.index).apply(myfunc)
estimates.head(20)
Out[39]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [113]:
%%time
etreg.fit(X_train,y_train)
Out[113]:
In [96]:
%%time
et_score = cross_val_score(etreg, XX, y, cv=5)
print("Score: %s\tMean: %.03f"%(et_score,et_score.mean()))
In [114]:
%%time
et_score = cross_val_score(etreg, XX, y, cv=5)
print("Score: %s\tMean: %.03f"%(et_score,et_score.mean()))
In [97]:
err = (etreg.predict(X_test)-y_test)**2
err.sum()/len(err)
Out[97]:
In [115]:
err = (etreg.predict(X_test)-y_test)**2
err.sum()/len(err)
Out[115]:
In [41]:
r = random.randrange(len(X_train))
print(r)
print(etreg.predict(X_train[r]))
print(y_train[r])
r = random.randrange(len(X_test))
print(r)
print(etreg.predict(X_test[r]))
print(y_test[r])
In [13]:
%%time
#filename = "data/reduced_test_5000.csv"
filename = "data/test.csv"
test = pd.read_csv(filename)
test = test.set_index('Id')
In [14]:
features_columns = list([u'Ref', u'Ref_5x5_10th',
u'Ref_5x5_50th', u'Ref_5x5_90th', u'RefComposite',
u'RefComposite_5x5_10th', u'RefComposite_5x5_50th',
u'RefComposite_5x5_90th', u'RhoHV', u'RhoHV_5x5_10th',
u'RhoHV_5x5_50th', u'RhoHV_5x5_90th', u'Zdr', u'Zdr_5x5_10th',
u'Zdr_5x5_50th', u'Zdr_5x5_90th', u'Kdp', u'Kdp_5x5_10th',
u'Kdp_5x5_50th', u'Kdp_5x5_90th'])
def getX(raw):
selected_columns = list([ u'radardist_km', u'Ref', u'Ref_5x5_10th',
u'Ref_5x5_50th', u'Ref_5x5_90th', u'RefComposite',
u'RefComposite_5x5_10th', u'RefComposite_5x5_50th',
u'RefComposite_5x5_90th', u'RhoHV', u'RhoHV_5x5_10th',
u'RhoHV_5x5_50th', u'RhoHV_5x5_90th', u'Zdr', u'Zdr_5x5_10th',
u'Zdr_5x5_50th', u'Zdr_5x5_90th', u'Kdp', u'Kdp_5x5_10th',
u'Kdp_5x5_50th', u'Kdp_5x5_90th'])
data = raw[selected_columns]
docX= []
for i in data.index.unique():
if isinstance(data.loc[i],pd.core.series.Series):
m = [data.loc[i].as_matrix()]
docX.append(m)
else:
m = data.loc[i].as_matrix()
docX.append(m)
X = np.array(docX)
return X
In [74]:
X=getX(test)
tmp = []
for i in X:
tmp.append(len(i))
tmp = np.array(tmp)
sns.countplot(tmp,order=range(tmp.min(),tmp.max()+1))
plt.title("Number of ID per number of observations\n(On test dataset)")
plt.plot()
Out[74]:
In [15]:
testFull = test.dropna()
In [16]:
%%time
X=getX(testFull) # 1min
XX = [np.array(t).mean(0) for t in X] # 10s
In [17]:
pd.DataFrame(etreg.predict(XX)).describe()
Out[17]:
In [18]:
predFull = zip(testFull.index.unique(),etreg.predict(XX))
In [36]:
testNan = test.drop(test[features_columns].dropna(how='all').index)
In [37]:
tmp = np.empty(len(testNan))
tmp.fill(0.445000) # 50th percentile of full Nan dataset
predNan = zip(testNan.index.unique(),tmp)
In [39]:
testLeft = test.drop(testNan.index.unique()).drop(testFull.index.unique())
In [40]:
tmp = np.empty(len(testLeft))
tmp.fill(1.27) # 50th percentile of full Nan dataset
predLeft = zip(testLeft.index.unique(),tmp)
In [41]:
len(testFull.index.unique())
Out[41]:
In [42]:
len(testNan.index.unique())
Out[42]:
In [43]:
len(testLeft.index.unique())
Out[43]:
In [44]:
pred = predFull + predNan + predLeft
In [45]:
pred.sort(key=lambda x: x[0], reverse=False)
In [46]:
submission = pd.DataFrame(pred)
submission.columns = ["Id","Expected"]
submission.head()
Out[46]:
In [47]:
submission.to_csv("first_submit.csv",index=False)
In [ ]:
In [34]:
filename = "data/sample_solution.csv"
sol = pd.read_csv(filename)
In [28]:
ss = np.array(sol)
In [30]:
%%time
for a,b in predFull:
ss[a-1][1]=b
In [41]:
ss
Out[41]:
In [45]:
sub = pd.DataFrame(ss)
sub.columns = ["Id","Expected"]
sub.Id = sub.Id.astype(int)
sub.head()
Out[45]:
In [46]:
sub.to_csv("submit_2.csv",index=False)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: