spec/math_extension_spec.rb in distribution-0.5.0 vs spec/math_extension_spec.rb in distribution-0.6.0
- old
+ new
@@ -1,16 +1,195 @@
require File.expand_path(File.dirname(__FILE__)+"/spec_helper.rb")
+include ExampleWithGSL
describe Distribution::MathExtension do
it "binomial coefficient should be correctly calculated" do
n=50
n.times do |k|
Math.binomial_coefficient(n,k).should eq(Math.factorial(n).quo(Math.factorial(k)*Math.factorial(n-k)))
end
end
+
+ it "ChebyshevSeries for :sin should return correct values" do
+ #Math::SIN_CS.evaluate()
+ end
+
+ it "log_1plusx_minusx should return correct values" do
+ # Tests from GSL-1.9
+ Math::Log.log_1plusx_minusx(1.0e-10).should be_within(1e-10).of(-4.999999999666666667e-21)
+ Math::Log.log_1plusx_minusx(1.0e-8).should be_within(1e-10).of(-4.999999966666666917e-17)
+ Math::Log.log_1plusx_minusx(1.0e-4).should be_within(1e-10).of(-4.999666691664666833e-09)
+ Math::Log.log_1plusx_minusx(0.1).should be_within(1e-10).of(-0.004689820195675139956)
+ Math::Log.log_1plusx_minusx(0.49).should be_within(1e-10).of(-0.09122388004263222704)
+
+ Math::Log.log_1plusx_minusx(-0.49).should be_within(1e-10).of(-0.18334455326376559639)
+ Math::Log.log_1plusx_minusx(1.0).should be_within(1e-10).of(Math::LN2 - 1.0)
+ Math::Log.log_1plusx_minusx(-0.99).should be_within(1e-10).of(-3.615170185988091368)
+ end
+
+ it "log_1plusx should return correct values" do
+ # Tests from GSL-1.9
+ Math::Log.log_1plusx(1.0e-10).should be_within(1e-10).of(9.999999999500000000e-11)
+ Math::Log.log_1plusx(1.0e-8).should be_within(1e-10).of(9.999999950000000333e-09)
+ Math::Log.log_1plusx(1.0e-4).should be_within(1e-10).of(0.00009999500033330833533)
+ Math::Log.log_1plusx(0.1).should be_within(1e-10).of(0.09531017980432486004)
+ Math::Log.log_1plusx(0.49).should be_within(1e-10).of(0.3987761199573677730)
+
+ Math::Log.log_1plusx(-0.49).should be_within(1e-10).of(-0.6733445532637655964)
+ Math::Log.log_1plusx(1.0).should be_within(1e-10).of(Math::LN2)
+ Math::Log.log_1plusx(-0.99).should be_within(1e-10).of(-4.605170185988091368)
+ end
+
+ it "log_beta should return correct values" do
+ Math::Beta.log_beta(1.0e-8, 1.0e-8).first.should be_within(1e-10).of(19.113827924512310617)
+ Math::Beta.log_beta(1.0e-8, 0.01).first.should be_within(1e-10).of(18.420681743788563403)
+ Math::Beta.log_beta(1.0e-8, 1.0).first.should be_within(1e-10).of(18.420680743952365472)
+ Math::Beta.log_beta(1.0e-8, 10.0).first.should be_within(1e-10).of(18.420680715662683009)
+ Math::Beta.log_beta(1.0e-8, 1000.0).first.should be_within(1e-10).of(18.420680669107656949)
+ Math::Beta.log_beta(0.1, 0.1).first.should be_within(1e-10).of(2.9813614810376273949)
+ Math::Beta.log_beta(0.1, 1.0).first.should be_within(1e-10).of(2.3025850929940456840)
+ Math::Beta.log_beta(0.1, 100.0).first.should be_within(1e-10).of(1.7926462324527931217)
+ Math::Beta.log_beta(0.1, 1000).first.should be_within(1e-10).of(1.5619821298353164928)
+ Math::Beta.log_beta(1.0, 1.00025).first.should be_within(1e-10).of(-0.0002499687552073570)
+ Math::Beta.log_beta(1.0, 1.01).first.should be_within(1e-10).of(-0.009950330853168082848)
+ Math::Beta.log_beta(1.0, 1000.0).first.should be_within(1e-10).of(-6.907755278982137052)
+ Math::Beta.log_beta(100.0, 100.0).first.should be_within(1e-10).of(-139.66525908670663927)
+ Math::Beta.log_beta(100.0, 1000.0).first.should be_within(1e-10).of(-336.4348576477366051)
+ Math::Beta.log_beta(100.0, 1.0e+8).first.should be_within(1e-10).of(-1482.9339185256447309)
+ end
+
+ it "regularized_beta should return correct values" do
+ Math.regularized_beta(0.0,1.0, 1.0).should be_within(1e-10).of(0.0)
+ Math.regularized_beta(1.0, 1.0, 1.0).should be_within(1e-10).of(1.0)
+ Math.regularized_beta(1.0, 0.1, 0.1).should be_within(1e-10).of(1.0)
+ Math.regularized_beta(0.5, 1.0, 1.0).should be_within(1e-10).of(0.5)
+ Math.regularized_beta(0.5, 0.1, 1.0).should be_within(1e-10).of(0.9330329915368074160)
+ Math.regularized_beta(0.5, 10.0, 1.0).should be_within(1e-10).of(0.0009765625000000000000)
+ Math.regularized_beta(0.5, 50.0, 1.0).should be_within(1e-10).of(8.881784197001252323e-16)
+ Math.regularized_beta(0.5, 1.0, 0.1).should be_within(1e-10).of(0.06696700846319258402)
+ Math.regularized_beta(0.5, 1.0, 10.0).should be_within(1e-10).of(0.99902343750000000000)
+ Math.regularized_beta(0.5, 1.0, 50.0).should be_within(1e-10).of(0.99999999999999911180)
+ Math.regularized_beta(0.1, 1.0, 1.0).should be_within(1e-10).of(0.10)
+ Math.regularized_beta(0.1, 1.0, 2.0).should be_within(1e-10).of(0.19)
+ Math.regularized_beta(0.9, 1.0, 2.0).should be_within(1e-10).of(0.99)
+ Math.regularized_beta(0.5, 50.0, 60.0).should be_within(1e-10).of(0.8309072939016694143)
+ Math.regularized_beta(0.5, 90.0, 90.0).should be_within(1e-10).of(0.5)
+ Math.regularized_beta(0.5, 500.0, 500.0).should be_within(1e-10).of(0.5)
+ Math.regularized_beta(0.4, 5000.0, 5000.0).should be_within(1e-10).of(4.518543727260666383e-91)
+ Math.regularized_beta(0.6, 5000.0, 5000.0).should be_within(1e-10).of(1.0)
+ Math.regularized_beta(0.6, 5000.0, 2000.0).should be_within(1e-10).of(8.445388773903332659e-89)
+ end
+ it_only_with_gsl "incomplete_beta should return correct values" do
+
+ a=rand()*10+1
+ b=rand()*10+1
+ ib = GSL::Function.alloc { |t| t**(a-1)*(1-t)**(b-1)}
+ w = GSL::Integration::Workspace.alloc(1000)
+ 1.upto(10) {|x|
+ inte=ib.qag([0,x / 10.0],w)
+ Math.incomplete_beta(x/10.0, a ,b).should be_within(1e-10).of(inte[0])
+ }
+ end
+
+ it "gammastar should return correct values" do
+ # Tests from GSL-1.9
+ Math::Gammastar.evaluate(1.0e-08).should be_within(1e-10).of(3989.423555759890865)
+ Math::Gammastar.evaluate(1.0e-05).should be_within(1e-10).of(126.17168469882690233)
+ Math::Gammastar.evaluate(0.001).should be_within(1e-10).of(12.708492464364073506)
+ Math::Gammastar.evaluate(1.5).should be_within(1e-10).of(1.0563442442685598666)
+ Math::Gammastar.evaluate(3.0).should be_within(1e-10).of(1.0280645179187893045)
+ Math::Gammastar.evaluate(9.0).should be_within(1e-10).of(1.0092984264218189715)
+ Math::Gammastar.evaluate(11.0).should be_within(1e-10).of(1.0076024283104962850)
+ Math::Gammastar.evaluate(100.0).should be_within(1e-10).of(1.0008336778720121418)
+ Math::Gammastar.evaluate(1.0e+05).should be_within(1e-10).of(1.0000008333336805529)
+ Math::Gammastar.evaluate(1.0e+20).should be_within(1e-10).of(1.0)
+ end
+
+ it "erfc_e should return correct values" do
+ # From GSL-1.9. For troubleshooting gammq.
+ Math::erfc_e(-10.0).should be_within(1e-10).of(2.0)
+ Math::erfc_e(-5.0000002).should be_within(1e-10).of(1.9999999999984625433)
+ Math::erfc_e(-5.0).should be_within(1e-10).of(1.9999999999984625402)
+ Math::erfc_e(-1.0).should be_within(1e-10).of(1.8427007929497148693)
+ Math::erfc_e(-0.5).should be_within(1e-10).of(1.5204998778130465377)
+ Math::erfc_e(1.0).should be_within(1e-10).of(0.15729920705028513066)
+ Math::erfc_e(3.0).should be_within(1e-10).of(0.000022090496998585441373)
+ Math::erfc_e(7.0).should be_within(1e-10).of(4.183825607779414399e-23)
+ Math::erfc_e(10.0).should be_within(1e-10).of(2.0884875837625447570e-45)
+ end
+
+
+ it "unnormalized_incomplete_gamma with x=0 should return correct values" do
+ Math.unnormalized_incomplete_gamma(-1.5, 0).should be_within(1e-10).of(4.0*Math.sqrt(Math::PI)/3.0)
+ Math.unnormalized_incomplete_gamma(-0.5, 0).should be_within(1e-10).of(-2*Math.sqrt(Math::PI))
+ Math.unnormalized_incomplete_gamma(0.5, 0).should be_within(1e-10).of(Math.sqrt(Math::PI))
+ Math.unnormalized_incomplete_gamma(1.0, 0).should eq 1.0
+ Math.unnormalized_incomplete_gamma(1.5, 0).should be_within(1e-10).of(Math.sqrt(Math::PI)/2.0)
+ Math.unnormalized_incomplete_gamma(2.0, 0).should eq 1.0
+ Math.unnormalized_incomplete_gamma(2.5, 0).should be_within(1e-10).of(0.75*Math.sqrt(Math::PI))
+ Math.unnormalized_incomplete_gamma(3.0, 0).should eq 2.0
+ Math.unnormalized_incomplete_gamma(3.5, 0).should be_within(1e-10).of(15.0*Math.sqrt(Math::PI) / 8.0)
+ Math.unnormalized_incomplete_gamma(4.0, 0).should eq 6.0
+ end
+
+ it "incomplete_gamma should return correct values" do
+ # Tests from GSL-1.9
+ Math.incomplete_gamma(1e-100, 0.001).should be_within(1e-10).of(1.0)
+ Math.incomplete_gamma(0.001, 0.001).should be_within(1e-10).of(0.9936876467088602902)
+ Math.incomplete_gamma(0.001, 1.0).should be_within(1e-10).of(0.9997803916424144436)
+ Math.incomplete_gamma(0.001, 10.0).should be_within(1e-10).of(0.9999999958306921828)
+ Math.incomplete_gamma(1.0, 0.001).should be_within(1e-10).of(0.0009995001666250083319)
+ Math.incomplete_gamma(1.0, 1.01).should be_within(1e-10).of(0.6357810204284766802)
+ Math.incomplete_gamma(1.0, 10.0).should be_within(1e-10).of(0.9999546000702375151)
+ Math.incomplete_gamma(10.0, 10.01).should be_within(1e-10).of(0.5433207586693410570)
+ Math.incomplete_gamma(10.0, 20.0).should be_within(1e-10).of(0.9950045876916924128)
+ Math.incomplete_gamma(1000.0, 1000.1).should be_within(1e-10).of(0.5054666401440661753)
+ Math.incomplete_gamma(1000.0, 2000.0).should be_within(1e-10).of(1.0)
+
+ # designed to trap the a-x=1 problem
+ # These next two are 1e-7 because they give the same output as GSL, but GSL is apparently not totally accurate here.
+ # It's a problem with log_1plusx_mx (log_1plusx_minusx in my code)
+ Math.incomplete_gamma(100, 99.0).should be_within(1e-7).of(0.4733043303994607)
+ Math.incomplete_gamma(200, 199.0).should be_within(1e-7).of(0.4811585880878718)
+
+ # Test for x86 cancellation problems
+ Math.incomplete_gamma(5670, 4574).should be_within(1e-10).of(3.063972328743934e-55)
+ end
+
+ it "gammq should return correct values" do
+ # Tests from GSL-1.9
+ Math.gammq(0.0, 0.001).should be_within(1e-10).of(0.0)
+ Math.gammq(0.001, 0.001).should be_within(1e-10).of(0.006312353291139709793)
+ Math.gammq(0.001, 1.0).should be_within(1e-10).of(0.00021960835758555639171)
+ Math.gammq(0.001, 2.0).should be_within(1e-10).of(0.00004897691783098147880)
+ Math.gammq(0.001, 5.0).should be_within(1e-10).of(1.1509813397308608541e-06)
+ Math.gammq(1.0, 0.001).should be_within(1e-10).of(0.9990004998333749917)
+ Math.gammq(1.0, 1.01).should be_within(1e-10).of(0.3642189795715233198)
+ Math.gammq(1.0, 10.0).should be_within(1e-10).of(0.00004539992976248485154)
+ Math.gammq(10.0, 10.01).should be_within(1e-10).of(0.4566792413306589430)
+ Math.gammq(10.0, 100.0).should be_within(1e-10).of(1.1253473960842733885e-31)
+ Math.gammq(1000.0, 1000.1).should be_within(1e-10).of(0.4945333598559338247)
+ Math.gammq(1000.0, 2000.0).should be_within(1e-10).of(6.847349459614753180e-136)
+
+ # designed to trap the a-x=1 problem
+ Math.gammq(100, 99.0).should be_within(1e-10).of(0.5266956696005394)
+ Math.gammq(200, 199.0).should be_within(1e-10).of(0.5188414119121281)
-
+ # Test for x86 cancellation problems
+ Math.gammq(5670, 4574).should be_within(1e-10).of(1.0000000000000000)
+
+
+ # test suggested by Michel Lespinasse [gsl-discuss Sat, 13 Nov 2004]
+ Math.gammq(1.0e+06-1.0, 1.0e+06-2.0).should be_within(1e-10).of(0.50026596175224547004)
+
+ # tests in asymptotic regime related to Lespinasse test
+ Math.gammq(1.0e+06+2.0, 1.0e+06+1.0).should be_within(1e-10).of(0.50026596135330304336)
+ Math.gammq(1.0e+06, 1.0e+06-2.0).should be_within(1e-10).of(0.50066490399940144811)
+ Math.gammq(1.0e+07, 1.0e+07-2.0).should be_within(1e-10).of(0.50021026104978614908)
+ end
+
+
it "rising_factorial should return correct values" do
x=rand(10)+1
Math.rising_factorial(x,0).should eq 1
Math.rising_factorial(x,1).should eq x
@@ -28,27 +207,23 @@
Math.permutations(n,n).should eq(Math.factorial(n) / Math.factorial(n-n))
end
- it "incomplete beta function should return similar results to R" do
- pending("Not working yet")
- Math.incomplete_beta(0.5,5,6).should be_within(1e-6).of(Math.beta(5,6)*0.6230469)
- Math.incomplete_beta(0.6,5,6).should be_within(1e-6).of(Math.beta(5,6)*0.0006617154)
- end
- it "regularized incomplete beta should behave properly" do
+
+ it "exact regularized incomplete beta should behave properly" do
- Math.regularized_beta_function(0.5,5,5).should eq 0.5
- Math.regularized_beta_function(0.5,5,6).should be_within(1e-6).of(0.6230469)
- Math.regularized_beta_function(0.5,5,7).should be_within(1e-6).of(0.725586)
+ Math.exact_regularized_beta(0.5,5,5).should be_within(1e-6).of(0.5)
+ Math.exact_regularized_beta(0.5,5,6).should be_within(1e-6).of(0.6230469)
+ Math.exact_regularized_beta(0.5,5,7).should be_within(1e-6).of(0.725586)
a=5
b=5
- Math.regularized_beta_function(0,a,b).should eq 0
- Math.regularized_beta_function(1,a,b).should eq 1
+ Math.exact_regularized_beta(0,a,b).should eq 0
+ Math.exact_regularized_beta(1,a,b).should eq 1
x=rand()
- Math.regularized_beta_function(x,a,b).should be_within(1e-6). of(1-Math.regularized_beta_function(1-x,b,a))
+ Math.exact_regularized_beta(x,a,b).should be_within(1e-6). of(1-Math.regularized_beta(1-x,b,a))
end
it "binomial coefficient(gamma) with n<=48 should be correct " do