Verify distribution uniformity/Chi-squared test: Difference between revisions

Content added Content deleted

Inline

@@ Line 410: / Line 410: @@
 {{libheader|GNU Scientific Library}}
-Instead of implementing the incomplete gamma function by ourselves, we bind to GNU Scientific Library; so we need a module to interface to the function we need (<tt>gsl_sf_gamma_inc</tt>)
+Instead of implementing the chi-squared distribution by ourselves, we bind to GNU Scientific Library; so we need a module to interface to the function we need (<tt>gsl_cdf_chisq_Q</tt>)
-<lang fortran>module GSLMiniBind
+<lang fortran>module gsl_mini_bind_m
-  implicit none
+    use iso_c_binding
-  interface
+    implicit none
-     real(c_double) function gsl_sf_gamma_inc(a, x) bind(C)
+    private
-       use iso_c_binding
-       real(c_double), value, intent(in) :: a, x
+    public :: p_value
-     end function gsl_sf_gamma_inc
-  end interface
+    interface
-end module GSLMiniBind</lang>
+        function gsl_cdf_chisq_q(x, nu) bind(c, name='gsl_cdf_chisq_Q')
+            import
+            real(c_double), value :: x
+            real(c_double), value :: nu
+            real(c_double) :: gsl_cdf_chisq_q
+        end function gsl_cdf_chisq_q
+    end interface
+contains
+    !> Get p-value from chi-square distribution
+    real function p_value(x, df)
+        real, intent(in) :: x
+        integer, intent(in) :: df
+        p_value = real(gsl_cdf_chisq_q(real(x, c_double), real(df, c_double)))
+    end function p_value
+end module gsl_mini_bind_m</lang>
 Now we're ready to complete the task.
-<lang fortran>program ChiTest
+<lang fortran>program chi2test
-  use GSLMiniBind
-  use iso_c_binding
-  implicit none
+    use gsl_mini_bind_m, only: p_value
-  real, dimension(5) :: dset1 = (/ 199809., 200665., 199607., 200270., 199649. /)
+    implicit none
-  real, dimension(5) :: dset2 = (/ 522573., 244456., 139979.,  71531.,  21461. /)
+    real :: dset1(5) = [199809., 200665., 199607., 200270., 199649.]
-  real :: dist, prob
+    real :: dset2(5) = [522573., 244456., 139979., 71531., 21461.]
-  integer :: dof
+    real :: dist, prob
-  print *, "Dataset 1:"
+    integer :: dof
-  print *, dset1
-  dist = chi2UniformDistance(dset1)
-  dof = size(dset1) - 1
-  write(*, '(A,I4,A,F12.4)') 'dof: ', dof, '   distance: ', dist
-  prob = chi2Probability(dof, dist)
-  write(*, '(A,F9.6)') 'probability: ', prob
-  write(*, '(A,L)') 'uniform? ', chiIsUniform(dset1, 0.05)
+    write (*, '(A)', advance='no') "Dataset 1:"
-  ! Lazy copy/past :|
+    write (*, '(5(F12.4,:,1x))') dset1
-  print *, "Dataset 2:"
-  print *, dset2
-  dist = chi2UniformDistance(dset2)
-  dof = size(dset2) - 1
-  write(*, '(A,I4,A,F12.4)') 'dof: ', dof, '   distance: ', dist
-  prob = chi2Probability(dof, dist)
-  write(*, '(A,F9.6)') 'probability: ', prob
-  write(*, '(A,L)') 'uniform? ', chiIsUniform(dset2, 0.05)
+    dist = chisq(dset1)
-contains
+    dof = size(dset1) - 1
+    write (*, '(A,I4,A,F12.4)') 'dof: ', dof, '   chisq: ', dist
+    prob = p_value(dist, dof)
+    write (*, '(A,F12.4)') 'probability: ', prob
+    write (*, '(A,L)') 'uniform? ', prob > 0.05
+    ! Lazy copy/past :|
-function chi2Probability(dof, distance)
+    write (*, '(/A)', advance='no') "Dataset 2:"
-  real :: chi2Probability
+    write (*, '(5(F12.4,:,1x))') dset2
-  integer, intent(in) :: dof
-  real, intent(in) :: distance
+    dist = chisq(dset2)
-  ! in order to make this work, we need linking with GSL library
+    dof = size(dset2) - 1
-  chi2Probability = gsl_sf_gamma_inc(real(0.5*dof, c_double), real(0.5*distance, c_double))
+    write (*, '(A,I4,A,F12.4)') 'dof: ', dof, '   chisq: ', dist
-end function chi2Probability
+    prob = p_value(dist, dof)
+    write (*, '(A,F12.4)') 'probability: ', prob
+    write (*, '(A,L)') 'uniform? ', prob > 0.05
+contains
-function chiIsUniform(dset, significance)
-  logical :: chiIsUniform
-  real, dimension(:), intent(in) :: dset
-  real, intent(in) :: significance
+    !> Get chi-square value for a set of data `ds`
-  integer :: dof
-  real :: dist
+    real function chisq(ds)
+        real, intent(in) :: ds(:)
+        real :: expected, summa
-  dof = size(dset) - 1
-  dist = chi2UniformDistance(dset)
-  chiIsUniform = chi2Probability(dof, dist) > significance
-end function chiIsUniform
-function chi2UniformDistance(ds)
-  real :: chi2UniformDistance
-  real, dimension(:), intent(in) :: ds
+        expected = sum(ds)/size(ds)
-  real :: expected, summa = 0.0
+        summa = sum((ds - expected)**2)
+        chisq = summa/expected
+    end function chisq
-  expected = sum(ds) / size(ds)
-  summa = sum( (ds - expected) ** 2 )
-  chi2UniformDistance = summa / expected
-end function chi2UniformDistance
-end program ChiTest</lang>
+end program chi2test</lang>
 =={{header|Go}}==