diff --git a/.editorconfig b/.editorconfig
index e7ef85b..bfea71f 100644
--- a/.editorconfig
+++ b/.editorconfig
@@ -9,7 +9,7 @@ insert_final_newline = true
 charset = utf-8
 
 [*.py]
-indent_style = tab
+indent_style = space
 indent_size = 4
 
 [Makefile]
diff --git a/Makefile b/Makefile
index 6133f0c..1c7a993 100644
--- a/Makefile
+++ b/Makefile
@@ -7,7 +7,7 @@ save:
 	pip freeze > requirements.txt
 
 run:
-	python src/main.py
+	python src/main.py $(MODULE)
 
 test:
 	python -m pytest -v
diff --git a/README.rst b/README.rst
index 1a210cf..c983338 100644
--- a/README.rst
+++ b/README.rst
@@ -2,3 +2,10 @@ Python by example
 =================
 
 The intention of this repository is to teach me the basics, along with the best or common practiced, of python development.
+
+Modules must include a "run" function in order to be executed as a dynamic import.
+
+To run
+======
+
+`$ make MODULE="module submodule" run`
diff --git a/requirements.txt b/requirements.txt
index e5ef09a..73816ef 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,8 +1,10 @@
 iniconfig==2.1.0
 mpmath==1.3.0
+narwhals==2.25.0
 numpy==2.3.3
 packaging==25.0
 pandas==2.3.3
+plotly==6.9.0
 pluggy==1.6.0
 Pygments==2.19.2
 pytest==8.4.2
diff --git a/src/main.py b/src/main.py
index d0d7d39..dcde173 100644
--- a/src/main.py
+++ b/src/main.py
@@ -1,14 +1,15 @@
-from sympy import diff, limit, oo, symbols
-import unittest
+import sys
 
-from modules.essential_math.examples.statistics_example import (
-	normal_distribution_example,
-	normal_distribution_exercise,
-	t_distribution_example,
-	basic_statistic_concepts_example,
-	z_scores_example,
-	final_exercises
-)
+if __name__ == "__main__":
+    if sys.argv.__len__() < 2:
+        print("[ERROR] Module name must be specified as an argument:")
+        print(" $ python src/main.py <module_name>")
+        exit(-1)
 
