--- a/probability-discrete.rst	Thu Apr 21 16:18:07 2016 +0300
+++ b/probability-discrete.rst	Thu Apr 21 16:19:13 2016 +0300
@@ -82,7 +82,7 @@
 
 .. math::
 
-  var[X] = E[(X - E[X])^2] = E[X^2] - E^2[X]
+  var[X] = E[(X - E[X])²] = E[X²] - (E[X])²
 
 :def:`Standard deviation` is a:
 
@@ -106,6 +106,8 @@
 
   p_X(x) = Σ_i P(A_i)·p_{X|A_i}(x)
 
+* https://en.wikipedia.org/wiki/Law_of_total_probability
+
 Conditional PMF on event
 ========================
 
@@ -496,3 +498,34 @@
        ------------
         2
        x  - 2 x + 1
+
+.. math::
+
+  E[(geom(p))²] = ∑_{x=1..∞} x²·p·(1-p)^(x-1)
+
+  = p·∑_{x=1..∞} x²·(1-p)^(x-1)
+
+  = p/(1-p)·∑_{x=0..∞} x²·(1-p)^x
+
+  = p/(1-p)·(1-p)·(1-p+1)/(1 - (1-p))³ = p·(2-p)/p³ = (2-p)/p²
+
+.. NOTE::
+
+   Maxima calculation::
+
+    load("simplify_sum");
+    (%i3) assume(x>0);
+    (%o3)                               [x > 0]
+    (%i4) assume(x<1);
+    (%o4)                               [x < 1]
+
+    (%i8) simplify_sum(sum(k^2 * x^k, k, 0, inf));
+                                              2
+                                         x + x
+    (%o8)                        - -------------------
+                                    3      2
+                                   x  - 3 x  + 3 x - 1
+
+So:
+
+.. math:: var(geom(p)) = E[(geom(p))²] - E[geom(p)]² = (2-p)/p² - 1/p² = (1-p)/p²