improve the t-distribution inv cdf

chen0040 · May 17, 2017 · c73bd2c · c73bd2c
1 parent b73079f
commit c73bd2c
Show file tree

Hide file tree

Showing 2 changed files with 63 additions and 8 deletions.
diff --git a/src/jsstats.js b/src/jsstats.js
@@ -263,11 +263,58 @@ var jsstats = jsstats || {};
 			df = this.df;
 		}
 		var delta = 0.005;
-        var row = df - 1;
-        var column = Math.round((p - this.percentiles[0]) / delta);
-        column = (column >= 0) ? column : 0;
-        column = (column < this.percentiles.Length) ? column : this.percentiles.length - 1;
-        return (row < this.ttable.length) ? this.ttable[row][column] : this.ntable[column];
+
+        if(p >= 0.5) {
+            var Z1 = 0;
+            for(Z = 0; Z < 100; Z++) {
+                if(this.cumulativeProbability(Z, df) >= p){
+                    break;
+                }
+                Z1 = Z;
+            }
+            var Z2 = Z1;
+            for(var Z = 0.0; Z < 100.0; Z+=1.0) {
+
+                if(this.cumulativeProbability(Z1 + Z / 100.0) >= p){
+                    break;
+                }
+                Z2 = Z1 + (Z)/100.0;
+            }
+            var Z3 = Z2;
+            for(var Z = 0.0; Z < 100.0; Z+=1.0) {
+
+                if(this.cumulativeProbability(Z2 + Z / 10000.0) >= p){
+                    break;
+                }
+                Z3 = Z2 + (Z)/10000.0;
+            }
+            return Z3;
+        } else {
+            var Z1 = 0;
+            for(var Z = 0; Z < 100; Z++) {
+                if(this.cumulativeProbability(-Z, df) <= p){
+                    break;
+                }
+                Z1 = Z;
+            }
+            var Z2 = Z1;
+            for(var Z = 0.0; Z < 100.0; Z+=1.0) {
+
+                if(this.cumulativeProbability(-Z1 - Z / 100.0) <= p){
+                    break;
+                }
+                Z2 = Z1 + (Z) / 100.0;
+            }
+            var Z3 = Z2;
+            for(var Z = 0.0; Z < 100.0; Z+=1.0) {
+
+                if(this.cumulativeProbability(-Z2 - Z / 10000.0) <= p){
+                    break;
+                }
+                Z3 = Z2 + (Z)/10000.0;
+            }
+            return -Z3;
+        }
     };
 
 	jsstats.TDistribution = TDistribution;

diff --git a/test/t-distribution-spec.js b/test/t-distribution-spec.js
@@ -11,14 +11,22 @@ describe("Create t distribution", function() {
   });
 
   describe('run cumulative probability', function(){
+    var distribution = new jsstats.TDistribution(10);
     it('has probability of 0.5 at t_df = 0', function(){
-      var distribution = new jsstats.TDistribution(10.0);
       expect(distribution.cumulativeProbability(0.0)).to.equal(0.5);
     });
     it('has t_df = 0 with probability of 0.5', function(){
-      var distribution = new jsstats.TDistribution(10);
-      expect(distribution.invCumulativeProbability(0.5)).to.equal(0.0);
+      expect(distribution.invCumulativeProbability(0.5)).to.above(-0.001);
+      expect(distribution.invCumulativeProbability(0.5)).to.below(+0.001);
     });
+    it('should run correctly and ascendingly for Z on range p = 0 to p = 1.0', function(){
+      var prevZ = -10000000;
+      for(var p = 0.0; p < 1.0; p += 0.01) {
+        var Z = distribution.invCumulativeProbability(p);
+        expect(Z).to.above(prevZ);
+        prevZ = Z;
+      }
+    })
   });