-if __name__=="__main__":
-	final_exercises()
+    module_name = sys.argv[1]
+    module_path = "modules."+module_name
+    if sys.argv.__len__() > 2:
+        for submodule_name in sys.argv[2:]:
+            module_path += "."+submodule_name
+    module_selected = __import__(module_path, fromlist=[""])
+    module_selected.run()
diff --git a/src/modules/backend/server.py b/src/modules/backend/server.py
new file mode 100644
index 0000000..57e4c9d
--- /dev/null
+++ b/src/modules/backend/server.py
@@ -0,0 +1,2 @@
+def run():
+    print("The server will start here")
diff --git a/src/modules/essential_math/examples/statistics_example.py b/src/modules/essential_math/examples/statistics_example.py
deleted file mode 100644
index 1f0ef06..0000000
--- a/src/modules/essential_math/examples/statistics_example.py
+++ /dev/null
@@ -1,168 +0,0 @@
-from modules.essential_math.statistics import (
-	mean,
-	median,
-	weighted_mean,
-	weighted_mean_inline,
-	population_variance,
-	population_variance_inline,
-	sample_variance,
-	standard_deviation,
-	normal_probability_density_function,
-	normal_cumulative_density_function,
-	inverse_cumulative_density_function,
-	z_score,
-	coeficient_of_variation,
-	test_central_limit_theorem,
-	generic_critical_z_value,
-	margin_of_error,
-	confidence_interval,
-	get_critical_value_range_t,
-)
-
-def basic_statistic_concepts_example():
-	print("=== Statistics module ===")
-
-	list = [ 1, 2, 3, 4, 5, 6]
-	print(">> The mean of {0} is {1}".format(list, mean(list)))
-
-	weights = [0.2, 0.5, 0.7, 1, 0, 0.9]
-	print(">> The weighted_mean of {0} is {1} and it is equivalent to {2}".format(list, weighted_mean(list, weights), weighted_mean_inline(list, weights)))
-
-	print(">> The median is {0}".format(median(list)))
-	
-	values = [ 0, 1, 5, 7, 9, 10, 14]
-	_population_variance = population_variance(values, sum(values) / len(values))
-	population_variance_calc_inline = population_variance_inline(values);
-	print("The population variance is", _population_variance, population_variance_calc_inline)
-
-	std_dev = standard_deviation(values, False)
-	print("The standard deviation is", std_dev)
-
-	sample = values.copy()
-	del sample[3]
-	del sample[1]
-	print("The sample variance for a population is", sample_variance(sample))
-	print("The standard deviation for a population is", standard_deviation(sample, True))
-
-def normal_distribution_example():
-	print("== Normal distribution ==")
-
-	values = [ 0, 1, 5, 7, 9, 10, 14]
-	mean = sum(values) / len(values)
-	std_dev = standard_deviation(values, False)
-	target_x = 1
-	
-	print(">> The probability_density_function for x = 1 over the example data is {0}".format(normal_probability_density_function(target_x, mean, std_dev)))
-
-	print(">> The probability for observing a value smaller than 1 is given by the cumulative density function and it is: {0}".format(normal_cumulative_density_function(target_x, mean, std_dev)))
-
-	target_probability = 0.5
-	expected_value = inverse_cumulative_density_function(target_probability, mean, std_dev);
-	print(">> For a probability of .5 we expect the value: ", expected_value)
-
-def normal_distribution_exercise():
-	# Population with cold MEAN recovery time of 18 days, with std_dev of 1.5 days. 
-	# Chances of recovery between 15 and 21 days
-	mean = 18
-	std_dev = 1.5
-	init = 15
-	end = 21
-	chances = normal_cumulative_density_function(end, mean, std_dev) - normal_cumulative_density_function(init, mean, std_dev)
-	print("Chances of recovering from a cold between 15 and 21 days: ", chances)
-	print("Chances of recovering before 15 days or after 21: ", 1.0 - chances)
-	# since its a normal distribution, the chances are equaly distributed
-	print("Chances of recovering before 15 days: ", (1.0 - chances) / 2)
-
-	# Apply rug (or drug) to 40 people and see a 16 MEAN recovery time. Test if rug improved mean or casuality
-	## One tailed atest: use inverse cdf in order to find the limit value for a given %.
-	new_mean = 16
-	min_target_percentage = 0.05 # this is a standard
-	min_mean = inverse_cumulative_density_function(min_target_percentage, mean, std_dev)
-	if (min_mean < new_mean):
-		print("The rug (drug) did nothing.")
-	else:
-		print("The rug (drug) worked.")
-	## One tailed test with a P value
-	p_value = normal_cumulative_density_function(new_mean, mean, std_dev)
-	if (p_value > min_target_percentage):
-		print("The rug (drug) did nothing.")
-	else:
-		print("The rug (drug) worked.")
-	
-	## Two tailed test (look for both sides of the normal distribution)
-	## Double the checks, harder to prove (x2) and checks if the rug(drug) makes the recovery time worse.
-	left_min_target = min_target_percentage / 2	
-	x1 = inverse_cumulative_density_function(left_min_target, mean, std_dev)
-	x2 = inverse_cumulative_density_function(1.0 - left_min_target, mean, std_dev)
-	if (new_mean < x1 or new_mean > x2):
-		print("The rug (drug) worked.", x1, x2)
-	else:
-		print("The rug (drug) did nothing.", x1, x2)
-	## Two tailed test with a P value
-	p1 = normal_cumulative_density_function(new_mean, mean, std_dev)
-	right_symetrical_mean = mean + (mean - new_mean)
-	p2 = 1.0 - normal_cumulative_density_function(right_symetrical_mean, mean, std_dev)
-	p_value = p1 + p2
-	if (p_value < min_target_percentage):
-		print("The rug (drug) worked.", p_value)
-	else:
-		print("The rug (drug) did nothing.", p_value)
-	#### CONCEPT: P-hacking, searching for data (in big data scenarios) that passes the p_value < 0.05 test and claiming for a relation.
-
-def z_scores_example():
-	print("== Z-scores ==")
-
-	print("A house (A) of 150K in a neighborhood of 140K mean and 3K std_dev has a Z-score: {0}".format(z_score(150000, 140000, 3000)))
-	print("A house (B) of 815K in a neighborhood of 800K mean and 10K std_dev has a Z-score: {0}".format(z_score(815000, 800000, 10000)))
-	print("The House A is much more expensive because its z-score is higher.")
-	print("The neighborhood of B has a coeficient of variation: {0}, and the one of A: {1}".format(coeficient_of_variation(3000, 140000), coeficient_of_variation(10000, 800000)))
-	print("This means that the neighborhood of A has more spread in its prices")
-
-def central_limit_theorem_example():
-	## Central limit theorem
-	test_central_limit_theorem(sample_size=1, sample_count=1000)
-	test_central_limit_theorem(sample_size=31, sample_count=1000)
-
-def t_distribution_example():
-	confidence = 0.95
-	sample_size = 25
-	(lower, upper) = get_critical_value_range_t(confidence, sample_size)
-	print("The confidence interval is: ", lower, upper)
-
-def final_exercises():
-	# 1.
-	pool_widths = (1.78, 1.75, 1.72, 1.74, 1.77)
-	pool_width_mean = mean(pool_widths)
-	pool_width_std = standard_deviation(pool_widths, True)
-	print("1: ", pool_width_mean, pool_width_std)
-
-	# 2. 
-	z_mean = 42
-	z_std_dev = 8
-	z_prob_init = 20
-	z_prob_end = 30
-	z_prob_final = normal_cumulative_density_function(z_prob_end,z_mean,z_std_dev) - normal_cumulative_density_function(z_prob_init, z_mean, z_std_dev)
-	print("2: ", z_prob_final)
-
-	# 3.
-	filament_value = 1.75
-	filament_sample_size = 34
-	filament_mean = 1.715588
-	filament_std_dev = 0.029252
-	filament_percentage_conficence = .99
-	filament_z_value = z_score(filament_value,filament_mean, filament_std_dev)
-	(filament_confidence_init, filament_conficence_end) = confidence_interval(filament_percentage_conficence, filament_sample_size, filament_std_dev, filament_z_value, filament_mean)
-	print("3: ", filament_confidence_init, filament_conficence_end)
-
-	# 4.
-	original_sales_average = 10345
-	original_sales_std_dev = 552
-	new_sales_average = 11641
-	min_sales_percentage = 0.05
-
-	sales_p1 = 1.0 - normal_cumulative_density_function(new_sales_average, original_sales_average, original_sales_std_dev)
-	sales_p = sales_p1 * 2 # take advantage of symmetry
-	if (sales_p < min_sales_percentage):
-		print("The sales campaing worked", sales_p)
-	else:
-		print("The sales campaing did NOT work")
diff --git a/src/modules/essential_math/math.py b/src/modules/essential_math/math.py
index c86052d..4ceb615 100644
--- a/src/modules/essential_math/math.py
+++ b/src/modules/essential_math/math.py
@@ -1,6 +1,6 @@
-## This module represents the first chapter of the book 
+## This module represents the first chapter of the book
 ##  "Essential Math for Data Science" - Thomas Nield
