PoissonMixture.h

#ifndef theplu_yat_statistics_poisson_mixture_
#define theplu_yat_statistics_poisson_mixture_

// $Id: PoissonMixture.h 4181 2022-06-14 23:45:33Z peter $

/*
  Copyright (C) 2022 Peter Johansson

  This file is part of the yat library, https://dev.thep.lu.se/yat

  The yat library 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.

  The yat library 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 yat. If not, see <https://www.gnu.org/licenses/>.
*/

#include <yat/utility/Vector.h>

#include <map>

namespace theplu {
namespace yat {
  namespace utility {
    class Matrix;
  }
namespace statistics {

  class PoissonMixture
  {
  public:
    void add(unsigned long int k, unsigned long int n=1);

    void clear(void);

    void fit(size_t n);

    void fit(const yat::utility::VectorBase& m,
             const yat::utility::VectorBase& tau);

    double logL(void);

    double mean(size_t i) const;

    double tau(size_t i) const;
  private:
    std::map<unsigned long int, unsigned long int> count_;
    yat::utility::Vector tau_;
    yat::utility::Vector m_;

    // guess the parameters using a Gaussian mixture model
    void init_fit(size_t n);
    bool expectation_step(yat::utility::Matrix& H) const;

    // Find the optimal parameters
    void optimize(void);

    // return true if sum squared difference is > 1e-9
    bool optimize_tau(const yat::utility::Matrix& H);
    // return true if sum squared difference is > 1e-9
    bool optimize_m(const yat::utility::Matrix& H);
  };

}}}
#endif