diff --git a/benchmarks.py b/benchmarks.py
index 57390fc..289f0a5 100644
--- a/benchmarks.py
+++ b/benchmarks.py
@@ -59,7 +59,14 @@ def getIntervalStatistics(original,models):
 		ret = ret + str( round(sharpness(forecasts),2)) + "		& " 
 		ret = ret + str( round(resolution(forecasts),2)) + "		& " 
 		ret = ret + str( round(coverage(original[fts.order-1 :],forecasts),2)) + "	\\ \n"
-	return ret 
+	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='Reds',edgecolors=None)
 
 def plotComparedSeries(original,models, colors):
 	fig = plt.figure(figsize=[25,10])
diff --git a/pifts.py b/pifts.py
index 60613f5..c362b32 100644
--- a/pifts.py
+++ b/pifts.py
@@ -1,4 +1,5 @@
 import numpy as np
+import pandas as pd
 from pyFTS import *
 
 class ProbabilisticFLRG(hofts.HighOrderFLRG):
@@ -79,6 +80,8 @@ class ProbabilisticIntervalFTS(ifts.IntervalFTS):
 		
 		ndata = np.array(data)
 		
+		#print(ndata)
+		
 		l = len(ndata)
 		
 		ret = []
@@ -171,7 +174,46 @@ class ProbabilisticIntervalFTS(ifts.IntervalFTS):
 				ret.append([np.min(lower),np.max(upper)])
 			
 		return ret
+		
+	def getGridClean(self,resolution):
+		grid = {}
+		for sbin in np.arange(self.sets[0].lower,self.sets[-1].upper,resolution):
+			grid[sbin] = 0
 			
+		return grid
+		
+	def gridCount(self, grid, resolution, interval):
+		for sbin in sorted(grid):
+			if sbin >= interval[0] and (sbin + resolution) <= interval[1]:
+				grid[sbin] = grid[sbin] + 1
+		return grid
+		
+	def forecastDistributionAhead(self,data,steps,resolution):
+		
+		ret = []
+		
+		intervals = self.forecastAhead(data,steps)
+		
+		for k in np.arange(self.order,steps):
+			
+			grid = self.getGridClean(resolution)
+			
+			qt1st = self.forecast([intervals[x][0] + (intervals[x][1]-intervals[x][0])/4 for x in np.arange(k-self.order,k)] )
+			qt2nd = self.forecast([intervals[x][0] + (intervals[x][1]-intervals[x][0])/2 for x in np.arange(k-self.order,k)] )
+			qt3rd = self.forecast([intervals[x][1] - (intervals[x][1]-intervals[x][0])/4 for x in np.arange(k-self.order,k)] )
+			grid = self.gridCount(grid,resolution, intervals[k])
+			grid = self.gridCount(grid,resolution, np.ravel(qt1st))
+			grid = self.gridCount(grid,resolution, np.ravel(qt2nd))
+			grid = self.gridCount(grid,resolution, np.ravel(qt3rd))
+			
+			tmp = np.array([ grid[k] for k in sorted(grid) ])
+			
+			ret.append( tmp/sum(tmp) )
+			
+		grid = self.getGridClean(resolution)
+		df = pd.DataFrame(ret, columns=sorted(grid))
+		return df
+		
 	def __str__(self):
 		tmp = self.name + ":\n"
 		for r in sorted(self.flrgs):