-##  Chaper 1 - Basic Math and Calculus Review 
+##  Chaper 1 - Basic Math and Calculus Review
 
 from cmath import log as complex_log  # used for complex numbers
 from math import e, exp, log
@@ -90,7 +90,7 @@ def t_approximate_integral(f, init, end, precission):
 def t_calculate_integral(f, init, end, symbol):
     return integrate(f, (symbol, init, end))
 
-def test_math_module():
+def run():
 	print("=== Math module ===")
 	t_exponent(2,8)
 	print(t_compound_interest(100, 20 / 100, 2, 12))
diff --git a/src/modules/essential_math/probability.py b/src/modules/essential_math/probability.py
index aa45717..261ed08 100644
--- a/src/modules/essential_math/probability.py
+++ b/src/modules/essential_math/probability.py
@@ -1,6 +1,6 @@
-## This module represents the second chapter of the book 
+## This module represents the second chapter of the book
 ##  "Essential Math for Data Science" - Thomas Nield
-##  Chapter 2 - Probability 
+##  Chapter 2 - Probability
 
 from scipy.stats import binom, beta
 from math import factorial
@@ -31,7 +31,7 @@ class BinomialDistribution:
 		# For each number calc the probability of that exact number of outcomes (no order)
 		for k in range(n + 1):
 			# 1. Simple combinatory with the binomial coeficient (combinations of k elements out of a pool of n without repetition without order)
