pontificate about tests with wide lognormals

README.md

@@ -176,6 +176,66 @@ int main(){
 }
 ```
 
+### Tests and the long tail of the lognormal
+
+Distribution functions can be tested with:
+
+```sh
+cd tests
+make && make run
+```
+
+`make verify` is an alias that runs all the tests and just displays the ones that are failing. 
+
+These tests are somewhat rudimentary: they get between 1M and 10M samples from a given sampling function, and check that their mean and standard deviations correspond to what they should theoretically should be.
+
+If you run `make run` (or `make verify`), you will see errors such as these:
+
+```
+[-] Mean test for normal(47211.047473, 682197.019012) NOT passed.
+Mean of normal(47211.047473, 682197.019012): 46933.673278, vs expected mean: 47211.047473
+delta: -277.374195, relative delta: -0.005910
+
+[-] Std test for lognormal(4.584666, 2.180816) NOT passed.
+Std of lognormal(4.584666, 2.180816): 11443.588861, vs expected std: 11342.434900
+delta: 101.153961, relative delta: 0.008839
+
+[-] Std test for to(13839.861856, 897828.354318) NOT passed.
+Std of to(13839.861856, 897828.354318): 495123.630575, vs expected std: 498075.002499
+delta: -2951.371925, relative delta: -0.005961
+```
+
+These tests I wouldn't worry about. Due to luck of the draw, their relative error is a bit over 0.005, or 0.5%, and so the test fails. But it would surprise me if that had some meaningful practical implication.
+
+The errors that should raise some worry are:
+
+```
+[-] Mean test for lognormal(1.210013, 4.766882) NOT passed.
+Mean of lognormal(1.210013, 4.766882): 342337.257677, vs expected mean: 288253.061628
+delta: 54084.196049, relative delta: 0.157985
+[-] Std test for lognormal(1.210013, 4.766882) NOT passed.
+Std of lognormal(1.210013, 4.766882): 208107782.972184, vs expected std: 24776840217.604111
+delta: -24568732434.631927, relative delta: -118.057730
+
+[-] Mean test for lognormal(-0.195240, 4.883106) NOT passed.
+Mean of lognormal(-0.195240, 4.883106): 87151.733198, vs expected mean: 123886.818303
+delta: -36735.085104, relative delta: -0.421507
+[-] Std test for lognormal(-0.195240, 4.883106) NOT passed.
+Std of lognormal(-0.195240, 4.883106): 33837426.331671, vs expected std: 18657000192.914921
+delta: -18623162766.583248, relative delta: -550.371727
+
+[-] Mean test for lognormal(0.644931, 4.795860) NOT passed.
+Mean of lognormal(0.644931, 4.795860): 125053.904456, vs expected mean: 188163.894101
+delta: -63109.989645, relative delta: -0.504662
+[-] Std test for lognormal(0.644931, 4.795860) NOT passed.
+Std of lognormal(0.644931, 4.795860): 39976300.711166, vs expected std: 18577298706.170452
+delta: -18537322405.459286, relative delta: -463.707799
+```
+
+What is happening in this case is that you are taking a normal, like `normal(-0.195240, 4.883106)`, and you are exponentiating it to arrive at a lognormal. But `normal(-0.195240, 4.883106)` is going to have some noninsignificant weight on, say, 18. But `exp(18) = 39976300`, and points like it are going to end up a nontrivial amount to the analytical mean and standard deviation, even though they have little probability mass.
+
+Fortunately, the reader can also check that for more plausible real-world values, like the 
+
 ## Related projects
 
 - [Squiggle](https://www.squiggle-language.com/)
@@ -185,7 +245,6 @@ int main(){
 
 ## To do list
 
-- [ ] Pontificate about lognormal tests
 - [ ] Have some more complicated & realistic example
 - [ ] Add summarization functions: 90% ci (or all c.i.?) 
 - [ ] Systematize references
@@ -232,3 +291,4 @@ int main(){
   - https://link.springer.com/article/10.1007/bf02293108
   - https://stats.stackexchange.com/questions/502146/how-does-numpy-generate-samples-from-a-beta-distribution
   - https://github.com/numpy/numpy/blob/5cae51e794d69dd553104099305e9f92db237c53/numpy/random/src/distributions/distributions.c
+- [x] Pontificate about lognormal tests