Permalink
Browse files

Removed Distribution module from statsample and moved to a gem. Chang…

…es on code to reflect new API
  • Loading branch information...
1 parent dd8a267 commit b2c233fa268cbba1143f7ce068d441b2aa163fba @clbustos committed Jan 27, 2011
View
2 Rakefile
@@ -40,7 +40,7 @@ h=Hoe.spec('statsample') do
#self.testlib=:minitest
self.rubyforge_name = "ruby-statsample"
self.developer('Claudio Bustos', 'clbustos@gmail.com')
- self.extra_deps << ["spreadsheet","~>0.6.5"] << ["reportbuilder", "~>1.4"] << ["minimization", "~>0.2.0"] << ["fastercsv", ">0"] << ["dirty-memoize", "~>0.0"] << ["extendmatrix","~>0.3.1"] << ["statsample-bivariate-extension", ">0"] << ["rserve-client", "~>0.2.5"] << ["rubyvis", "~>0.4.0"]
+ self.extra_deps << ["spreadsheet","~>0.6.5"] << ["reportbuilder", "~>1.4"] << ["minimization", "~>0.2.0"] << ["fastercsv", ">0"] << ["dirty-memoize", "~>0.0"] << ["extendmatrix","~>0.3.1"] << ["statsample-bivariate-extension", ">0"] << ["rserve-client", "~>0.2.5"] << ["rubyvis", "~>0.4.0"] << ["distribution", "~>0.2.0"]
self.extra_dev_deps << ["hoe","~>0"] << ["shoulda","~>0"] << ["minitest", "~>2.0"] << ["rserve-client", "~>0"] << ["gettext", "~>0"] << ["mocha", "~>0"] << ["hoe-git", "~>0"]
View
2 benchmarks/correlation_matrix_methods/correlation_matrix.rb
@@ -5,7 +5,7 @@
require 'benchmark'
def create_dataset(vars,cases)
- ran=Distribution::Normal.rng_ugaussian
+ ran=Distribution::Normal.rng
ds=vars.times.inject({}) {|ac,v|
ac["x#{v}"]=Statsample::Vector.new_scale(cases) {ran.call}
ac
View
2 examples/parallel_analysis.rb
@@ -7,7 +7,7 @@
iterations=50
Statsample::Analysis.store(Statsample::Factor::ParallelAnalysis) do
-rng = Distribution::Normal.rng_ugaussian()
+rng = Distribution::Normal.rng()
f1=rnorm(samples)
f2=rnorm(samples)
f3=rnorm(samples)
View
2 examples/velicer_map_test.rb
@@ -5,7 +5,7 @@
Statsample::Analysis.store(Statsample::Factor::MAP) do
- rng=Distribution::Normal.rng_ugaussian
+ rng=Distribution::Normal.rng
samples=100
variables=10
View
26 lib/distribution.rb
@@ -1,26 +0,0 @@
-begin
- require 'statistics2'
-rescue LoadError
- puts "You should install statistics2"
-end
-
-# Several distributions modules to calculate cdf, inverse cdf and pdf
-# See Distribution::Pdf for interface.
-#
-# Usage:
-# Distribution::Normal.cdf(1.96)
-# => 0.97500210485178
-# Distribution::Normal.p_value(0.95)
-# => 1.64485364660836
-module Distribution
- def self.has_gsl?
- Statsample.has_gsl?
- end
-
- autoload(:ChiSquare, 'distribution/chisquare')
- autoload(:T, 'distribution/t')
- autoload(:F, 'distribution/f')
- autoload(:Normal, 'distribution/normal')
- autoload(:NormalBivariate, 'distribution/normalbivariate')
- # autoload(:NormalMultivariate, 'distribution/normalmultivariate')
-end
View
23 lib/distribution/chisquare.rb
@@ -1,23 +0,0 @@
-module Distribution
- # Calculate cdf and inverse cdf for Chi Square Distribution.
- #
- # Based on Statistics2 module
- #
- module ChiSquare
- class << self
- # Return the P-value of the corresponding integral with
- # k degrees of freedom
- def p_value(pr,k)
- Statistics2.pchi2X_(k.to_i, pr)
- end
- # Chi-square cumulative distribution function (cdf).
- #
- # Returns the integral of Chi-squared distribution
- # with k degrees of freedom over [0, x]
- #
- def cdf(x, k)
- Statistics2.chi2dist(k.to_i,x)
- end
- end
- end
-end
View
35 lib/distribution/f.rb
@@ -1,35 +0,0 @@
-module Distribution
- # Calculate cdf and inverse cdf for Fisher Distribution.
- # Uses Statistics2 module
- module F
- class << self
- # Return the P-value of the corresponding integral with
- # k degrees of freedom
- #
- # Distribution::F.p_value(0.95,1,2)
- def p_value(pr,k1,k2)
- # Statistics2 has some troubles with extreme f values
- if Distribution.has_gsl?
- GSL::Cdf.fdist_Pinv(pr,k1,k2)
- else
- #puts "F:#{k1}, #{k2},#{pr}"
- Statistics2.pfdist(k1,k2, pr)
- end
- end
- # F cumulative distribution function (cdf).
- #
- # Returns the integral of F-distribution
- # with k1 and k2 degrees of freedom
- # over [0, x].
- # Distribution::F.cdf(20,3,2)
- #
- def cdf(x, k1, k2)
- if Distribution.has_gsl?
- GSL::Cdf.fdist_P(x.to_f,k1,k2)
- else
- Statistics2.fdist(k1, k2,x)
- end
- end
- end
- end
-end
View
64 lib/distribution/normal.rb
@@ -1,64 +0,0 @@
-module Distribution
- # Calculate cdf and inverse cdf for Normal Distribution.
- # Uses Statistics2 module
- module Normal
- class << self
- def rng_ugaussian
- rng_gaussian(0,1,nil)
- end
- # Return a proc which return a random number within a gaussian distribution -> N(+mean+,+sigma+)
- # == Reference:
- # * http://www.taygeta.com/random/gaussian.html
- def rng_gaussian(mean=0,sigma=1,seed=nil)
- seed||=rand(1e8)
- if Distribution.has_gsl?
- rng=GSL::Rng.alloc(GSL::Rng::MT19937,seed)
- lambda { mean+rng.gaussian(sigma)}
- else
- returned,y1,y2=0,0,0
- lambda {
- if returned==0
- begin
- x1 = 2.0 * rand - 1.0
- x2 = 2.0 * rand - 1.0
- w = x1 * x1 + x2 * x2
- end while ( w >= 1.0 )
- w = Math::sqrt( (-2.0 * Math::log( w ) ) / w )
- y1 = x1 * w
- y2 = x2 * w
- returned=1
- y1*sigma + mean
- else
- returned=0
- y2 * sigma + mean
- end
- }
- end
- end
- # Return the P-value of the corresponding integral
- def p_value(pr)
- Statistics2.pnormaldist(pr)
- end
- # Normal cumulative distribution function (cdf).
- #
- # Returns the integral of normal distribution
- # over (-Infty, x].
- #
- def cdf(x)
- Statistics2.normaldist(x)
- end
-
- if Distribution.has_gsl?
- alias :cdf_ruby :cdf
- def cdf(x) # :nodoc:
- GSL::Cdf::gaussian_P(x)
- end
- end
- # Normal probability density function (pdf)
- # With x=0 and sigma=1
- def pdf(x)
- (1.0/Math::sqrt(2*Math::PI))*Math::exp(-(x**2/2.0))
- end
- end
- end
-end
View
284 lib/distribution/normalbivariate.rb
@@ -1,284 +0,0 @@
-module Distribution
- # Calculate pdf and cdf for bivariate normal distribution.
- #
- # Pdf if easy to calculate, but CDF is not trivial. Several papers
- # describe methods to calculate the integral.
- #
- # Three methods are implemented on this module:
- # * Genz:: Used by default, with improvement to calculate p on rho > 0.95
- # * Hull:: Port from a C++ code
- # * Jantaravareerat:: Iterative (and slow)
- #
-
- module NormalBivariate
-
- class << self
- SIDE=0.1 # :nodoc:
- LIMIT=5 # :nodoc:
- # Return the partial derivative of cdf over x, with y and rho constant
- # Reference:
- # * Tallis, 1962, p.346, cited by Olsson, 1979
- def partial_derivative_cdf_x(x,y,rho)
- Distribution::Normal.pdf(x) * Distribution::Normal.cdf((y-rho*x).quo( Math::sqrt( 1 - rho**2 )))
- end
- alias :pd_cdf_x :partial_derivative_cdf_x
- # Probability density function for a given x, y and rho value.
- #
- # Source: http://en.wikipedia.org/wiki/Multivariate_normal_distribution
- def pdf(x,y, rho, s1=1.0, s2=1.0)
- 1.quo(2 * Math::PI * s1 * s2 * Math::sqrt( 1 - rho**2 )) * (Math::exp(-(1.quo(2*(1-rho**2))) *
- ((x**2.quo(s1)) + (y**2.quo(s2)) - (2*rho*x*y).quo(s1*s2))))
- end
-
- def f(x,y,aprime,bprime,rho)
- r=aprime*(2*x-aprime)+bprime*(2*y-bprime)+2*rho*(x-aprime)*(y-bprime)
- Math::exp(r)
- end
-
- # CDF for a given x, y and rho value.
- # Uses Genz algorithm (cdf_genz method).
- #
- def cdf(a,b,rho)
- cdf_genz(a,b,rho)
- end
-
- def sgn(x)
- if(x>=0)
- 1
- else
- -1
- end
- end
-
- # Normal cumulative distribution function (cdf) for a given x, y and rho.
- # Based on Hull (1993, cited by Arne, 2003)
- #
- # References:
- # * Arne, B.(2003). Financial Numerical Recipes in C ++. Available on http://finance.bi.no/~bernt/gcc_prog/recipes/recipes/node23.html
- def cdf_hull(a,b,rho)
- #puts "a:#{a} - b:#{b} - rho:#{rho}"
- if (a<=0 and b<=0 and rho<=0)
- # puts "ruta 1"
- aprime=a.quo(Math::sqrt(2.0*(1.0-rho**2)))
- bprime=b.quo(Math::sqrt(2.0*(1.0-rho**2)))
- aa=[0.3253030, 0.4211071, 0.1334425, 0.006374323]
- bb=[0.1337764, 0.6243247, 1.3425378, 2.2626645]
- sum=0
- 4.times do |i|
- 4.times do |j|
- sum+=aa[i]*aa[j] * f(bb[i], bb[j], aprime, bprime,rho)
- end
- end
- sum=sum*(Math::sqrt(1.0-rho**2).quo(Math::PI))
- return sum
- elsif(a*b*rho<=0.0)
-
- #puts "ruta 2"
- if(a<=0 and b>=0 and rho>=0)
- return Distribution::Normal.cdf(a) - cdf(a,-b,-rho)
- elsif (a>=0.0 and b<=0.0 and rho>=0)
- return Distribution::Normal.cdf(b) - cdf(-a,b,-rho)
- elsif (a>=0.0 and b>=0.0 and rho<=0)
- return Distribution::Normal.cdf(a) + Distribution::Normal.cdf(b) - 1.0 + cdf(-a,-b,rho)
- end
- elsif (a*b*rho>=0.0)
- #puts "ruta 3"
- denum=Math::sqrt(a**2 - 2*rho*a*b + b**2)
- rho1=((rho*a-b)*sgn(a)).quo(denum)
- rho2=((rho*b-a)*sgn(b)).quo(denum)
- delta=(1.0-sgn(a)*sgn(b)).quo(4)
- #puts "#{rho1} - #{rho2}"
- return cdf(a, 0.0, rho1) + cdf(b, 0.0, rho2) - delta
- end
- raise "Should'nt be here! #{a} - #{b} #{rho}"
- end
-
- # CDF. Iterative method by Jantaravareerat (n/d)
- #
- # Reference:
- # * Jantaravareerat, M. & Thomopoulos, N. (n/d). Tables for standard bivariate normal distribution
-
- def cdf_jantaravareerat(x,y,rho,s1=1,s2=1)
- # Special cases
- return 1 if x>LIMIT and y>LIMIT
- return 0 if x<-LIMIT or y<-LIMIT
- return Distribution::Normal.cdf(y) if x>LIMIT
- return Distribution::Normal.cdf(x) if y>LIMIT
-
- #puts "x:#{x} - y:#{y}"
- x=-LIMIT if x<-LIMIT
- x=LIMIT if x>LIMIT
- y=-LIMIT if y<-LIMIT
- y=LIMIT if y>LIMIT
-
- x_squares=((LIMIT+x) / SIDE).to_i
- y_squares=((LIMIT+y) / SIDE).to_i
- sum=0
- x_squares.times do |i|
- y_squares.times do |j|
- z1=-LIMIT+(i+1)*SIDE
- z2=-LIMIT+(j+1)*SIDE
- #puts " #{z1}-#{z2}"
- h=(pdf(z1,z2,rho,s1,s2)+pdf(z1-SIDE,z2,rho,s1,s2)+pdf(z1,z2-SIDE,rho,s1,s2) + pdf(z1-SIDE,z2-SIDE,rho,s1,s2)).quo(4)
- sum+= (SIDE**2)*h # area
- end
- end
- sum
- end
- # Normal cumulative distribution function (cdf) for a given x, y and rho.
- # Ported from Fortran code by Alan Genz
- #
- # Original documentation
- # DOUBLE PRECISION FUNCTION BVND( DH, DK, R )
- # A function for computing bivariate normal probabilities.
- #
- # Alan Genz
- # Department of Mathematics
- # Washington State University
- # Pullman, WA 99164-3113
- # Email : alangenz_AT_wsu.edu
- #
- # This function is based on the method described by
- # Drezner, Z and G.O. Wesolowsky, (1989),
- # On the computation of the bivariate normal integral,
- # Journal of Statist. Comput. Simul. 35, pp. 101-107,
- # with major modifications for double precision, and for |R| close to 1.
- #
- # Original location:
- # * http://www.math.wsu.edu/faculty/genz/software/fort77/tvpack.f
- def cdf_genz(x,y,rho)
- dh=-x
- dk=-y
- r=rho
- twopi = 6.283185307179586
-
- w=11.times.collect {[nil]*4};
- x=11.times.collect {[nil]*4}
-
- data=[
- 0.1713244923791705E+00, -0.9324695142031522E+00,
- 0.3607615730481384E+00, -0.6612093864662647E+00,
- 0.4679139345726904E+00, -0.2386191860831970E+00]
-
- (1..3).each {|i|
- w[i][1]=data[(i-1)*2]
- x[i][1]=data[(i-1)*2+1]
-
- }
- data=[
- 0.4717533638651177E-01,-0.9815606342467191E+00,
- 0.1069393259953183E+00,-0.9041172563704750E+00,
- 0.1600783285433464E+00,-0.7699026741943050E+00,
- 0.2031674267230659E+00,-0.5873179542866171E+00,
- 0.2334925365383547E+00,-0.3678314989981802E+00,
- 0.2491470458134029E+00,-0.1252334085114692E+00]
- (1..6).each {|i|
- w[i][2]=data[(i-1)*2]
- x[i][2]=data[(i-1)*2+1]
-
-
- }
- data=[
- 0.1761400713915212E-01,-0.9931285991850949E+00,
- 0.4060142980038694E-01,-0.9639719272779138E+00,
- 0.6267204833410906E-01,-0.9122344282513259E+00,
- 0.8327674157670475E-01,-0.8391169718222188E+00,
- 0.1019301198172404E+00,-0.7463319064601508E+00,
- 0.1181945319615184E+00,-0.6360536807265150E+00,
- 0.1316886384491766E+00,-0.5108670019508271E+00,
- 0.1420961093183821E+00,-0.3737060887154196E+00,
- 0.1491729864726037E+00,-0.2277858511416451E+00,
- 0.1527533871307259E+00,-0.7652652113349733E-01]
-
- (1..10).each {|i|
- w[i][3]=data[(i-1)*2]
- x[i][3]=data[(i-1)*2+1]
-
-
- }
-
-
- if ( r.abs < 0.3 )
- ng = 1
- lg = 3
- elsif ( r.abs < 0.75 )
- ng = 2
- lg = 6
- else
- ng = 3
- lg = 10
- end
-
-
- h = dh
- k = dk
- hk = h*k
- bvn = 0
- if ( r.abs < 0.925 )
- if ( r.abs > 0 )
- hs = ( h*h + k*k ).quo(2)
- asr = Math::asin(r)
- (1..lg).each do |i|
- [-1,1].each do |is|
- sn = Math::sin(asr*(is* x[i][ng]+1).quo(2) )
- bvn = bvn + w[i][ng] * Math::exp( ( sn*hk-hs ).quo( 1-sn*sn ) )
- end # do
- end # do
- bvn = bvn*asr.quo( 2*twopi )
- end # if
- bvn = bvn + Distribution::Normal.cdf(-h) * Distribution::Normal.cdf(-k)
-
-
- else # r.abs
- if ( r < 0 )
- k = -k
- hk = -hk
- end
-
- if ( r.abs < 1 )
- as = ( 1 - r )*( 1 + r )
- a = Math::sqrt(as)
- bs = ( h - k )**2
- c = ( 4 - hk ).quo(8)
- d = ( 12 - hk ).quo(16)
- asr = -( bs.quo(as) + hk ).quo(2)
- if ( asr > -100 )
- bvn = a*Math::exp(asr) * ( 1 - c*( bs - as )*( 1 - d*bs.quo(5) ).quo(3) + c*d*as*as.quo(5) )
- end
- if ( -hk < 100 )
- b = Math::sqrt(bs)
- bvn = bvn - Math::exp( -hk.quo(2) ) * Math::sqrt(twopi)*Distribution::Normal.cdf(-b.quo(a))*b *
- ( 1 - c*bs*( 1 - d*bs.quo(5) ).quo(3) )
- end
-
-
- a = a.quo(2)
- (1..lg).each do |i|
- [-1,1].each do |is|
- xs = (a*( is*x[i][ng] + 1 ) )**2
- rs = Math::sqrt( 1 - xs )
- asr = -( bs/xs + hk ).quo(2)
- if ( asr > -100 )
- bvn = bvn + a*w[i][ng] * Math::exp( asr ) *
- ( Math::exp( -hk*( 1 - rs ).quo(2*( 1 + rs ) ) ) .quo(rs) - ( 1 + c*xs*( 1 + d*xs ) ) )
- end
- end
- end
- bvn = -bvn/twopi
- end
-
- if ( r > 0 )
- bvn = bvn + Distribution::Normal.cdf(-[h,k].max)
- else
- bvn = -bvn
- if ( k > h )
- bvn = bvn + Distribution::Normal.cdf(k) - Distribution::Normal.cdf(h)
- end
- end
- end
- bvn
- end
- private :f, :sgn
- end
- end
-end
View
73 lib/distribution/normalmultivariate.rb
@@ -1,73 +0,0 @@
-module Distribution
- # Calculate cdf and inverse cdf for Multivariate Distribution.
- module NormalMultivariate
- class << self
- # Returns multivariate cdf distribution
- # * a is the array of lower values
- # * b is the array of higher values
- # * s is an symmetric positive definite covariance matrix
- def cdf(aa,bb,sigma, epsilon=0.0001, alpha=2.5, max_iterations=100) # :nodoc:
- raise "Doesn't work yet"
- a=[nil]+aa
- b=[nil]+bb
- m=aa.size
- sigma=sigma.to_gsl if sigma.respond_to? :to_gsl
-
- cc=GSL::Linalg::Cholesky.decomp(sigma)
- c=cc.lower
- intsum=0
- varsum=0
- n=0
- d=Array.new(m+1,nil)
- e=Array.new(m+1,nil)
- f=Array.new(m+1,nil)
- (1..m).each {|i|
- d[i]=0.0 if a[i].nil?
- e[i]=1.0 if b[i].nil?
- }
- d[1]=uPhi(a[1].quo( c[0,0])) unless d[1]==0
- e[1]=uPhi(b[1].quo( c[0,0])) unless e[1]==1
- f[1]=e[1]-d[1]
-
- error=1000
- begin
- w=(m+1).times.collect {|i| rand*epsilon}
- y=[]
- (2..m).each do |i|
- y[i-1]=iPhi(d[i-1] + w[i-1] * (e[i-1] - d[i-1]))
- sumc=0
- (1..(i-1)).each do |j|
- sumc+=c[i-1, j-1]*y[j]
- end
-
- if a[i]!=nil
- d[i]=uPhi((a[i]-sumc).quo(c[i-1,i-1]))
- end
- # puts "sumc:#{sumc}"
-
- if b[i]!=nil
- #puts "e[#{i}] :#{c[i-1,i-1]}"
- e[i]=uPhi((b[i]-sumc).quo(c[i-1, i-1]))
- end
- f[i]=(e[i]-d[i])*f[i-1]
- end
- intsum+=intsum+f[m]
- varsum=varsum+f[m]**2
- n+=1
- error=alpha*Math::sqrt((varsum.quo(n) - (intsum.quo(n))**2).quo(n))
- end while(error>epsilon and n<max_iterations)
-
- f=intsum.quo(n)
- #p intsum
- #puts "f:#{f}, n:#{n}, error:#{error}"
- f
- end
- def iPhi(pr)
- Distribution::Normal.p_value(pr)
- end
- def uPhi(x)
- Distribution::Normal.cdf(x)
- end
- end
- end
-end
View
55 lib/distribution/t.rb
@@ -1,55 +0,0 @@
-require 'rbconfig'
-module Distribution
-
- # Calculate cdf and inverse cdf for T Distribution.
- # Uses Statistics2 Module.
- module T
- class << self
- # Return the P-value of the corresponding integral with
- # k degrees of freedom
- def p_value(pr,k)
- Statistics2.ptdist(k, pr)
- end
- # T cumulative distribution function (cdf).
- #
- # Returns the integral of t-distribution
- # with n degrees of freedom over (-Infty, x].
- #
- def cdf(x,k)
- if RbConfig::CONFIG['arch']=~/i686/
- tdist(k, x)
- else
- Statistics2.tdist(k,x)
- end
- end
-
- # Returns the integral of t-distribution with n degrees of freedom over (-Infty, x].
- def tdist(n, t)
- p_t(n, t)
- end
-
- # t-distribution ([1])
- # (-\infty, x]
- def p_t(df, t)
- c2 = df.to_f / (df + t * t);
- s = Math.sqrt(1.0 - c2)
- s = -s if t < 0.0
- p = 0.0;
- i = df % 2 + 2
- while i <= df
- p += s
- s *= (i - 1) * c2 / i
- i += 2
- end
- if df.is_a? Float or df & 1 != 0
- 0.5+(p*Math.sqrt(c2)+Math.atan(t/Math.sqrt(df)))/Math::PI
- else
- (1.0 + p) / 2.0
- end
- end
-
-
-
- end
- end
-end
View
2 lib/statsample/factor/parallelanalysis.rb
@@ -126,7 +126,7 @@ def compute
@ds_eigenvalues=Statsample::Dataset.new((1..@n_variables).map{|v| "ev_%05d" % v})
@ds_eigenvalues.fields.each {|f| @ds_eigenvalues[f].type=:scale}
if bootstrap_method==:parameter or bootstrap_method==:random
- rng = Distribution::Normal.rng_ugaussian
+ rng = Distribution::Normal.rng
end
@iterations.times do |i|
View
2 lib/statsample/shorthand.rb
@@ -24,7 +24,7 @@ def c(*args)
end
# Random generation for the normal distribution
def rnorm(n,mean=0,sd=1)
- rng=Distribution::Normal.rng_gaussian(mean,sd)
+ rng=Distribution::Normal.rng(mean,sd)
Statsample::Vector.new_scale(n) { rng.call}
end
# Creates a new Statsample::Dataset
View
2 lib/statsample/test/chisquare.rb
@@ -26,7 +26,7 @@ def chi_square
@value
end
def probability
- 1-Distribution::ChiSquare.cdf(@value,@df)
+ 1-Distribution::ChiSquare.cdf(@value.to_f,@df)
end
def compute_chi
sum=0
View
72 test/test_distribution.rb
@@ -1,72 +0,0 @@
-require(File.expand_path(File.dirname(__FILE__)+'/helpers_tests.rb'))
-require 'distribution'
-
-
-class DistributionTestCase < MiniTest::Unit::TestCase
- def test_chi
- if Distribution.has_gsl?
- [2,3,4,5].each{|k|
- chis=rand()*10
- area=Distribution::ChiSquare.cdf(chis, k)
- assert_in_delta(area, GSL::Cdf.chisq_P(chis,k),0.0001)
- assert_in_delta(chis, Distribution::ChiSquare.p_value(area,k),0.0001,"Error on prob #{area} and k #{k}")
- }
- end
- end
- def test_t
- if Distribution.has_gsl?
- [-2,0.1,0.5,1,2].each{|t|
- [2,5,10].each{|n|
- area=Distribution::T.cdf(t,n)
- assert_in_delta(area, GSL::Cdf.tdist_P(t,n),0.0001)
- assert_in_delta(Distribution::T.p_value(area,n), GSL::Cdf.tdist_Pinv(area,n),0.0001)
-
- }
- }
- end
- end
- def test_normal
- if Distribution.has_gsl?
- [-2,0.1,0.5,1,2].each{|x|
- area=Distribution::Normal.cdf(x)
- assert_in_delta(area, GSL::Cdf.ugaussian_P(x),0.0001)
- assert_in_delta(Distribution::Normal.p_value(area), GSL::Cdf.ugaussian_Pinv(area),0.0001)
- assert_in_delta(Distribution::Normal.pdf(x), GSL::Ran::ugaussian_pdf(x),0.0001)
- }
- end
- end
- def test_normal_bivariate
- if Distribution.has_gsl?
- [0.2,0.4,0.6,0.8,0.9, 0.99,0.999,0.999999].each {|rho|
- assert_in_delta(GSL::Ran::bivariate_gaussian_pdf(0, 0, 1,1,rho), Distribution::NormalBivariate.pdf(0,0, rho , 1,1),1e-8)
-
- }
- end
-
- [-3,-2,-1,0,1,1.5].each {|x|
- assert_in_delta(Distribution::NormalBivariate.cdf_hull(x,x,0.5), Distribution::NormalBivariate.cdf_genz(x,x,0.5), 0.001)
- #assert_in_delta(Distribution::NormalBivariate.cdf_genz(x,x,0.5), Distribution::NormalBivariate.cdf_jantaravareerat(x,x,0.5), 0.001)
- }
-
- assert_in_delta(0.686, Distribution::NormalBivariate.cdf(2,0.5,0.5), 0.001)
- assert_in_delta(0.498, Distribution::NormalBivariate.cdf(2,0.0,0.5), 0.001)
- assert_in_delta(0.671, Distribution::NormalBivariate.cdf(1.5,0.5,0.5), 0.001)
-
- assert_in_delta(Distribution::Normal.cdf(0), Distribution::NormalBivariate.cdf(10,0,0.9), 0.001)
- end
- def test_f
- if Distribution.has_gsl?
- [0.1,0.5,1,2,10,20,30].each{|f|
- [2,5,10].each{|n2|
- [2,5,10].each{|n1|
- area=Distribution::F.cdf(f,n1,n2)
- assert_in_delta(area, GSL::Cdf.fdist_P(f,n1,n2),0.0001)
- assert_in_delta(Distribution::F.p_value(area,n1,n2), GSL::Cdf.fdist_Pinv(area,n1,n2),0.0001)
-
- }
- }
- }
- end
- end
-
-end
View
4 test/test_factor.rb
@@ -35,7 +35,7 @@ def test_covariance_matrix
end
def test_principalcomponents_ruby_gsl
- ran=Distribution::Normal.rng_ugaussian
+ ran=Distribution::Normal.rng
# @r=::Rserve::Connection.new
@@ -83,7 +83,7 @@ def test_principalcomponents()
end
def principalcomponents(gsl)
- ran=Distribution::Normal.rng_ugaussian
+ ran=Distribution::Normal.rng
samples=50
x1=samples.times.map { ran.call()}.to_scale
x2=samples.times.map {|i| ran.call()*0.5+x1[i]*0.5}.to_scale
View
2 test/test_factor_pa.rb
@@ -13,7 +13,7 @@ def test_parallelanalysis_with_data
samples=100
variables=10
iterations=50
- rng = Distribution::Normal.rng_ugaussian
+ rng = Distribution::Normal.rng
f1=samples.times.collect {rng.call}.to_scale
f2=samples.times.collect {rng.call}.to_scale
vectors={}
View
2 test/test_logit.rb
@@ -26,7 +26,7 @@ class StatsampleLogitTestCase < MiniTest::Unit::TestCase
should "return same similat values to as R gml" do
r=Rserve::Connection.new
- ran=Distribution::Normal.rng_ugaussian
+ ran=Distribution::Normal.rng
samples=100
a,b,c=ran.call,ran.call,ran.call
logit=lambda {|x| Math.exp(x) / (1+Math.exp(x))}

0 comments on commit b2c233f

Please sign in to comment.