829 lines
31 KiB
Python
829 lines
31 KiB
Python
#!/usr/bin/python
|
|
# -*- coding: utf8 -*-
|
|
|
|
import numpy as np
|
|
import pandas as pd
|
|
import matplotlib as plt
|
|
import matplotlib.colors as pltcolors
|
|
import matplotlib.cm as cmx
|
|
import matplotlib.pyplot as plt
|
|
from mpl_toolkits.mplot3d import Axes3D
|
|
# from sklearn.cross_validation import KFold
|
|
from pyFTS.benchmarks import Measures, naive, ResidualAnalysis, ProbabilityDistribution
|
|
from pyFTS.partitioners import Grid
|
|
from pyFTS.common import Membership, FuzzySet, FLR, Transformations, Util
|
|
from pyFTS import fts, chen, yu, ismailefendi, sadaei, hofts, hwang, pfts, ifts
|
|
|
|
colors = ['grey', 'rosybrown', 'maroon', 'red','orange', 'yellow', 'olive', 'green',
|
|
'cyan', 'blue', 'darkblue', 'purple', 'darkviolet']
|
|
|
|
ncol = len(colors)
|
|
|
|
styles = ['-','--','-.',':','.']
|
|
|
|
nsty = len(styles)
|
|
|
|
def allPointForecasters(data_train, data_test, partitions, max_order=3, statistics=True, residuals=True,
|
|
series=True, save=False, file=None, tam=[20, 5], models=None, transformation=None,
|
|
distributions=False):
|
|
|
|
if models is None:
|
|
models = [naive.Naive, chen.ConventionalFTS, yu.WeightedFTS, ismailefendi.ImprovedWeightedFTS,
|
|
sadaei.ExponentialyWeightedFTS, hofts.HighOrderFTS, pfts.ProbabilisticFTS]
|
|
|
|
objs = []
|
|
|
|
if transformation is not None:
|
|
data_train_fs = Grid.GridPartitionerTrimf(transformation.apply(data_train),partitions)
|
|
else:
|
|
data_train_fs = Grid.GridPartitionerTrimf(data_train, partitions)
|
|
|
|
count = 1
|
|
|
|
lcolors = []
|
|
|
|
for count, model in enumerate(models, start=0):
|
|
#print(model)
|
|
mfts = model("")
|
|
if not mfts.isHighOrder:
|
|
if transformation is not None:
|
|
mfts.appendTransformation(transformation)
|
|
mfts.train(data_train, data_train_fs)
|
|
objs.append(mfts)
|
|
lcolors.append( colors[count % ncol] )
|
|
else:
|
|
for order in np.arange(1,max_order+1):
|
|
if order >= mfts.minOrder:
|
|
mfts = model(" n = " + str(order))
|
|
if transformation is not None:
|
|
mfts.appendTransformation(transformation)
|
|
mfts.train(data_train, data_train_fs, order=order)
|
|
objs.append(mfts)
|
|
lcolors.append(colors[(count + order) % ncol])
|
|
|
|
if statistics:
|
|
print(getPointStatistics(data_test, objs))
|
|
|
|
if residuals:
|
|
print(ResidualAnalysis.compareResiduals(data_test, objs))
|
|
ResidualAnalysis.plotResiduals2(data_test, objs, save=save, file=file, tam=tam)
|
|
|
|
if series:
|
|
plotComparedSeries(data_test, objs, lcolors, typeonlegend=False, save=save, file=file, tam=tam,
|
|
intervals=False)
|
|
|
|
if distributions:
|
|
lcolors.insert(0,'black')
|
|
pmfs = []
|
|
pmfs.append(
|
|
ProbabilityDistribution.ProbabilityDistribution("Original", 100, [min(data_train), max(data_train)], data=data_train) )
|
|
|
|
for m in objs:
|
|
forecasts = m.forecast(data_train)
|
|
pmfs.append(
|
|
ProbabilityDistribution.ProbabilityDistribution(m.shortname, 100, [min(data_train), max(data_train)],
|
|
data=forecasts))
|
|
print(getProbabilityDistributionStatistics(pmfs,data_train))
|
|
|
|
plotProbabilityDistributions(pmfs, lcolors)
|
|
|
|
|
|
def getPointStatistics(data, models, externalmodels = None, externalforecasts = None, indexers=None):
|
|
ret = "Model & Order & RMSE & SMAPE & Theil's U \\\\ \n"
|
|
for count,model in enumerate(models,start=0):
|
|
if indexers is not None and indexers[count] is not None:
|
|
ndata = np.array(indexers[count].get_data(data[model.order:]))
|
|
else:
|
|
ndata = np.array(data[model.order:])
|
|
|
|
if model.isMultivariate or indexers is None:
|
|
forecasts = model.forecast(data)
|
|
elif not model.isMultivariate and indexers is not None:
|
|
forecasts = model.forecast( indexers[count].get_data(data) )
|
|
|
|
if model.hasSeasonality:
|
|
nforecasts = np.array(forecasts)
|
|
else:
|
|
nforecasts = np.array(forecasts[:-1])
|
|
ret += model.shortname + " & "
|
|
ret += str(model.order) + " & "
|
|
ret += str(round(Measures.rmse(ndata, nforecasts), 2)) + " & "
|
|
ret += str(round(Measures.smape(ndata, nforecasts), 2))+ " & "
|
|
ret += str(round(Measures.UStatistic(ndata, nforecasts), 2))
|
|
#ret += str(round(Measures.TheilsInequality(np.array(data[fts.order:]), np.array(forecasts[:-1])), 4))
|
|
ret += " \\\\ \n"
|
|
if externalmodels is not None:
|
|
l = len(externalmodels)
|
|
for k in np.arange(0,l):
|
|
ret += externalmodels[k] + " & "
|
|
ret += " 1 & "
|
|
ret += str(round(Measures.rmse(data, externalforecasts[k][:-1]), 2)) + " & "
|
|
ret += str(round(Measures.smape(data, externalforecasts[k][:-1]), 2))+ " & "
|
|
ret += str(round(Measures.UStatistic(data, externalforecasts[k][:-1]), 2))
|
|
ret += " \\\\ \n"
|
|
return ret
|
|
|
|
|
|
def getProbabilityDistributionStatistics(pmfs, data):
|
|
ret = "Model & Entropy & Empirical Likelihood & Pseudo Likelihood \\\\ \n"
|
|
for k in pmfs:
|
|
ret += k.name + " & "
|
|
ret += str(k.entropy()) + " & "
|
|
ret += str(k.empiricalloglikelihood())+ " & "
|
|
ret += str(k.pseudologlikelihood(data))
|
|
ret += " \\\\ \n"
|
|
return ret
|
|
|
|
|
|
def allIntervalForecasters(data_train, data_test, partitions, max_order=3,save=False, file=None, tam=[20, 5],
|
|
models=None, transformation=None):
|
|
if models is None:
|
|
models = [ifts.IntervalFTS, pfts.ProbabilisticFTS]
|
|
|
|
objs = []
|
|
|
|
if transformation is not None:
|
|
data_train_fs = Grid.GridPartitionerTrimf(transformation.apply(data_train),partitions)
|
|
else:
|
|
data_train_fs = Grid.GridPartitionerTrimf(data_train, partitions)
|
|
|
|
lcolors = []
|
|
|
|
for count, model in Util.enumerate2(models, start=0, step=2):
|
|
mfts = model("")
|
|
if not mfts.isHighOrder:
|
|
if transformation is not None:
|
|
mfts.appendTransformation(transformation)
|
|
mfts.train(data_train, data_train_fs)
|
|
objs.append(mfts)
|
|
lcolors.append( colors[count % ncol] )
|
|
else:
|
|
for order in np.arange(1,max_order+1):
|
|
if order >= mfts.minOrder:
|
|
mfts = model(" n = " + str(order))
|
|
if transformation is not None:
|
|
mfts.appendTransformation(transformation)
|
|
mfts.train(data_train, data_train_fs, order=order)
|
|
objs.append(mfts)
|
|
lcolors.append(colors[count % ncol])
|
|
|
|
print(getIntervalStatistics(data_test, objs))
|
|
|
|
plotComparedSeries(data_test, objs, lcolors, typeonlegend=False, save=save, file=file, tam=tam, intervals=True)
|
|
|
|
|
|
def getIntervalStatistics(original, models):
|
|
ret = "Model & Order & Sharpness & Resolution & Coverage \\\\ \n"
|
|
for fts in models:
|
|
forecasts = fts.forecastInterval(original)
|
|
ret += fts.shortname + " & "
|
|
ret += str(fts.order) + " & "
|
|
ret += str(round(Measures.sharpness(forecasts), 2)) + " & "
|
|
ret += str(round(Measures.resolution(forecasts), 2)) + " & "
|
|
ret += str(round(Measures.coverage(original[fts.order:], forecasts[:-1]), 2)) + " \\\\ \n"
|
|
return ret
|
|
|
|
|
|
def plotDistribution(dist):
|
|
for k in dist.index:
|
|
alpha = np.array([dist[x][k] for x in dist]) * 100
|
|
x = [k for x in np.arange(0, len(alpha))]
|
|
y = dist.columns
|
|
plt.scatter(x, y, c=alpha, marker='s', linewidths=0, cmap='Oranges', norm=pltcolors.Normalize(vmin=0, vmax=1),
|
|
vmin=0, vmax=1, edgecolors=None)
|
|
|
|
|
|
def plotComparedSeries(original, models, colors, typeonlegend=False, save=False, file=None, tam=[20, 5],
|
|
intervals=True):
|
|
fig = plt.figure(figsize=tam)
|
|
ax = fig.add_subplot(111)
|
|
|
|
mi = []
|
|
ma = []
|
|
|
|
legends = []
|
|
|
|
ax.plot(original, color='black', label="Original", linewidth=1.5)
|
|
|
|
for count, fts in enumerate(models, start=0):
|
|
if fts.hasPointForecasting and not intervals:
|
|
forecasted = fts.forecast(original)
|
|
mi.append(min(forecasted) * 0.95)
|
|
ma.append(max(forecasted) * 1.05)
|
|
for k in np.arange(0, fts.order):
|
|
forecasted.insert(0, None)
|
|
lbl = fts.shortname
|
|
if typeonlegend: lbl += " (Point)"
|
|
ax.plot(forecasted, color=colors[count], label=lbl, ls="-")
|
|
|
|
if fts.hasIntervalForecasting and intervals:
|
|
forecasted = fts.forecastInterval(original)
|
|
lower = [kk[0] for kk in forecasted]
|
|
upper = [kk[1] for kk in forecasted]
|
|
mi.append(min(lower) * 0.95)
|
|
ma.append(max(upper) * 1.05)
|
|
for k in np.arange(0, fts.order):
|
|
lower.insert(0, None)
|
|
upper.insert(0, None)
|
|
lbl = fts.shortname
|
|
if typeonlegend: lbl += " (Interval)"
|
|
ax.plot(lower, color=colors[count], label=lbl, ls="-")
|
|
ax.plot(upper, color=colors[count], ls="-")
|
|
|
|
handles0, labels0 = ax.get_legend_handles_labels()
|
|
lgd = ax.legend(handles0, labels0, loc=2, bbox_to_anchor=(1, 1))
|
|
legends.append(lgd)
|
|
|
|
# ax.set_title(fts.name)
|
|
ax.set_ylim([min(mi), max(ma)])
|
|
ax.set_ylabel('F(T)')
|
|
ax.set_xlabel('T')
|
|
ax.set_xlim([0, len(original)])
|
|
|
|
Util.showAndSaveImage(fig, file, save, lgd=legends)
|
|
|
|
|
|
def plotProbabilityDistributions(pmfs,lcolors):
|
|
fig = plt.figure(figsize=[15, 7])
|
|
ax = fig.add_subplot(111)
|
|
|
|
for k,m in enumerate(pmfs,start=0):
|
|
m.plot(ax, color=lcolors[k])
|
|
|
|
handles0, labels0 = ax.get_legend_handles_labels()
|
|
ax.legend(handles0, labels0)
|
|
|
|
|
|
def allAheadForecasters(data_train, data_test, partitions, start, steps, resolution = None, max_order=3,save=False, file=None, tam=[20, 5],
|
|
models=None, transformation=None, option=2):
|
|
if models is None:
|
|
models = [pfts.ProbabilisticFTS]
|
|
|
|
if resolution is None: resolution = (max(data_train) - min(data_train)) / 100
|
|
|
|
objs = []
|
|
|
|
if transformation is not None:
|
|
data_train_fs = Grid.GridPartitionerTrimf(transformation.apply(data_train),partitions)
|
|
else:
|
|
data_train_fs = Grid.GridPartitionerTrimf(data_train, partitions)
|
|
|
|
lcolors = []
|
|
|
|
for count, model in Util.enumerate2(models, start=0, step=2):
|
|
mfts = model("")
|
|
if not mfts.isHighOrder:
|
|
if transformation is not None:
|
|
mfts.appendTransformation(transformation)
|
|
mfts.train(data_train, data_train_fs)
|
|
objs.append(mfts)
|
|
lcolors.append( colors[count % ncol] )
|
|
else:
|
|
for order in np.arange(1,max_order+1):
|
|
if order >= mfts.minOrder:
|
|
mfts = model(" n = " + str(order))
|
|
if transformation is not None:
|
|
mfts.appendTransformation(transformation)
|
|
mfts.train(data_train, data_train_fs, order=order)
|
|
objs.append(mfts)
|
|
lcolors.append(colors[count % ncol])
|
|
|
|
distributions = [False for k in objs]
|
|
|
|
distributions[0] = True
|
|
|
|
print(getDistributionStatistics(data_test[start:], objs, steps, resolution))
|
|
|
|
#plotComparedIntervalsAhead(data_test, objs, lcolors, distributions=, save=save, file=file, tam=tam, intervals=True)
|
|
|
|
|
|
def getDistributionStatistics(original, models, steps, resolution):
|
|
ret = "Model & Order & Interval & Distribution \\\\ \n"
|
|
for fts in models:
|
|
densities1 = fts.forecastAheadDistribution(original,steps,resolution, parameters=3)
|
|
densities2 = fts.forecastAheadDistribution(original, steps, resolution, parameters=2)
|
|
ret += fts.shortname + " & "
|
|
ret += str(fts.order) + " & "
|
|
ret += str(round(Measures.crps(original, densities1), 3)) + " & "
|
|
ret += str(round(Measures.crps(original, densities2), 3)) + " \\\\ \n"
|
|
return ret
|
|
|
|
|
|
def plotComparedIntervalsAhead(original, models, colors, distributions, time_from, time_to,
|
|
interpol=False, save=False, file=None, tam=[20, 5], resolution=None,
|
|
cmap='Blues',option=2):
|
|
fig = plt.figure(figsize=tam)
|
|
ax = fig.add_subplot(111)
|
|
|
|
cm = plt.get_cmap(cmap)
|
|
cNorm = pltcolors.Normalize(vmin=0, vmax=1)
|
|
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm)
|
|
|
|
if resolution is None: resolution = (max(original) - min(original)) / 100
|
|
|
|
mi = []
|
|
ma = []
|
|
|
|
for count, fts in enumerate(models, start=0):
|
|
if fts.hasDistributionForecasting and distributions[count]:
|
|
density = fts.forecastAheadDistribution(original[time_from - fts.order:time_from],
|
|
time_to, resolution, parameters=option)
|
|
|
|
Y = []
|
|
X = []
|
|
C = []
|
|
S = []
|
|
y = density.columns
|
|
t = len(y)
|
|
|
|
ss = time_to ** 2
|
|
|
|
for k in density.index:
|
|
#alpha = [scalarMap.to_rgba(density[col][k]) for col in density.columns]
|
|
col = [density[col][k]*5 for col in density.columns]
|
|
|
|
x = [time_from + k for x in np.arange(0, t)]
|
|
|
|
s = [ss for x in np.arange(0, t)]
|
|
|
|
ic = resolution/10
|
|
|
|
for cc in np.arange(0, resolution, ic):
|
|
Y.append(y + cc)
|
|
X.append(x)
|
|
C.append(col)
|
|
S.append(s)
|
|
|
|
Y = np.hstack(Y)
|
|
X = np.hstack(X)
|
|
C = np.hstack(C)
|
|
S = np.hstack(S)
|
|
|
|
s = ax.scatter(X, Y, c=C, marker='s',s=S, linewidths=0, edgecolors=None, cmap=cmap)
|
|
s.set_clim([0, 1])
|
|
cb = fig.colorbar(s)
|
|
|
|
cb.set_label('Density')
|
|
|
|
|
|
if fts.hasIntervalForecasting:
|
|
forecasts = fts.forecastAheadInterval(original[time_from - fts.order:time_from], time_to)
|
|
lower = [kk[0] for kk in forecasts]
|
|
upper = [kk[1] for kk in forecasts]
|
|
mi.append(min(lower))
|
|
ma.append(max(upper))
|
|
for k in np.arange(0, time_from - fts.order):
|
|
lower.insert(0, None)
|
|
upper.insert(0, None)
|
|
ax.plot(lower, color=colors[count], label=fts.shortname)
|
|
ax.plot(upper, color=colors[count])
|
|
|
|
else:
|
|
forecasts = fts.forecast(original)
|
|
mi.append(min(forecasts))
|
|
ma.append(max(forecasts))
|
|
for k in np.arange(0, time_from):
|
|
forecasts.insert(0, None)
|
|
ax.plot(forecasts, color=colors[count], label=fts.shortname)
|
|
|
|
ax.plot(original, color='black', label="Original")
|
|
handles0, labels0 = ax.get_legend_handles_labels()
|
|
ax.legend(handles0, labels0, loc=2)
|
|
# ax.set_title(fts.name)
|
|
_mi = min(mi)
|
|
if _mi < 0:
|
|
_mi *= 1.1
|
|
else:
|
|
_mi *= 0.9
|
|
_ma = max(ma)
|
|
if _ma < 0:
|
|
_ma *= 0.9
|
|
else:
|
|
_ma *= 1.1
|
|
|
|
ax.set_ylim([_mi, _ma])
|
|
ax.set_ylabel('F(T)')
|
|
ax.set_xlabel('T')
|
|
ax.set_xlim([0, len(original)])
|
|
|
|
#plt.colorbar()
|
|
|
|
Util.showAndSaveImage(fig, file, save)
|
|
|
|
|
|
def plotCompared(original, forecasts, labels, title):
|
|
fig = plt.figure(figsize=[13, 6])
|
|
ax = fig.add_subplot(111)
|
|
ax.plot(original, color='k', label="Original")
|
|
for c in range(0, len(forecasts)):
|
|
ax.plot(forecasts[c], label=labels[c])
|
|
handles0, labels0 = ax.get_legend_handles_labels()
|
|
ax.legend(handles0, labels0)
|
|
ax.set_title(title)
|
|
ax.set_ylabel('F(T)')
|
|
ax.set_xlabel('T')
|
|
ax.set_xlim([0, len(original)])
|
|
ax.set_ylim([min(original), max(original)])
|
|
|
|
|
|
def SelecaoKFold_MenorRMSE(original, parameters, modelo):
|
|
nfolds = 5
|
|
ret = []
|
|
errors = np.array([[0 for k in parameters] for z in np.arange(0, nfolds)])
|
|
forecasted_best = []
|
|
print("Série Original")
|
|
fig = plt.figure(figsize=[18, 10])
|
|
fig.suptitle("Comparação de modelos ")
|
|
ax0 = fig.add_axes([0, 0.5, 0.65, 0.45]) # left, bottom, width, height
|
|
ax0.set_xlim([0, len(original)])
|
|
ax0.set_ylim([min(original), max(original)])
|
|
ax0.set_title('Série Temporal')
|
|
ax0.set_ylabel('F(T)')
|
|
ax0.set_xlabel('T')
|
|
ax0.plot(original, label="Original")
|
|
min_rmse_fold = 100000.0
|
|
best = None
|
|
fc = 0 # Fold count
|
|
kf = KFold(len(original), n_folds=nfolds)
|
|
for train_ix, test_ix in kf:
|
|
train = original[train_ix]
|
|
test = original[test_ix]
|
|
min_rmse = 100000.0
|
|
best_fold = None
|
|
forecasted_best_fold = []
|
|
errors_fold = []
|
|
pc = 0 # Parameter count
|
|
for p in parameters:
|
|
sets = Grid.GridPartitionerTrimf(train, p)
|
|
fts = modelo(str(p) + " particoes")
|
|
fts.train(train, sets)
|
|
forecasted = [fts.forecast(xx) for xx in test]
|
|
error = Measures.rmse(np.array(forecasted), np.array(test))
|
|
errors_fold.append(error)
|
|
print(fc, p, error)
|
|
errors[fc, pc] = error
|
|
if error < min_rmse:
|
|
min_rmse = error
|
|
best_fold = fts
|
|
forecasted_best_fold = forecasted
|
|
pc = pc + 1
|
|
forecasted_best_fold = [best_fold.forecast(xx) for xx in original]
|
|
ax0.plot(forecasted_best_fold, label=best_fold.name)
|
|
if np.mean(errors_fold) < min_rmse_fold:
|
|
min_rmse_fold = np.mean(errors)
|
|
best = best_fold
|
|
forecasted_best = forecasted_best_fold
|
|
fc = fc + 1
|
|
handles0, labels0 = ax0.get_legend_handles_labels()
|
|
ax0.legend(handles0, labels0)
|
|
ax1 = Axes3D(fig, rect=[0.7, 0.5, 0.3, 0.45], elev=30, azim=144)
|
|
# ax1 = fig.add_axes([0.6, 0.0, 0.45, 0.45], projection='3d')
|
|
ax1.set_title('Comparação dos Erros Quadráticos Médios')
|
|
ax1.set_zlabel('RMSE')
|
|
ax1.set_xlabel('K-fold')
|
|
ax1.set_ylabel('Partições')
|
|
X, Y = np.meshgrid(np.arange(0, nfolds), parameters)
|
|
surf = ax1.plot_surface(X, Y, errors.T, rstride=1, cstride=1, antialiased=True)
|
|
ret.append(best)
|
|
ret.append(forecasted_best)
|
|
|
|
# Modelo diferencial
|
|
print("\nSérie Diferencial")
|
|
errors = np.array([[0 for k in parameters] for z in np.arange(0, nfolds)])
|
|
forecastedd_best = []
|
|
ax2 = fig.add_axes([0, 0, 0.65, 0.45]) # left, bottom, width, height
|
|
ax2.set_xlim([0, len(original)])
|
|
ax2.set_ylim([min(original), max(original)])
|
|
ax2.set_title('Série Temporal')
|
|
ax2.set_ylabel('F(T)')
|
|
ax2.set_xlabel('T')
|
|
ax2.plot(original, label="Original")
|
|
min_rmse = 100000.0
|
|
min_rmse_fold = 100000.0
|
|
bestd = None
|
|
fc = 0
|
|
diff = Transformations.differential(original)
|
|
kf = KFold(len(original), n_folds=nfolds)
|
|
for train_ix, test_ix in kf:
|
|
train = diff[train_ix]
|
|
test = diff[test_ix]
|
|
min_rmse = 100000.0
|
|
best_fold = None
|
|
forecasted_best_fold = []
|
|
errors_fold = []
|
|
pc = 0
|
|
for p in parameters:
|
|
sets = Grid.GridPartitionerTrimf(train, p)
|
|
fts = modelo(str(p) + " particoes")
|
|
fts.train(train, sets)
|
|
forecasted = [fts.forecastDiff(test, xx) for xx in np.arange(len(test))]
|
|
error = Measures.rmse(np.array(forecasted), np.array(test))
|
|
print(fc, p, error)
|
|
errors[fc, pc] = error
|
|
errors_fold.append(error)
|
|
if error < min_rmse:
|
|
min_rmse = error
|
|
best_fold = fts
|
|
pc = pc + 1
|
|
forecasted_best_fold = [best_fold.forecastDiff(original, xx) for xx in np.arange(len(original))]
|
|
ax2.plot(forecasted_best_fold, label=best_fold.name)
|
|
if np.mean(errors_fold) < min_rmse_fold:
|
|
min_rmse_fold = np.mean(errors)
|
|
best = best_fold
|
|
forecasted_best = forecasted_best_fold
|
|
fc = fc + 1
|
|
handles0, labels0 = ax2.get_legend_handles_labels()
|
|
ax2.legend(handles0, labels0)
|
|
ax3 = Axes3D(fig, rect=[0.7, 0, 0.3, 0.45], elev=30, azim=144)
|
|
# ax1 = fig.add_axes([0.6, 0.0, 0.45, 0.45], projection='3d')
|
|
ax3.set_title('Comparação dos Erros Quadráticos Médios')
|
|
ax3.set_zlabel('RMSE')
|
|
ax3.set_xlabel('K-fold')
|
|
ax3.set_ylabel('Partições')
|
|
X, Y = np.meshgrid(np.arange(0, nfolds), parameters)
|
|
surf = ax3.plot_surface(X, Y, errors.T, rstride=1, cstride=1, antialiased=True)
|
|
ret.append(best)
|
|
ret.append(forecasted_best)
|
|
return ret
|
|
|
|
|
|
def SelecaoSimples_MenorRMSE(original, parameters, modelo):
|
|
ret = []
|
|
errors = []
|
|
forecasted_best = []
|
|
print("Série Original")
|
|
fig = plt.figure(figsize=[20, 12])
|
|
fig.suptitle("Comparação de modelos ")
|
|
ax0 = fig.add_axes([0, 0.5, 0.65, 0.45]) # left, bottom, width, height
|
|
ax0.set_xlim([0, len(original)])
|
|
ax0.set_ylim([min(original), max(original)])
|
|
ax0.set_title('Série Temporal')
|
|
ax0.set_ylabel('F(T)')
|
|
ax0.set_xlabel('T')
|
|
ax0.plot(original, label="Original")
|
|
min_rmse = 100000.0
|
|
best = None
|
|
for p in parameters:
|
|
sets = Grid.GridPartitionerTrimf(original, p)
|
|
fts = modelo(str(p) + " particoes")
|
|
fts.train(original, sets)
|
|
# print(original)
|
|
forecasted = fts.forecast(original)
|
|
forecasted.insert(0, original[0])
|
|
# print(forecasted)
|
|
ax0.plot(forecasted, label=fts.name)
|
|
error = Measures.rmse(np.array(forecasted), np.array(original))
|
|
print(p, error)
|
|
errors.append(error)
|
|
if error < min_rmse:
|
|
min_rmse = error
|
|
best = fts
|
|
forecasted_best = forecasted
|
|
handles0, labels0 = ax0.get_legend_handles_labels()
|
|
ax0.legend(handles0, labels0)
|
|
ax1 = fig.add_axes([0.7, 0.5, 0.3, 0.45]) # left, bottom, width, height
|
|
ax1.set_title('Comparação dos Erros Quadráticos Médios')
|
|
ax1.set_ylabel('RMSE')
|
|
ax1.set_xlabel('Quantidade de Partições')
|
|
ax1.set_xlim([min(parameters), max(parameters)])
|
|
ax1.plot(parameters, errors)
|
|
ret.append(best)
|
|
ret.append(forecasted_best)
|
|
# Modelo diferencial
|
|
print("\nSérie Diferencial")
|
|
difffts = Transformations.differential(original)
|
|
errors = []
|
|
forecastedd_best = []
|
|
ax2 = fig.add_axes([0, 0, 0.65, 0.45]) # left, bottom, width, height
|
|
ax2.set_xlim([0, len(difffts)])
|
|
ax2.set_ylim([min(difffts), max(difffts)])
|
|
ax2.set_title('Série Temporal')
|
|
ax2.set_ylabel('F(T)')
|
|
ax2.set_xlabel('T')
|
|
ax2.plot(difffts, label="Original")
|
|
min_rmse = 100000.0
|
|
bestd = None
|
|
for p in parameters:
|
|
sets = Grid.GridPartitionerTrimf(difffts, p)
|
|
fts = modelo(str(p) + " particoes")
|
|
fts.train(difffts, sets)
|
|
forecasted = fts.forecast(difffts)
|
|
forecasted.insert(0, difffts[0])
|
|
ax2.plot(forecasted, label=fts.name)
|
|
error = Measures.rmse(np.array(forecasted), np.array(difffts))
|
|
print(p, error)
|
|
errors.append(error)
|
|
if error < min_rmse:
|
|
min_rmse = error
|
|
bestd = fts
|
|
forecastedd_best = forecasted
|
|
handles0, labels0 = ax2.get_legend_handles_labels()
|
|
ax2.legend(handles0, labels0)
|
|
ax3 = fig.add_axes([0.7, 0, 0.3, 0.45]) # left, bottom, width, height
|
|
ax3.set_title('Comparação dos Erros Quadráticos Médios')
|
|
ax3.set_ylabel('RMSE')
|
|
ax3.set_xlabel('Quantidade de Partições')
|
|
ax3.set_xlim([min(parameters), max(parameters)])
|
|
ax3.plot(parameters, errors)
|
|
ret.append(bestd)
|
|
ret.append(forecastedd_best)
|
|
return ret
|
|
|
|
|
|
def compareModelsPlot(original, models_fo, models_ho):
|
|
fig = plt.figure(figsize=[13, 6])
|
|
fig.suptitle("Comparação de modelos ")
|
|
ax0 = fig.add_axes([0, 0, 1, 1]) # left, bottom, width, height
|
|
rows = []
|
|
for model in models_fo:
|
|
fts = model["model"]
|
|
ax0.plot(model["forecasted"], label=model["name"])
|
|
for model in models_ho:
|
|
fts = model["model"]
|
|
ax0.plot(model["forecasted"], label=model["name"])
|
|
handles0, labels0 = ax0.get_legend_handles_labels()
|
|
ax0.legend(handles0, labels0)
|
|
|
|
|
|
def compareModelsTable(original, models_fo, models_ho):
|
|
fig = plt.figure(figsize=[12, 4])
|
|
fig.suptitle("Comparação de modelos ")
|
|
columns = ['Modelo', 'Ordem', 'Partições', 'RMSE', 'MAPE (%)']
|
|
rows = []
|
|
for model in models_fo:
|
|
fts = model["model"]
|
|
error_r = Measures.rmse(model["forecasted"], original)
|
|
error_m = round(Measures.mape(model["forecasted"], original) * 100, 2)
|
|
rows.append([model["name"], fts.order, len(fts.sets), error_r, error_m])
|
|
for model in models_ho:
|
|
fts = model["model"]
|
|
error_r = Measures.rmse(model["forecasted"][fts.order:], original[fts.order:])
|
|
error_m = round(Measures.mape(model["forecasted"][fts.order:], original[fts.order:]) * 100, 2)
|
|
rows.append([model["name"], fts.order, len(fts.sets), error_r, error_m])
|
|
ax1 = fig.add_axes([0, 0, 1, 1]) # left, bottom, width, height
|
|
ax1.set_xticks([])
|
|
ax1.set_yticks([])
|
|
ax1.table(cellText=rows,
|
|
colLabels=columns,
|
|
cellLoc='center',
|
|
bbox=[0, 0, 1, 1])
|
|
sup = "\\begin{tabular}{"
|
|
header = ""
|
|
body = ""
|
|
footer = ""
|
|
|
|
for c in columns:
|
|
sup = sup + "|c"
|
|
if len(header) > 0:
|
|
header = header + " & "
|
|
header = header + "\\textbf{" + c + "} "
|
|
sup = sup + "|} \\hline\n"
|
|
header = header + "\\\\ \\hline \n"
|
|
|
|
for r in rows:
|
|
lin = ""
|
|
for c in r:
|
|
if len(lin) > 0:
|
|
lin = lin + " & "
|
|
lin = lin + str(c)
|
|
|
|
body = body + lin + "\\\\ \\hline \n"
|
|
|
|
return sup + header + body + "\\end{tabular}"
|
|
|
|
|
|
def simpleSearch_RMSE(train, test, model, partitions, orders, save=False, file=None, tam=[10, 15],
|
|
plotforecasts=False, elev=30, azim=144, intervals=False,parameters=None):
|
|
_3d = len(orders) > 1
|
|
ret = []
|
|
errors = np.array([[0 for k in range(len(partitions))] for kk in range(len(orders))])
|
|
forecasted_best = []
|
|
fig = plt.figure(figsize=tam)
|
|
# fig.suptitle("Comparação de modelos ")
|
|
if plotforecasts:
|
|
ax0 = fig.add_axes([0, 0.4, 0.9, 0.5]) # left, bottom, width, height
|
|
ax0.set_xlim([0, len(train)])
|
|
ax0.set_ylim([min(train) * 0.9, max(train) * 1.1])
|
|
ax0.set_title('Forecasts')
|
|
ax0.set_ylabel('F(T)')
|
|
ax0.set_xlabel('T')
|
|
min_rmse = 1000000.0
|
|
best = None
|
|
|
|
for pc, p in enumerate(partitions, start=0):
|
|
|
|
sets = Grid.GridPartitionerTrimf(train, p)
|
|
for oc, o in enumerate(orders, start=0):
|
|
fts = model("q = " + str(p) + " n = " + str(o))
|
|
fts.train(train, sets, o,parameters=parameters)
|
|
if not intervals:
|
|
forecasted = fts.forecast(test)
|
|
if not fts.hasSeasonality:
|
|
error = Measures.rmse(np.array(test[o:]), np.array(forecasted[:-1]))
|
|
else:
|
|
error = Measures.rmse(np.array(test[o:]), np.array(forecasted))
|
|
for kk in range(o):
|
|
forecasted.insert(0, None)
|
|
if plotforecasts: ax0.plot(forecasted, label=fts.name)
|
|
else:
|
|
forecasted = fts.forecastInterval(test)
|
|
error = 1.0 - Measures.rmse_interval(np.array(test[o:]), np.array(forecasted[:-1]))
|
|
errors[oc, pc] = error
|
|
if error < min_rmse:
|
|
min_rmse = error
|
|
best = fts
|
|
forecasted_best = forecasted
|
|
|
|
# print(min_rmse)
|
|
if plotforecasts:
|
|
# handles0, labels0 = ax0.get_legend_handles_labels()
|
|
# ax0.legend(handles0, labels0)
|
|
ax0.plot(test, label="Original", linewidth=3.0, color="black")
|
|
if _3d: ax1 = Axes3D(fig, rect=[0, 1, 0.9, 0.9], elev=elev, azim=azim)
|
|
if not plotforecasts: ax1 = Axes3D(fig, rect=[0, 1, 0.9, 0.9], elev=elev, azim=azim)
|
|
# ax1 = fig.add_axes([0.6, 0.5, 0.45, 0.45], projection='3d')
|
|
if _3d:
|
|
ax1.set_title('Error Surface')
|
|
ax1.set_ylabel('Model order')
|
|
ax1.set_xlabel('Number of partitions')
|
|
ax1.set_zlabel('RMSE')
|
|
X, Y = np.meshgrid(partitions, orders)
|
|
surf = ax1.plot_surface(X, Y, errors, rstride=1, cstride=1, antialiased=True)
|
|
else:
|
|
ax1 = fig.add_axes([0, 1, 0.9, 0.9])
|
|
ax1.set_title('Error Curve')
|
|
ax1.set_ylabel('Number of partitions')
|
|
ax1.set_xlabel('RMSE')
|
|
ax0.plot(errors,partitions)
|
|
ret.append(best)
|
|
ret.append(forecasted_best)
|
|
|
|
# plt.tight_layout()
|
|
|
|
Util.showAndSaveImage(fig, file, save)
|
|
|
|
return ret
|
|
|
|
|
|
def pftsExploreOrderAndPartitions(data,save=False, file=None):
|
|
fig, axes = plt.subplots(nrows=4, ncols=1, figsize=[6, 8])
|
|
data_fs1 = Grid.GridPartitionerTrimf(data, 10)
|
|
mi = []
|
|
ma = []
|
|
|
|
axes[0].set_title('Point Forecasts by Order')
|
|
axes[2].set_title('Interval Forecasts by Order')
|
|
|
|
for order in np.arange(1, 6):
|
|
fts = pfts.ProbabilisticFTS("")
|
|
fts.shortname = "n = " + str(order)
|
|
fts.train(data, data_fs1, order=order)
|
|
point_forecasts = fts.forecast(data)
|
|
interval_forecasts = fts.forecastInterval(data)
|
|
lower = [kk[0] for kk in interval_forecasts]
|
|
upper = [kk[1] for kk in interval_forecasts]
|
|
mi.append(min(lower) * 0.95)
|
|
ma.append(max(upper) * 1.05)
|
|
for k in np.arange(0, order):
|
|
point_forecasts.insert(0, None)
|
|
lower.insert(0, None)
|
|
upper.insert(0, None)
|
|
axes[0].plot(point_forecasts, label=fts.shortname)
|
|
axes[2].plot(lower, label=fts.shortname)
|
|
axes[2].plot(upper)
|
|
|
|
axes[1].set_title('Point Forecasts by Number of Partitions')
|
|
axes[3].set_title('Interval Forecasts by Number of Partitions')
|
|
|
|
for partitions in np.arange(5, 11):
|
|
data_fs = Grid.GridPartitionerTrimf(data, partitions)
|
|
fts = pfts.ProbabilisticFTS("")
|
|
fts.shortname = "q = " + str(partitions)
|
|
fts.train(data, data_fs, 1)
|
|
point_forecasts = fts.forecast(data)
|
|
interval_forecasts = fts.forecastInterval(data)
|
|
lower = [kk[0] for kk in interval_forecasts]
|
|
upper = [kk[1] for kk in interval_forecasts]
|
|
mi.append(min(lower) * 0.95)
|
|
ma.append(max(upper) * 1.05)
|
|
point_forecasts.insert(0, None)
|
|
lower.insert(0, None)
|
|
upper.insert(0, None)
|
|
axes[1].plot(point_forecasts, label=fts.shortname)
|
|
axes[3].plot(lower, label=fts.shortname)
|
|
axes[3].plot(upper)
|
|
|
|
for ax in axes:
|
|
ax.set_ylabel('F(T)')
|
|
ax.set_xlabel('T')
|
|
ax.plot(data, label="Original", color="black", linewidth=1.5)
|
|
handles, labels = ax.get_legend_handles_labels()
|
|
ax.legend(handles, labels, loc=2, bbox_to_anchor=(1, 1))
|
|
ax.set_ylim([min(mi), max(ma)])
|
|
ax.set_xlim([0, len(data)])
|
|
|
|
plt.tight_layout()
|
|
|
|
Util.showAndSaveImage(fig, file, save)
|
|
|