diff --git a/fit_statistics.py b/fit_statistics.py
new file mode 100644
index 0000000..ea825ba
--- /dev/null
+++ b/fit_statistics.py
@@ -0,0 +1,64 @@
+
+from ModelFit import get_best_fit
+import matplotlib.pyplot as plt
+
+import os
+import numpy as np
+
+
+def get_all_cell_folders(dir):
+    paths = []
+    for item in os.listdir(dir):
+        cell_path = os.path.join(dir, item)
+        if os.path.isdir(cell_path):
+            paths.append(cell_path)
+
+    return paths
+
+
+def main():
+
+    values_models = {}
+    values_cells = {}
+
+    for cell in get_all_cell_folders("results/invivo_results/"):
+        best_fit = get_best_fit(cell)
+        c_behaviour, m_behaviour = best_fit.get_behaviour_values()
+
+        for key in c_behaviour.keys():
+            if key not in values_cells.keys():
+                values_cells[key] = []
+            if key not in values_models.keys():
+                values_models[key] = []
+
+            values_cells[key].append(c_behaviour[key])
+            values_models[key].append(m_behaviour[key])
+
+    percentage_error_data = []
+    labels = []
+    for key in values_models.keys():
+        errors = []
+        for i in range(len(values_models[key])):
+            if values_cells[key][i] == 0:
+                if values_models[key][i] == 0:
+                    errors.append(0)
+                else:
+                    print("Cannot calc % error if reference is 0")
+                continue
+            errors.append((values_models[key][i] / values_cells[key][i]) -1)
+
+        percentage_error_data.append(errors)
+        labels.append(key)
+
+    plt.boxplot(percentage_error_data)
+    plt.xticks(range(1, len(percentage_error_data)+1, 1), labels=labels, rotation=45)
+    plt.show()
+    plt.close()
+
+    # plt.plot([np.median(x) for x in percentage_error_data], 'o')
+    # plt.xticks(range(len(percentage_error_data)), labels=labels, rotation=45)
+    # plt.show()
+
+
+if __name__ == '__main__':
+    main()