diff --git a/fermi.go b/fermi.go
index 6029a06..aff812e 100644
--- a/fermi.go
+++ b/fermi.go
@@ -122,14 +122,15 @@ func prettyPrintDist(dist Dist) {
 		pretty.PrettyPrintFloat(w)
 		fmt.Println()
 	case FilledSamples:
-		sorted_xs := make([]float64, N_SAMPLES)
+		n := len(v.xs)
+		sorted_xs := make([]float64, n)
 		copy(sorted_xs, v.xs)
 		sort.Slice(sorted_xs, func(i, j int) bool {
 			return sorted_xs[i] < sorted_xs[j]
 		})
 
-		low := sorted_xs[N_SAMPLES/20]
-		high := sorted_xs[N_SAMPLES*19/20]
+		low := sorted_xs[int(math.Round(float64(n)*0.05))]
+		high := sorted_xs[int(math.Round(float64(n)*0.05))]
 		fmt.Printf("=> ")
 		pretty.PrettyPrint2Floats(low, high)
 
@@ -149,6 +150,50 @@ func printAndReturnErr(err_msg string) error {
 	return errors.New(err_msg)
 }
 
+func prettyPrintStats(dist Dist) {
+	xs := dist.Samples()
+	n := len(xs)
+
+	mean := 0.0
+	for i := 0; i < n; i++ {
+		mean += xs[i]
+	}
+	mean /= float64(n)
+	fmt.Printf("Mean:   %f\n", mean)
+
+	stdev := 0.0
+	for i := 0; i < n; i++ {
+		stdev += math.Pow(xs[i]-mean, 2)
+	}
+	stdev = math.Sqrt(stdev / float64(n))
+	fmt.Printf("Stdev:  %f\n", stdev)
+
+	sorted_xs := make([]float64, n)
+	copy(sorted_xs, xs)
+	sort.Slice(sorted_xs, func(i, j int) bool {
+		return sorted_xs[i] < sorted_xs[j]
+	})
+	ci_01 := sorted_xs[int(math.Round(float64(n)*0.01))]
+	ci_05 := sorted_xs[int(math.Round(float64(n)*0.05))]
+	ci_10 := sorted_xs[int(math.Round(float64(n)*0.10))]
+	ci_25 := sorted_xs[int(math.Round(float64(n)*0.25))]
+	ci_50 := sorted_xs[int(math.Round(float64(n)*0.50))]
+	ci_75 := sorted_xs[int(math.Round(float64(n)*0.75))]
+	ci_90 := sorted_xs[int(math.Round(float64(n)*0.90))]
+	ci_95 := sorted_xs[int(math.Round(float64(n)*0.95))]
+	ci_99 := sorted_xs[int(math.Round(float64(n)*0.99))]
+	fmt.Printf("ci 1%%:  %f\n", ci_01)
+	fmt.Printf("ci 5%%:  %f\n", ci_05)
+	fmt.Printf("ci 10%%: %f\n", ci_10)
+	fmt.Printf("ci 25%%: %f\n", ci_25)
+	fmt.Printf("ci 50%%: %f\n", ci_50)
+	fmt.Printf("ci 75%%: %f\n", ci_75)
+	fmt.Printf("ci 90%%: %f\n", ci_90)
+	fmt.Printf("ci 95%%: %f\n", ci_95)
+	fmt.Printf("ci 99%%: %f\n", ci_99)
+
+}
+
 /* Operations */
 // Generic operations with samples
 func operateDistsAsSamples(dist1 Dist, dist2 Dist, op string) (Dist, error) {
@@ -381,6 +426,8 @@ replForLoop:
 		case words[0] == "clear" || words[0] == "c" || words[0] == ".":
 			stack.old_dist = INIT_DIST
 			fmt.Println()
+		case words[0] == "stats" || words[0] == "s":
+			prettyPrintStats(stack.old_dist)
 		/* Variable assignment */
 		case words[0] == "=:" && len(words) == 2:
 			stack.vars[words[1]] = stack.old_dist