Changeset 3743


Ignore:
Timestamp:
Jul 12, 2018, 2:43:25 AM (5 years 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.