Changeset 3743


Ignore:
Timestamp:
Jul 12, 2018, 2:43:25 AM (3 months ago)
Author:
peter
Message:

merge release 0.15.1 into trunk

Location:
trunk
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • trunk

  • trunk/NEWS

    r3713 r3743  
    99
    1010yat 0.15.x series from http://dev.thep.lu.se/yat/svn/branches/0.15-stable
     11
     12version 0.15.1 (released 12 July 2018)
     13  - Fixed bug in calculation of exact p_value when one element was
     14    zero (bug #908)
     15  - Fixed usage of GSL_error, so 'status' is propagated correctly (see r3737)
     16
     17  A complete list of closed tickets can be found here [[br]]
     18  http://dev.thep.lu.se/yat/query?status=closed&milestone=yat+0.15.1
    1119
    1220version 0.15 (released 9 November 2017)
     
    461469Copyright (C) 2010, 2011 Peter Johansson
    462470Copyright (C) 2012 Jari Häkkinen, Peter Johansson
    463 Copyright (C) 2013, 2014, 2015, 2016, 2017 Peter Johansson
     471Copyright (C) 2013, 2014, 2015, 2016, 2017, 2018 Peter Johansson
    464472
    465473This file is part of yat library, http://dev.thep.lu.se/yat
  • trunk/m4/version.m4

    r3713 r3743  
    22#
    33# Copyright (C) 2008, 2009 Jari Häkkinen, Peter Johansson
    4 # Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2017 Peter Johansson
     4# Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2017, 2018 Peter Johansson
    55#
    66# This file is part of the yat library, http://dev.thep.lu.se/yat
     
    8787# yat-0.14.4 11:4:0
    8888# yat-0.15   12:0:0
     89# yat-0.15.1 12:1:0
    8990#
    9091# *Accidently, the libtool number was not updated for yat 0.5
  • trunk/test/fisher.cc

    r3455 r3743  
    22
    33/*
    4         Copyright (C) 2008, 2009, 2010, 2012, 2013, 2014, 2015 Peter Johansson
     4        Copyright (C) 2008, 2009, 2010, 2012, 2013, 2014, 2015, 2018 Peter Johansson
    55
    66        This file is part of the yat library, http://dev.thep.lu.se/yat
     
    2626#include "yat/statistics/Fisher.h"
    2727
     28#include <gsl/gsl_cdf.h>
     29#include <gsl/gsl_randist.h>
    2830#include <climits>
    2931
     
    3436void test_large_numbers(test::Suite&);
    3537void test_yates(test::Suite&);
     38void test_cdf_hypergeometric_Q(test::Suite&);
     39void test_ticket908(test::Suite&);
    3640
    3741int main(int argc, char* argv[])
     
    7377        test_large_numbers(suite);
    7478        test_yates(suite);
     79        test_cdf_hypergeometric_Q(suite);
     80        test_ticket908(suite);
    7581        return suite.return_value();
    7682}
     
    168174        }
    169175}
     176
     177
     178void test_cdf_hypergeometric_Q(test::Suite& suite)
     179{
     180        suite.out() << "\ntest " << __func__ << "\n";
     181        unsigned int n1 = 16;
     182        unsigned int n2 = 3;
     183        unsigned int t = 16;
     184        for (unsigned int k=0; k<=t; ++k) {
     185                // P(X=k)
     186                double p = gsl_ran_hypergeometric_pdf(k, n1, n2, t);
     187                // P(X<=k)
     188                double P = gsl_cdf_hypergeometric_P(k, n1, n2, t);
     189                // P(X>k)
     190                double Q = gsl_cdf_hypergeometric_Q(k, n1, n2, t);
     191                suite.out() << k << " " << p << " " << P << " " << Q << "\n";
     192                double sum = 0;
     193                for (unsigned int x=0; x<=k; ++x)
     194                        sum += gsl_ran_hypergeometric_pdf(x, n1, n2, t);
     195                if (!suite.equal(sum, P, 20)) {
     196                        suite.add(false);
     197                        suite.err() << "error: gsl_cdf_hypergeometric_P for k=" << k << "\n";
     198                }
     199                if (!suite.equal(Q, 1-P, 1000)) {
     200                        suite.add(false);
     201                        suite.err() << "failed: Q = 1-P for k=" << k << "\n";
     202                }
     203        }
     204}
     205
     206
     207// test bug #908 reported in http://dev.thep.lu.se/yat/ticket/908
     208void test_ticket908(test::Suite& suite)
     209{
     210        suite.out() << "===\ntesting ticket 908\n";
     211        statistics::Fisher f;
     212        int a1 = 13;
     213        int a2 = 3;
     214        int b1 = 3;
     215        int b2 = 0;
     216
     217        double correct_p = 0;
     218        unsigned int n1 = a1 + a2;
     219        unsigned int n2 = b1 + b2;
     220        unsigned int t = a1+b1;
     221        for (unsigned int k=0; k<=n1; ++k) {
     222                double P = gsl_ran_hypergeometric_pdf(k, n1, n2, t);
     223                if (P <=        gsl_ran_hypergeometric_pdf(a1, n1, n2, t)) {
     224                        suite.out() << "counting " << k << " " << P << "\n";
     225                        correct_p += P;
     226                }
     227        }
     228        statistics::Fisher fisher;
     229        std::vector<double> p;
     230        fisher.oddsratio(a1, a2, b1, b2);
     231        p.push_back(fisher.p_value());
     232        fisher.oddsratio(a2, a1, b2, b1);
     233        p.push_back(fisher.p_value());
     234        fisher.oddsratio(b1, b2, a1, a2);
     235        p.push_back(fisher.p_value());
     236        fisher.oddsratio(b2, b1, a2, a1);
     237        p.push_back(fisher.p_value());
     238        for (size_t i=0; i<p.size(); ++i) {
     239                if (!suite.add(suite.equal(correct_p, p[i], 20)))
     240                        suite.err() << "error: incorrect p: " << p[i] << " for case "
     241                                                                        << i << "; expected: " << correct_p << "\n";
     242        }
     243}
  • trunk/yat/statistics/Fisher.cc

    r3455 r3743  
    55        Copyright (C) 2005 Peter Johansson
    66        Copyright (C) 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
    7         Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015 Peter Johansson
     7        Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2018 Peter Johansson
    88
    99        This file is part of the yat library, http://dev.thep.lu.se/yat
     
    8787
    8888        /*
     89         (n)       n!
     90         ( )  = --------
     91         (k)    k!(n-k)!
     92
     93
     94
    8995                ( n )
    90                 (k+1)   (k+1)!(n-k)!
    91     ----- = ------------
    92                  (n)    k!(n-k-1)!
     96                (k+1)   k!(n-k)!
     97    ----- = ------------ = (n-k) / (k+1)
     98                 (n)    (k+1)!(n-k-1)!
    9399     (k)
    94100         */
    95101        double Fisher::choose_ratio(unsigned int n, unsigned int k) const
    96102        {
    97                 if (k==0)
    98                         return n;
    99                 if (k==n-1)
    100                         return k+1;
     103                assert(k+1 <= n);
    101104                return static_cast<double>(n-k) / (k+1);
    102105        }
     
    106109                                                                                                                                                        unsigned int n2, unsigned int t) const
    107110        {
     111                // P(X = k+1) / P(X = k) =
     112                // (choose(n1, k+1) * choose(n2, t-k-1) / choose(n1+n2, t) ) /
     113                // (choose(n1, k) * choose(n2, t-k) / choose(n1+n2, t) ) =
     114                // choose(n1, k+1) / choose(n1, k) /( choose(n2, t-k)/choose(n2,t-k-1) ) =
     115                // choose(n1,k+1)/choose(n1,k) / (choose(n2,(t-k-1)+1)/choose(n2,t-k-1))
     116                // choose_ratio(n1,k) / choose_ratio(n2, t-k-1)
     117                // where choose_ratio(a, b) = choose(a, b+1) / choose(a,b)
    108118                return choose_ratio(n1, k) / choose_ratio(n2, t-k-1);
    109119        }
     
    240250        {
    241251                // The hypergeometric distribution is unimodal (only one peak) and the
    242                 // peak (mode) occurs at: (t+1)*(n1+1)/(n+1)
    243 
    244                 double mode = static_cast<double>(a_+c_+1)*static_cast<double>(a_+b_+1)/
    245                         static_cast<double>(a_+b_+c_+d_+1);
     252                // peak (mode) occurs at: floor((t+1)*(n1+1)/(n+1))
     253
     254                double mode = std::floor((a_+c_+1.0)*(a_+b_+1.0)/(a_+b_+c_+d_+1.0));
    246255                if (a_ >= mode)
    247256                        return p_value_exact(a_, a_+b_, c_+d_, a_+c_);
     
    267276                assert(n1);
    268277                assert(n2);
     278
     279                // special case when k=0 and mode (peak of distribution) is at
     280                // zero as well. If mode is not zero 2x2 table should have been
     281                // mirrored before calling this function.
     282                if (k == 0) {
     283                        assert(std::floor((a_+c_+1.0)*(a_+b_+1.0)/(a_+b_+c_+d_+1.0)) == 0.0);
     284                        return 1.0;
     285                }
    269286
    270287                assert(k);
  • trunk/yat/utility/Matrix.cc

    r3691 r3743  
    66        Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
    77        Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
    8         Copyright (C) 2009, 2010, 2011, 2012, 2017 Peter Johansson
     8        Copyright (C) 2009, 2010, 2011, 2012, 2017, 2018 Peter Johansson
    99
    1010        This file is part of the yat library, http://dev.thep.lu.se/yat
     
    210210                int status=gsl_matrix_div_elements(m_, other.gsl_matrix_p());
    211211                if (status)
    212                         throw utility::GSL_error(std::string("matrix::div_elements",status));
     212                        throw utility::GSL_error("matrix::div_elements",status);
    213213        }
    214214
     
    280280                int status=gsl_matrix_mul_elements(m_, other.gsl_matrix_p());
    281281                if (status)
    282                         throw utility::GSL_error(std::string("Matrix::mul_elements",status));
     282                        throw utility::GSL_error("Matrix::mul_elements",status);
    283283        }
    284284
     
    323323                int status=gsl_matrix_swap_columns(m_, i, j);
    324324                if (status)
    325                         throw utility::GSL_error(std::string("Matrix::swap_columns",status));
     325                        throw utility::GSL_error("Matrix::swap_columns",status);
    326326        }
    327327
     
    332332                int status=gsl_matrix_swap_rowcol(m_, i, j);
    333333                if (status)
    334                         throw utility::GSL_error(std::string("Matrix::swap_rowcol",status));
     334                        throw utility::GSL_error("Matrix::swap_rowcol",status);
    335335        }
    336336
     
    341341                int status=gsl_matrix_swap_rows(m_, i, j);
    342342                if (status)
    343                         throw utility::GSL_error(std::string("Matrix::swap_rows",status));
     343                        throw utility::GSL_error("Matrix::swap_rows",status);
    344344        }
    345345
     
    419419                int status=gsl_matrix_add(m_, other.m_);
    420420                if (status)
    421                         throw utility::GSL_error(std::string("Matrix::operator+=", status));
     421                        throw utility::GSL_error("Matrix::operator+=", status);
    422422                return *this;
    423423        }
     
    437437                int status=gsl_matrix_sub(m_, other.m_);
    438438                if (status)
    439                         throw utility::GSL_error(std::string("Matrix::operator-=", status));
     439                        throw utility::GSL_error("Matrix::operator-=", status);
    440440                return *this;
    441441        }
     
    532532                int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p());
    533533                if (status)
    534                         throw utility::GSL_error(std::string("swap(Matrix&,Matrix&)",status));
     534                        throw utility::GSL_error("swap(Matrix&,Matrix&)",status);
    535535        }
    536536
  • trunk/yat/utility/SVD.cc

    r3736 r3743  
    88        Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
    99        Copyright (C) 2009 Jari Häkkinen
    10         Copyright (C) 2011, 2012 Peter Johansson
     10        Copyright (C) 2011, 2012, 2018 Peter Johansson
    1111
    1212        This file is part of the yat library, http://dev.thep.lu.se/yat
     
    8282                }
    8383                if (status)
    84                         throw GSL_error(std::string("SVD::decompose",status));
     84                        throw GSL_error("SVD::decompose",status);
    8585        }
    8686
     
    123123                                                                                                                                         x.gsl_vector_p());
    124124                if (status)
    125                         throw GSL_error(std::string("SVD::solve",status));
     125                        throw GSL_error("SVD::solve",status);
    126126        }
    127127
  • trunk/yat/utility/Vector.cc

    r3691 r3743  
    66        Copyright (C) 2005, 2006, 2007 Jari Häkkinen, Peter Johansson, Markus Ringnér
    77        Copyright (C) 2008 Jari Häkkinen, Peter Johansson
    8         Copyright (C) 2009, 2010, 2012, 2014, 2017 Peter Johansson
     8        Copyright (C) 2009, 2010, 2012, 2014, 2017, 2018 Peter Johansson
    99
    1010        This file is part of the yat library, http://dev.thep.lu.se/yat
     
    211211                int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
    212212                if (status)
    213                         throw utility::GSL_error(std::string("swap(Vector&,Vector&)",status));
     213                        throw utility::GSL_error("swap(Vector&,Vector&)",status);
    214214        }
    215215
  • trunk/yat/utility/VectorBase.cc

    r3661 r3743  
    66        Copyright (C) 2005, 2006, 2007 Jari Häkkinen, Peter Johansson, Markus Ringnér
    77        Copyright (C) 2008 Jari Häkkinen, Peter Johansson
    8         Copyright (C) 2009, 2011, 2012, 2016, 2017 Peter Johansson
     8        Copyright (C) 2009, 2011, 2012, 2016, 2017, 2018 Peter Johansson
    99
    1010        This file is part of the yat library, http://dev.thep.lu.se/trac/yat
     
    179179                if (status) {
    180180                        gsl_permutation_free(p);
    181                         throw utility::GSL_error(std::string("sort_index(vector&,const VectorBase&)",status));
     181                        throw utility::GSL_error("sort_index(vector&,const VectorBase&)",status);
    182182                }
    183183                std::vector<size_t> tmp(p->data,p->data+p->size);
  • trunk/yat/utility/VectorMutable.cc

    r3623 r3743  
    66        Copyright (C) 2005, 2006, 2007 Jari Häkkinen, Peter Johansson, Markus Ringnér
    77        Copyright (C) 2008, 2009 Jari Häkkinen, Peter Johansson
    8         Copyright (C) 2010, 2012, 2016, 2017 Peter Johansson
     8        Copyright (C) 2010, 2012, 2016, 2017, 2018 Peter Johansson
    99
    1010        This file is part of the yat library, http://dev.thep.lu.se/trac/yat
     
    9292                int status=gsl_vector_div(vec_,other.gsl_vector_p());
    9393                if (status)
    94                         throw utility::GSL_error(std::string("VectorMutable::div",status));
     94                        throw utility::GSL_error("VectorMutable::div",status);
    9595        }
    9696
     
    121121                int status=gsl_vector_mul(vec_,other.gsl_vector_p());
    122122                if (status)
    123                         throw utility::GSL_error(std::string("VectorMutable::div",status));
     123                        throw utility::GSL_error("VectorMutable::div",status);
    124124        }
    125125
     
    148148                int status=gsl_vector_add(vec_, other.gsl_vector_p());
    149149                if (status)
    150                         throw utility::GSL_error(std::string("VectorMutable::add", status));
     150                        throw utility::GSL_error("VectorMutable::add", status);
    151151                return *this;
    152152        }
     
    173173                int status=gsl_vector_sub(vec_, other.gsl_vector_p());
    174174                if (status)
    175                         throw utility::GSL_error(std::string("VectorMutable::sub", status));
     175                        throw utility::GSL_error("VectorMutable::sub", status);
    176176                return *this;
    177177        }
Note: See TracChangeset for help on using the changeset viewer.