berry.hpp

Go to the documentation of this file.

/*
 * Copyright 2008,2009,2010,2011,2012 Stephan Ossowski, Korbinian Schneeberger,
 * Felix Ott, Joerg Hagmann, Alf Scotland, Sebastian Bender
 * 
 * This file is part of SHORE.
 * 
 * SHORE is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * SHORE is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with SHORE.  If not, see <http://www.gnu.org/licenses/>.
 */


#ifndef SHORE_STATISTICS_BERRY_HPP__
#define SHORE_STATISTICS_BERRY_HPP__

#include "shore/statistics/glm.hpp"
#include "shore/base/mathops.hpp"

namespace shore {

class berry
{
 public:

        typedef long double obs_t;

 private:

        glm& m_lm;

        struct workspace
        {
                obs_t* Cs;
                obs_t* g0;

                obs_t* log_obs;
                obs_t* residuals;

                obs_t skewkurt[2];

                workspace(const int g,const size_t rows);
                ~workspace();

        } m_ws;

        int m_depth;
        int m_grid_size;
        obs_t m_conv_crit;
        obs_t m_min_num;

 public:

        berry(glm* lm,
          const int d=5,
          const int g=100,
          const obs_t cc=.001,
          const obs_t mn=1);

        template<class Iter>
        obs_t search(Iter observations)
        {
                const int DF=m_lm.get_residual_df();

                const obs_t min=
                        *std::min_element(observations,observations+m_lm.nrows());
                const obs_t max=
                        *std::max_element(observations,observations+m_lm.nrows());

                shore::seq(m_min_num-min,9.999*max,m_grid_size,m_ws.Cs);

                obs_t last_g0=0.0;
                obs_t last_c=0.0;

                for(int i=0;i<m_depth;++i)
                {
                        for(int j=0;j<m_grid_size;++j)
                        {
                                const obs_t c=m_ws.Cs[j];

                                for(size_t k=0;k<m_lm.nrows();++k)
                                {
                                        m_ws.log_obs[k]=shore::log2(observations[k]+c);
                                }

                                m_lm.fit(m_ws.log_obs);
                                m_lm.get_residuals(m_ws.residuals);

                                shore::skewness_kurtosis(m_ws.residuals,
                                        m_ws.residuals+m_lm.nrows(),m_ws.skewkurt,true);

                                m_ws.g0[j]=std::abs(m_ws.skewkurt[0])
                                        +std::abs(m_ws.skewkurt[1]+6.0/(DF+2));
                        }

                        const obs_t*const min_ptr=
                                std::min_element(m_ws.g0,m_ws.g0+m_grid_size);

                        const obs_t min_g0= *min_ptr;
                        const obs_t min_c=m_ws.Cs[min_ptr-m_ws.g0];

                        // convergence crit. reached?
                        if((i!=0)&&(std::abs(last_g0-min_g0)<=m_conv_crit))
                        {
                                return (last_g0<min_g0)?last_c:min_c;
                        }

                        last_g0=min_g0;
                        last_c=min_c;

                        if(min_ptr==m_ws.g0)
                        {
                                shore::seq(m_ws.Cs[0],m_ws.Cs[1],m_grid_size,m_ws.Cs);
                        }
                        else if(min_ptr==(&m_ws.g0[m_grid_size-1]))
                        {
                                shore::seq(m_ws.Cs[m_grid_size-2],
                                           m_ws.Cs[m_grid_size-1],m_grid_size,m_ws.Cs);
                        }
                        else
                        {
                                shore::seq(m_ws.Cs[min_ptr-m_ws.g0-1],
                                           m_ws.Cs[min_ptr-m_ws.g0+1],m_grid_size,m_ws.Cs);
                        }
                }

                return last_c;
        }
};

} // namespace

#endif // SHORE_STATISTICS_BERRY_HPP__