-			combinatory = BinomialDistribution.binomial_coeficient(n, k)	
+			combinatory = BinomialDistribution.binomial_coeficient(n, k)
 			# 2. Probability of success, the probability of making it k times
 			probability_of_success = p ** k
 			# 3. Probability of failure, inverse of the success
@@ -44,14 +44,14 @@ class BinomialDistribution:
 # > Returns a continuous function, so the final probability of X or better must be calculated using integrals (the area under the curve)
 class BetaDistribution:
 	@staticmethod
-	def calc(probability, success_count, failure_count):	
+	def calc(probability, success_count, failure_count):
 		# Only calcs the rpbability to the left
 		return beta.cdf(probability, success_count, failure_count)
 
 	@staticmethod
 	def calc_right(probability, success_count, failure_count):
 		return 1.0 - BetaDistribution.calc(probability, success_count, failure_count)
-	
+
 	@staticmethod
 	def calc_region(init_probability, end_probability, success_count, failure_count):
 		return BetaDistribution.calc(end_probability, success_count, failure_count) - BetaDistribution.calc(init_probability, success_count, failure_count)
@@ -63,7 +63,7 @@ class Exercises:
 	def one():
 		# 30% change of rain and 40% change your umbrella order will arrive. P(R AND U)
 		return p_r * p_u
-	
+
 	@staticmethod
 	def two():
 		# Same Ps as previous. P(!R OR U)
@@ -81,7 +81,7 @@ class Exercises:
 		p_bail = 0.4
 		p_at_least_50_bail = 0.0
 		for x in range(50, n + 1):
-			p_at_least_50_bail += binom.pmf(x, n, p_bail) 
+			p_at_least_50_bail += binom.pmf(x, n, p_bail)
 		return p_at_least_50_bail
 
 	@staticmethod
@@ -91,7 +91,7 @@ class Exercises:
 		t = 2
 		return  1.0 - BetaDistribution.calc(0.5, 8, 2)
 
-	
+
 	@staticmethod
 	def test():
 		print("P(R AND U) = {0}".format(Exercises.one()))
@@ -101,7 +101,7 @@ class Exercises:
 		print("P(fair) = {0}".format(Exercises.five()))
 
 
-def test_probability_module():
+def run():
 	print("=== Probability module ===")
 	print(">> Binomial distribution")
 	BinomialDistribution.example()
diff --git a/src/modules/essential_math/statistics.py b/src/modules/essential_math/statistics.py
index 77ae442..aefc9d2 100644
--- a/src/modules/essential_math/statistics.py
+++ b/src/modules/essential_math/statistics.py
@@ -1,6 +1,6 @@
-## This module represents the third chapter of the book 
+## This module represents the third chapter of the book
 ##  "Essential Math for Data Science" - Thomas Nield
-##  Chapter 3 - Statistics 
+##  Chapter 3 - Statistics
 
 from math import sqrt, pi, e, exp
 from scipy.stats import norm, t
@@ -8,101 +8,268 @@ from scipy.stats import norm, t
 import random
 import plotly.express as px
 
