yat/regression/detail/Cox.h

Code
Comments
Other
Rev Date Author Line
4198 19 Aug 22 peter 1 #ifndef theplu_yat_regression_detail_cox
4198 19 Aug 22 peter 2 #define theplu_yat_regression_detail_cox
4198 19 Aug 22 peter 3
4198 19 Aug 22 peter 4 // $Id$
4198 19 Aug 22 peter 5
4198 19 Aug 22 peter 6 /*
4198 19 Aug 22 peter 7   Copyright (C) 2022 Peter Johansson
4198 19 Aug 22 peter 8
4198 19 Aug 22 peter 9   This file is part of the yat library, https://dev.thep.lu.se/yat
4198 19 Aug 22 peter 10
4198 19 Aug 22 peter 11   The yat library is free software; you can redistribute it and/or
4198 19 Aug 22 peter 12   modify it under the terms of the GNU General Public License as
4198 19 Aug 22 peter 13   published by the Free Software Foundation; either version 3 of the
4198 19 Aug 22 peter 14   License, or (at your option) any later version.
4198 19 Aug 22 peter 15
4198 19 Aug 22 peter 16   The yat library is distributed in the hope that it will be useful,
4198 19 Aug 22 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
4198 19 Aug 22 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4198 19 Aug 22 peter 19   General Public License for more details.
4198 19 Aug 22 peter 20
4198 19 Aug 22 peter 21   You should have received a copy of the GNU General Public License
4198 19 Aug 22 peter 22   along with yat. If not, see <https://www.gnu.org/licenses/>.
4198 19 Aug 22 peter 23 */
4198 19 Aug 22 peter 24
4198 19 Aug 22 peter 25 #include <yat/utility/VectorBase.h>
4198 19 Aug 22 peter 26
4198 19 Aug 22 peter 27 #include <gsl/gsl_cdf.h>
4198 19 Aug 22 peter 28
4198 19 Aug 22 peter 29 #include <boost/math/tools/roots.hpp>
4198 19 Aug 22 peter 30
4198 19 Aug 22 peter 31 #include <algorithm>
4198 19 Aug 22 peter 32 #include <memory>
4198 19 Aug 22 peter 33
4198 19 Aug 22 peter 34 namespace theplu {
4198 19 Aug 22 peter 35 namespace yat {
4198 19 Aug 22 peter 36 namespace regression {
4198 19 Aug 22 peter 37 namespace cox {
4198 19 Aug 22 peter 38
4198 19 Aug 22 peter 39   template<typename T>
4198 19 Aug 22 peter 40   struct traits
4198 19 Aug 22 peter 41   {
4198 19 Aug 22 peter 42     typedef const T& const_reference;
4198 19 Aug 22 peter 43     //typedef T& reference;
4198 19 Aug 22 peter 44   };
4198 19 Aug 22 peter 45
4198 19 Aug 22 peter 46
4198 19 Aug 22 peter 47   template<>
4198 19 Aug 22 peter 48   struct traits<utility::Vector>
4198 19 Aug 22 peter 49   {
4198 19 Aug 22 peter 50     typedef const utility::VectorBase& const_reference;
4198 19 Aug 22 peter 51     //typedef utility::VectorMutable& reference;
4198 19 Aug 22 peter 52   };
4198 19 Aug 22 peter 53
4198 19 Aug 22 peter 54
4198 19 Aug 22 peter 55   template<typename T>
4198 19 Aug 22 peter 56   class Implementation
4198 19 Aug 22 peter 57   {
4198 19 Aug 22 peter 58   public:
4198 19 Aug 22 peter 59     void add(typename traits<T>::const_reference x, double time, bool event)
4198 19 Aug 22 peter 60     {
4198 19 Aug 22 peter 61       data_.resize(data_.size()+1);
4198 19 Aug 22 peter 62       data_.back().time = time;
4198 19 Aug 22 peter 63       data_.back().x = x;
4198 19 Aug 22 peter 64       data_.back().event = event;
4198 19 Aug 22 peter 65     }
4198 19 Aug 22 peter 66
4198 19 Aug 22 peter 67     void clear(void)
4198 19 Aug 22 peter 68     {
4198 19 Aug 22 peter 69       data_.clear();
4198 19 Aug 22 peter 70     }
4198 19 Aug 22 peter 71
4198 19 Aug 22 peter 72   protected:
4198 19 Aug 22 peter 73     struct Data
4198 19 Aug 22 peter 74     {
4198 19 Aug 22 peter 75       double time;
4198 19 Aug 22 peter 76       T x;
4198 19 Aug 22 peter 77       bool event; // false if censored
4198 19 Aug 22 peter 78       double theta(typename traits<T>::const_reference beta) const
4198 19 Aug 22 peter 79       { return std::exp(beta * x); }
4198 19 Aug 22 peter 80     };
4198 19 Aug 22 peter 81     std::vector<Data> data_;
4198 19 Aug 22 peter 82
4198 19 Aug 22 peter 83     struct TimePoint
4198 19 Aug 22 peter 84     {
4198 19 Aug 22 peter 85       typedef typename std::vector<Data>::const_iterator iterator;
4198 19 Aug 22 peter 86       iterator events_begin(void) const { return begin; }
4198 19 Aug 22 peter 87       iterator events_end(void) const { return middle; }
4198 19 Aug 22 peter 88       iterator censored_begin(void) const { return middle; }
4198 19 Aug 22 peter 89       iterator censored_end(void) const { return end; }
4198 19 Aug 22 peter 90       double time(void) const {return begin->time; }
4198 19 Aug 22 peter 91       // number of events
4198 19 Aug 22 peter 92       size_t size(void) const
4198 19 Aug 22 peter 93       { return events_end() - events_begin(); }
4198 19 Aug 22 peter 94
4198 19 Aug 22 peter 95       iterator begin;
4198 19 Aug 22 peter 96       iterator middle;
4198 19 Aug 22 peter 97       iterator end;
4198 19 Aug 22 peter 98     };
4198 19 Aug 22 peter 99     std::vector<TimePoint> times_;
4198 19 Aug 22 peter 100
4198 19 Aug 22 peter 101
4198 19 Aug 22 peter 102     void prepare_times(void)
4198 19 Aug 22 peter 103     {
4198 19 Aug 22 peter 104       times_.clear();
4198 19 Aug 22 peter 105       std::sort(data_.begin(), data_.end(),
4198 19 Aug 22 peter 106                 [](const Data& lhs, const Data& rhs) {
4198 19 Aug 22 peter 107                   if (lhs.time != rhs.time)
4198 19 Aug 22 peter 108                     return lhs.time < rhs.time;
4198 19 Aug 22 peter 109                   // sort events before non-events
4198 19 Aug 22 peter 110                   return lhs.event && !rhs.event;
4198 19 Aug 22 peter 111                 });
4198 19 Aug 22 peter 112
4198 19 Aug 22 peter 113       for (auto it = data_.begin(); it!=data_.end(); ) {
4198 19 Aug 22 peter 114         TimePoint time;
4198 19 Aug 22 peter 115         time.begin = it;
4198 19 Aug 22 peter 116         while (it!=data_.end() && it->event && it->time == time.begin->time)
4198 19 Aug 22 peter 117           ++it;
4198 19 Aug 22 peter 118         time.middle = it;
4198 19 Aug 22 peter 119         while (it!=data_.end() && !it->event &&
4198 19 Aug 22 peter 120                it->time == time.begin->time)
4198 19 Aug 22 peter 121           ++it;
4198 19 Aug 22 peter 122         time.end = it;
4198 19 Aug 22 peter 123         times_.push_back(std::move(time));
4198 19 Aug 22 peter 124       }
4198 19 Aug 22 peter 125     }
4198 19 Aug 22 peter 126
4198 19 Aug 22 peter 127   };
4198 19 Aug 22 peter 128 }}}}
4198 19 Aug 22 peter 129 #endif

Generated on Sat Oct 14 02:23:04 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)