diff --git a/pylot/core/pick/compare.py b/pylot/core/pick/compare.py
index 156cffd7..5eaccb4c 100644
--- a/pylot/core/pick/compare.py
+++ b/pylot/core/pick/compare.py
@@ -3,6 +3,10 @@
 
 import copy
 
+import numpy as np
+
+import matplotlib.pyplot as plt
+
 from obspy import read_events
 
 from pylot.core.read.io import picks_from_evt
@@ -20,16 +24,18 @@ class Comparison(object):
     compared pick sets. The results can be displayed as histograms showing its
     properties.
     """
+
     def __init__(self, **kwargs):
         names = list()
         self._pdfs = dict()
-        for name, fn in kwargs:
+        for name, fn in kwargs.items():
             self._pdfs[name] = PDFDictionary.from_quakeml(fn)
             names.append(name)
         if len(names) > 2:
             raise ValueError('Comparison is only defined for two '
-                                 'arguments!')
+                             'arguments!')
         self._names = names
+        self._compare = self.compare_picksets()
 
     def __nonzero__(self):
         if not len(self.names) == 2 or not self._pdfs:
@@ -49,10 +55,23 @@ class Comparison(object):
                                                             ' is either not a' \
                                                             ' list or its ' \
                                                             'length is not 2:' \
-                                                            'names : {names}'.format(names=names)
+                                                            'names : {names}'.format(
+            names=names)
         self._names = names
 
-    def compare_picksets(self):
+    @property
+    def comparison(self):
+        return self._compare
+
+    @property
+    def stations(self):
+        return self.comparison.keys()
+
+    @property
+    def nstations(self):
+        return len(self.stations)
+
+    def compare_picksets(self, type='exp'):
         """
         Compare two picksets A and B and return a dictionary compiling the results.
         Comparison is carried out with the help of pdf representation of the picks
@@ -67,8 +86,8 @@ class Comparison(object):
         """
         compare_pdfs = dict()
 
-        pdf_a = self.get(self.names[0])
-        pdf_b = self.get(self.names[1])
+        pdf_a = self.get(self.names[0]).pdf_data(type)
+        pdf_b = self.get(self.names[1]).pdf_data(type)
 
         for station, phases in pdf_a.items():
             if station in pdf_b.keys():
@@ -82,12 +101,42 @@ class Comparison(object):
 
         return compare_pdfs
 
+    def plot(self):
+        nstations = self.nstations
+        stations = self.stations
+        istations = range(nstations)
+        fig, axarr = plt.subplots(nstations, 2, sharex='col', sharey='row')
+
+        for n in istations:
+            station = stations[n]
+            compare_pdf = self.comparison[station]
+            for l, phase in enumerate(compare_pdf.keys()):
+                axarr[n, l].plot(compare_pdf[phase].axis,
+                                 compare_pdf[phase].data)
+                if n is 0:
+                    axarr[n, l].set_title(phase)
+                if l is 0:
+                    axann = axarr[n, l].annotate(station, xy=(.05, .5),
+                                                 xycoords='axes fraction')
+                    bbox_props = dict(boxstyle='round', facecolor='lightgrey',
+                                      alpha=.7)
+                    axann.set_bbox(bbox_props)
+                if n == int(np.median(istations)) and l is 0:
+                    label = 'probability density (qualitative)'
+                    axarr[n, l].set_ylabel(label)
+        plt.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
+        plt.setp([a.get_yticklabels() for a in axarr[:, 1]], visible=False)
+        plt.setp([a.get_yticklabels() for a in axarr[:, 0]], visible=False)
+
+        plt.show()
+
 
 class PDFDictionary(object):
     """
     A PDFDictionary is a dictionary like object containing structured data on
     the probability density function of seismic phase onsets.
     """
+
     def __init__(self, data):
         self._pickdata = data
 
@@ -111,7 +160,7 @@ class PDFDictionary(object):
         if len(cat) > 1:
             raise NotImplementedError('reading more than one event at the same '
                                       'time is not implemented yet! Sorry!')
-        self.pick_data = picks_from_evt(cat[0])
+        return PDFDictionary(picks_from_evt(cat[0]))
 
     def pdf_data(self, type='exp'):
         """
@@ -127,10 +176,10 @@ class PDFDictionary(object):
 
         for station, phases in pdf_picks.items():
             for phase, values in phases.items():
-                phases[phase] = ProbabilityDensityFunction.fromPick(values['epp'],
-                                                                    values['mpp'],
-                                                                    values['lpp'],
-                                                                    type=type)
+                phases[phase] = ProbabilityDensityFunction.fromPick(
+                    values['epp'],
+                    values['mpp'],
+                    values['lpp'],
+                    type=type)
 
         return pdf_picks
-