-def mean(list):
-	return sum(list) / len(list)
+
+def my_mean(list):
+    return sum(list) / len(list)
+
 
 def weighted_mean(items, weights):
-	if (len(items) != len(weights)):
-		return
-	total = 0
-	for i in range(len(items)):
-		total += items[i] * weights[i]
-	return total / sum(weights)
+    if (len(items) != len(weights)):
+        return
+    total = 0
+    for i in range(len(items)):
+        total += items[i] * weights[i]
+    return total / sum(weights)
+
 
 def weighted_mean_inline(items, weights):
-	return sum(s * w for s, w in zip(items, weights)) / sum(weights)
+    return sum(s * w for s, w in zip(items, weights)) / sum(weights)
+
+
 # also called 50% quantile
 def median(items):
-	ordered = sorted(items)
-	length = len(ordered)
-	pair = length % 2 == 0
-	mid = int(length / 2) - 1 if pair else int(n/2)
+    ordered = sorted(items)
+    length = len(ordered)
+    pair = length % 2 == 0
+    mid = int(length / 2) - 1 if pair else int(length/2)
+
+    if pair:
+        return (ordered[mid] + ordered[mid+1]) / 2
+    else:
+        return ordered[mid]
 
-	if pair:
-		return (ordered[mid] + ordered[mid+1]) / 2
-	else:
-		return ordered[mid]
 
 def mode(items):
-	sums = []
+    sums = []
+
 
 def population_variance(value_list, mean):
-	summatory = 0.0
-	for value in value_list:
-		summatory += (value - mean) ** 2
-	return summatory / len(value_list)
+    summatory = 0.0
+    for value in value_list:
+        summatory += (value - mean) ** 2
+    return summatory / len(value_list)
+
 
 def population_variance_inline(value_list):
-	return sum((v - (sum(value_list) / len(value_list))) ** 2 for v in value_list) / len(value_list)
+    return sum((v - (sum(value_list) / len(value_list))) ** 2 for v in value_list) / len(value_list)
+
 
 def sample_variance(value_list):
-	mean = sum(value_list) / len(value_list)
-	return sum((value - mean) ** 2 for value in value_list) / (len(value_list) - 1)
+    mean = sum(value_list) / len(value_list)
+    return sum((value - mean) ** 2 for value in value_list) / (len(value_list) - 1)
+
 
 def population_standard_deviation(value_list):
-	return sqrt(population_variance_inline(value_list))
+    return sqrt(population_variance_inline(value_list))
+
 
 def sample_standard_deviation(value_list):
-	return sqrt(sample_variance(value_list))
+    return sqrt(sample_variance(value_list))
+
 
 def standard_deviation(value_list, is_sample):
-	return sample_standard_deviation(value_list) if is_sample else population_standard_deviation(value_list)
+    return sample_standard_deviation(value_list) if is_sample else population_standard_deviation(value_list)
 
-## Normal distribution
+
+# Normal distribution
 # PDF generates the Normal Distribution (symetric arround the mean)
 def normal_probability_density_function(x: float, mean: float, standard_deviation: float):
-	return (1.0 / (2.0 * pi * standard_deviation ** 2) ** 0.5) * exp(-1.0 * ((x - mean) ** 2 / (2.0 * standard_deviation ** 2)))
+    return (1.0 / (2.0 * pi * standard_deviation ** 2) ** 0.5) * exp(-1.0 * ((x - mean) ** 2 / (2.0 * standard_deviation ** 2)))
+
 
 def normal_cumulative_density_function(x, mean, std_deviation):
-	return norm.cdf(x, mean, std_deviation)
+    return norm.cdf(x, mean, std_deviation)
+
 
 # Check exected value for a given probability
 def inverse_cumulative_density_function(prob, mean, std_dev):
-	x = norm.ppf(prob, mean, std_dev)
-	return x
+    x = norm.ppf(prob, mean, std_dev)
+    return x
+
 
 # Z-scores are valuable in order to normalize 2 pieces of data
 def z_score(value, data_mean, std_deviation):
-	return (value - data_mean) / std_deviation
+    return (value - data_mean) / std_deviation
+
 
 def coeficient_of_variation(std_deviation, mean):
-	return (std_deviation / mean)
+    return (std_deviation / mean)
+
 
 def test_central_limit_theorem(sample_size, sample_count):
-	x_values = [(sum([random.uniform(0.0,1.0) for i in range(sample_size)]) / sample_size) for _ in range(sample_count)]
-	y_values = [1 for _ in range(sample_count)]
-	px.histogram(x=x_values, y=y_values, nbins=20).show()
+    x_values = [(sum([random.uniform(0.0,1.0) for i in range(sample_size)]) / sample_size) for _ in range(sample_count)]
+    y_values = [1 for _ in range(sample_count)]
+    px.histogram(x=x_values, y=y_values, nbins=20).show()
+
 
 def generic_critical_z_value(probability):
-	norm_dist = norm(loc=0.0, scale=1.0)
-	left_tail_area = (1.0 - probability) / 2.0
-	upper_area = 1.0 - ((1.0 - probability) / 2.0)
-	return norm_dist.ppf(left_tail_area), norm_dist.ppf(upper_area)
+    norm_dist = norm(loc=0.0, scale=1.0)
+    left_tail_area = (1.0 - probability) / 2.0
+    upper_area = 1.0 - ((1.0 - probability) / 2.0)
+    return norm_dist.ppf(left_tail_area), norm_dist.ppf(upper_area)
+
 
 def margin_of_error(sample_size, standard_deviation, z_value):
-	return z_value * (standard_deviation / sqrt(sample_size)) # +-, we return the one provided by the z_value (tail or upper)
+    return z_value * (standard_deviation / sqrt(sample_size)) # +-, we return the one provided by the z_value (tail or upper)
+
 
 # How confident we are at a population metric given a sample (the interval we are "probability" sure the value will be)
 def confidence_interval(probability, sample_size, standard_deviation, z_value, mean):
-	critical_z = generic_critical_z_value(probability)
-	margin_error = margin_of_error(sample_size, standard_deviation, z_value)
-	return mean + margin_error, mean - margin_error
+    critical_z = generic_critical_z_value(probability)
+    margin_error = margin_of_error(sample_size, standard_deviation, z_value)
+    return mean + margin_error, mean - margin_error
 
-## T Distribution
-## Similar to the normal distribution but made for smaller sample-sizes (30 or less)
-## When we get close to the 31 items, both are identical
+
+# T Distribution
+# Similar to the normal distribution but made for smaller sample-sizes (30 or less)
+# When we get close to the 31 items, both are identical
 def get_critical_value_range_t(conficence_percentage: float, sample_size: int):
-	untrusted_percentage = 1.0 - conficence_percentage	
-	lower = t.ppf(untrusted_percentage / 2, df=sample_size-1)
-	upper = t.ppf(conficence_percentage + (untrusted_percentage / 2), df=sample_size-1)
-	return (lower, upper)
+    untrusted_percentage = 1.0 - conficence_percentage
+    lower = t.ppf(untrusted_percentage / 2, df=sample_size-1)
+    upper = t.ppf(conficence_percentage + (untrusted_percentage / 2), df=sample_size-1)
+    return (lower, upper)
+
+
+def run():
+    print("=== Statistics module ===")
+
+    list = [ 1, 2, 3, 4, 5, 6]
+    print(">> The mean of {0} is {1}".format(list, my_mean(list)))
+
+    weights = [0.2, 0.5, 0.7, 1, 0, 0.9]
+    print(">> The weighted_mean of {0} is {1} and it is equivalent to {2}".format(list, weighted_mean(list, weights), weighted_mean_inline(list, weights)))
+
+    print(">> The median is {0}".format(median(list)))
+
+    values = [ 0, 1, 5, 7, 9, 10, 14]
+    _population_variance = population_variance(values, sum(values) / len(values))
+    population_variance_calc_inline = population_variance_inline(values);
+    print("The population variance is", _population_variance, population_variance_calc_inline)
+
+    std_dev = standard_deviation(values, False)
+    print("The standard deviation is", std_dev)
+
+    sample = values.copy()
+    del sample[3]
+    del sample[1]
+    print("The sample variance for a population is", sample_variance(sample))
+    print("The standard deviation for a population is", standard_deviation(sample, True))
+
+    print("== Normal distribution ==")
+
+    values = [ 0, 1, 5, 7, 9, 10, 14]
+    mean = sum(values) / len(values)
+    std_dev = standard_deviation(values, False)
+    target_x = 1
+
+    print(">> The probability_density_function for x = 1 over the example data is {0}".format(normal_probability_density_function(target_x, mean, std_dev)))
+
+    print(">> The probability for observing a value smaller than 1 is given by the cumulative density function and it is: {0}".format(normal_cumulative_density_function(target_x, mean, std_dev)))
+
+    target_probability = 0.5
+    expected_value = inverse_cumulative_density_function(target_probability, mean, std_dev);
+    print(">> For a probability of .5 we expect the value: ", expected_value)
+
+    # Population with cold MEAN recovery time of 18 days, with std_dev of 1.5 days.
+    # Chances of recovery between 15 and 21 days
+    mean = 18
+    std_dev = 1.5
+    init = 15
+    end = 21
+    chances = normal_cumulative_density_function(end, mean, std_dev) - normal_cumulative_density_function(init, mean, std_dev)
+    print("Chances of recovering from a cold between 15 and 21 days: ", chances)
+    print("Chances of recovering before 15 days or after 21: ", 1.0 - chances)
+    # since its a normal distribution, the chances are equaly distributed
+    print("Chances of recovering before 15 days: ", (1.0 - chances) / 2)
+
+    # Apply rug (or drug) to 40 people and see a 16 MEAN recovery time. Test if rug improved mean or casuality
+    ## One tailed atest: use inverse cdf in order to find the limit value for a given %.
+    new_mean = 16
+    min_target_percentage = 0.05 # this is a standard
+    min_mean = inverse_cumulative_density_function(min_target_percentage, mean, std_dev)
+    if (min_mean < new_mean):
+        print("The rug (drug) did nothing.")
+    else:
+        print("The rug (drug) worked.")
+    ## One tailed test with a P value
+    p_value = normal_cumulative_density_function(new_mean, mean, std_dev)
+    if (p_value > min_target_percentage):
+        print("The rug (drug) did nothing.")
+    else:
+        print("The rug (drug) worked.")
+
+    ## Two tailed test (look for both sides of the normal distribution)
+    ## Double the checks, harder to prove (x2) and checks if the rug(drug) makes the recovery time worse.
+    left_min_target = min_target_percentage / 2
+    x1 = inverse_cumulative_density_function(left_min_target, mean, std_dev)
+    x2 = inverse_cumulative_density_function(1.0 - left_min_target, mean, std_dev)
+    if (new_mean < x1 or new_mean > x2):
+        print("The rug (drug) worked.", x1, x2)
+    else:
+        print("The rug (drug) did nothing.", x1, x2)
+    ## Two tailed test with a P value
+    p1 = normal_cumulative_density_function(new_mean, mean, std_dev)
+    right_symetrical_mean = mean + (mean - new_mean)
+    p2 = 1.0 - normal_cumulative_density_function(right_symetrical_mean, mean, std_dev)
+    p_value = p1 + p2
+    if (p_value < min_target_percentage):
+        print("The rug (drug) worked.", p_value)
+    else:
+        print("The rug (drug) did nothing.", p_value)
+    #### CONCEPT: P-hacking, searching for data (in big data scenarios) that passes the p_value < 0.05 test and claiming for a relation.
+
+    print("== Z-scores ==")
+
+    print("A house (A) of 150K in a neighborhood of 140K mean and 3K std_dev has a Z-score: {0}".format(z_score(150000, 140000, 3000)))
+    print("A house (B) of 815K in a neighborhood of 800K mean and 10K std_dev has a Z-score: {0}".format(z_score(815000, 800000, 10000)))
+    print("The House A is much more expensive because its z-score is higher.")
+    print("The neighborhood of B has a coeficient of variation: {0}, and the one of A: {1}".format(coeficient_of_variation(3000, 140000), coeficient_of_variation(10000, 800000)))
+    print("This means that the neighborhood of A has more spread in its prices")
+
+    print("== Central Limit Theorem ==")
+    test_central_limit_theorem(sample_size=1, sample_count=1000)
+    test_central_limit_theorem(sample_size=31, sample_count=1000)
+
+    print("== T Distribution ==")
+    confidence = 0.95
+    sample_size = 25
+    (lower, upper) = get_critical_value_range_t(confidence, sample_size)
+    print("The confidence interval is: ", lower, upper)
+
+    print("== Final Exercises ==")
+    # 1.
+    pool_widths = (1.78, 1.75, 1.72, 1.74, 1.77)
+    pool_width_mean = my_mean(pool_widths)
+    pool_width_std = standard_deviation(pool_widths, True)
+    print("1: ", pool_width_mean, pool_width_std)
+
+    # 2.
+    z_mean = 42
+    z_std_dev = 8
+    z_prob_init = 20
+    z_prob_end = 30
+    z_prob_final = normal_cumulative_density_function(z_prob_end,z_mean,z_std_dev) - normal_cumulative_density_function(z_prob_init, z_mean, z_std_dev)
+    print("2: ", z_prob_final)
+
+    # 3.
+    filament_value = 1.75
+    filament_sample_size = 34
+    filament_mean = 1.715588
+    filament_std_dev = 0.029252
+    filament_percentage_conficence = .99
+    filament_z_value = z_score(filament_value,filament_mean, filament_std_dev)
+    (filament_confidence_init, filament_conficence_end) = confidence_interval(filament_percentage_conficence, filament_sample_size, filament_std_dev, filament_z_value, filament_mean)
+    print("3: ", filament_confidence_init, filament_conficence_end)
+
+    # 4.
+    original_sales_average = 10345
+    original_sales_std_dev = 552
+    new_sales_average = 11641
+    min_sales_percentage = 0.05
+
+    sales_p1 = 1.0 - normal_cumulative_density_function(new_sales_average, original_sales_average, original_sales_std_dev)
+    sales_p = sales_p1 * 2 # take advantage of symmetry
+    if (sales_p < min_sales_percentage):
+        print("The sales campaing worked", sales_p)
+    else:
+        print("The sales campaing did NOT work")
diff --git a/src/modules/exercises/essentials.py b/src/modules/exercises/essentials.py
index 4959c65..ae49cf6 100644
--- a/src/modules/exercises/essentials.py
+++ b/src/modules/exercises/essentials.py
@@ -1,3 +1,7 @@
+def run():
+    print("! This module is ran through tests, do not expect an output.")
+
+
 def maximum_subarray_sum(input_array: list[int]):
 	max_sum = input_array[0]
 	subarray = [input_array[0]]
@@ -31,7 +35,7 @@ def trap_rain_water(input_aray: list[int]):
 			else:
 				wall.insert(coord_index, "a")
 		container.insert(wall_index, wall)
-	
+
 	total_water = 0
 	# Step 2: fill the air with water
 	for wall_x in range(0, len(input_aray)):
diff --git a/src/modules/primitives/strings.py b/src/modules/primitives/strings.py
index 9463f2a..832b6ce 100644
--- a/src/modules/primitives/strings.py
+++ b/src/modules/primitives/strings.py
@@ -1,3 +1,7 @@
 def t_strings():
-	my_template = "This is a value: {test}"
-	print(my_template.format(test="another"))
+    my_template = "This is a value: {test}"
+    print(my_template.format(test="another"))
+
+
+def run():
+    t_strings()