aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--distributions/component.hpp9
-rw-r--r--distributions/kmm_component.hpp205
-rw-r--r--misc/optvec.hpp18
-rw-r--r--statistics/kmm_stat.hpp101
4 files changed, 323 insertions, 10 deletions
diff --git a/distributions/component.hpp b/distributions/component.hpp
index d761bf8..ee13adb 100644
--- a/distributions/component.hpp
+++ b/distributions/component.hpp
@@ -46,7 +46,14 @@ namespace opt_utilities
{
for(int i=0;i<rhs.components.size();++i)
{
- add_component(*rhs.components[i],rhs.get_param_info(rhs.weight_num[i]).get_value());
+ if(i>0)
+ {
+ add_component(*rhs.components[i],rhs.get_param_info(rhs.weight_num[i-1]).get_value());
+ }
+ else
+ {
+ add_component(*rhs.components[i]);
+ }
}
for(int i=0;i<rhs.get_num_params();++i)
{
diff --git a/distributions/kmm_component.hpp b/distributions/kmm_component.hpp
new file mode 100644
index 0000000..0e77508
--- /dev/null
+++ b/distributions/kmm_component.hpp
@@ -0,0 +1,205 @@
+/**
+ \file kmm_component.hpp
+ \brief rpresents a distribution composed of more than one kmm_components
+ \author Junhua Gu
+ */
+
+
+#ifndef KMM_COMPONENT_MODEL_H_
+#define KMM_COMPONENT_MODEL_H_
+#define OPT_HEADER
+#include <core/fitter.hpp>
+#include <cmath>
+#include <misc/optvec.hpp>
+#include <sstream>
+#include <iostream>
+
+
+
+namespace opt_utilities
+{
+ template <typename T>
+ class kmm_component
+ :public model<optvec<T>,optvec<T>,optvec<T>,std::string>
+ {
+ private:
+ std::vector<model<optvec<T>,optvec<T>,optvec<T>,std::string>*> components;
+ std::vector<int> weight_num;
+ private:
+ kmm_component* do_clone()const
+ {
+ return new kmm_component<T>(*this);
+ }
+
+ const char* do_get_type_name()const
+ {
+ return "Multi-distribution";
+ }
+ public:
+ kmm_component()
+ {
+ //this->push_param_info(param_info<optvec<T> >("x0",0));
+ //this->push_param_info(param_info<optvec<T> >("sigma",1));
+ }
+
+ kmm_component(const kmm_component& rhs)
+ {
+ for(int i=0;i<rhs.components.size();++i)
+ {
+ if(i>=1)
+ {
+ add_component(*rhs.components[i],rhs.get_param_info(rhs.weight_num[i-1]).get_value());
+ }
+ else
+ {
+ add_component(*rhs.components[i]);
+ }
+ }
+ for(int i=0;i<rhs.get_num_params();++i)
+ {
+ set_param_info(rhs.get_param_info(i));
+ }
+
+ }
+
+ kmm_component& operator=(const kmm_component& rhs)
+ {
+ if(this==&rhs)
+ {
+ return rhs;
+ }
+ for(int i=0;i<components.size();++i)
+ {
+ components[i]->destroy();
+ }
+ components.clear();
+ for(int i=0;i<rhs.components.size();++i)
+ {
+ add_component(*rhs.components[i]);
+ }
+ for(int i=0;i<rhs.get_num_params();++i)
+ {
+ set_param_info(rhs.get_param_info(i));
+ }
+
+ }
+
+ ~kmm_component()
+ {
+ for(int i=0;i<components.size();++i)
+ {
+ components[i]->destroy();
+ }
+ }
+
+
+ public:
+ void add_component(const model<optvec<T>,optvec<T>,optvec<T>,std::string>& m,const T& w=0)
+ {
+ int morder=components.size();
+ components.push_back(m.clone());
+ std::ostringstream oss;
+ oss<<morder;
+ std::string smorder=oss.str();
+ if(components.size()>1)
+ {
+ param_info<optvec<T>,std::string> p1(std::string("_w")+smorder,w);
+ this->push_param_info(p1);
+ weight_num.push_back(this->get_num_params()-1);
+ }
+ int np=m.get_num_params();
+ for(int i=0;i<np;++i)
+ {
+ param_info<optvec<T>,std::string> p(m.get_param_info(i));
+ param_info<optvec<T>,std::string> p1(p.get_name()+smorder,p.get_value());
+ this->push_param_info(p1);
+ }
+ }
+
+ optvec<T> convert_unit_sph(const optvec<T>& p)
+ {
+ int ndim=p.size()+1;
+ optvec<double> result(ndim);
+ result[0]=1;
+ for(int i=0;i<p.size();++i)
+ {
+ for(int j=0;j<=i;++j)
+ {
+ result[j]*=cos(p[i]);
+ }
+ result[i+1]=sin(p[i]);
+ }
+ for(int i=0;i<result.size();++i)
+ {
+ result[i]=result[i]*result[i];
+ }
+ return result;
+ }
+
+
+ optvec<T> do_eval(const optvec<T>& x,const optvec<T>& param)
+ {
+ T result(0);
+ optvec<T> weight_angle;
+ for(int i=0;i<weight_num.size();++i)
+ {
+ weight_angle.push_back(param[weight_num[i]]);
+ }
+ optvec<T> weight(convert_unit_sph(weight_angle));
+
+ int pnum=0;
+ for(int i=0;i<components.size();++i)
+ {
+ T temp_result=0;
+ optvec<T> p(components[i]->get_num_params());
+ for(int j=0;j<p.size();++j)
+ {
+ p[j]=param[pnum++];
+ }
+ result+=(components[i]->eval(x,p)[0]*weight[i]);
+ ++pnum;
+ }
+ optvec<T> result1(1);
+ result1[0]=result;;
+ return result1;
+ }
+
+ optvec<T> eval_unsumed(const optvec<T>& x,const optvec<T>& param)
+ {
+ optvec<T> weight_angle;
+ for(int i=0;i<weight_num.size();++i)
+ {
+ weight_angle.push_back(param[weight_num[i]]);
+ }
+ optvec<T> weight(convert_unit_sph(weight_angle));
+ optvec<T> result(components.size());
+ int pnum=0;
+ for(int i=0;i<components.size();++i)
+ {
+ T temp_result=0;
+ optvec<T> p(components[i]->get_num_params());
+ for(int j=0;j<p.size();++j)
+ {
+ p[j]=param[pnum++];
+ }
+ result[i]=(components[i]->eval(x,p)[0]*weight[i]);
+ ++pnum;
+ }
+
+ return result;
+ }
+
+
+
+ private:
+ std::string do_get_information()const
+ {
+ return "";
+ }
+ };
+}
+
+
+
+#endif
+//EOF
diff --git a/misc/optvec.hpp b/misc/optvec.hpp
index e52e853..7967e5c 100644
--- a/misc/optvec.hpp
+++ b/misc/optvec.hpp
@@ -53,7 +53,7 @@ namespace opt_utilities
optvec<T> operator+(const optvec<T>& x1,const optvec<T>& x2)
{
- optvec<T> result(min(x1.size(),x2.size()));
+ optvec<T> result(std::min(x1.size(),x2.size()));
for(size_t i=0;i!=result.size();++i)
{
result[i]=x1[i]+x2[i];
@@ -76,7 +76,7 @@ namespace opt_utilities
template<typename T>
optvec<T>& operator+=(optvec<T>& x1,const optvec<T>& x2)
{
- for(size_t i=0;i!=min(x1.size(),x2.size());++i)
+ for(size_t i=0;i!=std::min(x1.size(),x2.size());++i)
{
x1[i]+=x2[i];
}
@@ -87,7 +87,7 @@ namespace opt_utilities
optvec<T> operator-(const optvec<T>& x1,const optvec<T>& x2)
{
- optvec<T> result(min(x1.size(),x2.size()));
+ optvec<T> result(std::min(x1.size(),x2.size()));
for(size_t i=0;i!=result.size();++i)
{
result[i]=x1[i]-x2[i];
@@ -123,7 +123,7 @@ namespace opt_utilities
template<typename T>
optvec<T>& operator-=(optvec<T>& x1,const optvec<T>& x2)
{
- for(size_t i=0;i!=min(x1.size(),x2.size());++i)
+ for(size_t i=0;i!=std::min(x1.size(),x2.size());++i)
{
x1[i]-=x2[i];
}
@@ -134,7 +134,7 @@ namespace opt_utilities
optvec<T> operator*(const optvec<T>& x1,const optvec<T>& x2)
{
- optvec<T> result(min(x1.size(),x2.size()));
+ optvec<T> result(std::min(x1.size(),x2.size()));
for(size_t i=0;i!=result.size();++i)
{
result[i]=x1[i]*x2[i];
@@ -170,7 +170,7 @@ namespace opt_utilities
template<typename T>
optvec<T>& operator*=(optvec<T>& x1,const optvec<T>& x2)
{
- for(size_t i=0;i!=min(x1.size(),x2.size());++i)
+ for(size_t i=0;i!=std::min(x1.size(),x2.size());++i)
{
x1[i]*=x2[i];
}
@@ -192,7 +192,7 @@ namespace opt_utilities
optvec<T> operator/(const optvec<T>& x1,const optvec<T>& x2)
{
- optvec<T> result(min(x1.size(),x2.size()));
+ optvec<T> result(std::min(x1.size(),x2.size()));
for(size_t i=0;i!=result.size();++i)
{
result[i]=x1[i]/x2.at(i);
@@ -227,7 +227,7 @@ namespace opt_utilities
template<typename T>
optvec<T>& operator/=(optvec<T>& x1,const optvec<T>& x2)
{
- for(size_t i=0;i!=min(x1.size(),x2.size());++i)
+ for(size_t i=0;i!=std::min(x1.size(),x2.size());++i)
{
x1[i]/=x2.at(i);
}
@@ -258,7 +258,7 @@ namespace opt_utilities
template<typename T>
bool operator<(const optvec<T>& x1,const optvec<T>& x2)
{
- for(size_t i=0;i!=min(x1.size(),x2.size());++i)
+ for(size_t i=0;i!=std::min(x1.size(),x2.size());++i)
{
if(x1[i]!=x2[i])
{
diff --git a/statistics/kmm_stat.hpp b/statistics/kmm_stat.hpp
new file mode 100644
index 0000000..0a205bd
--- /dev/null
+++ b/statistics/kmm_stat.hpp
@@ -0,0 +1,101 @@
+/**
+ \file kmm_stat.hpp
+ \brief maximum-liklihood statistic
+ \author Junhua Gu
+ */
+
+#ifndef KMM_STAT_HPP
+#define KMM_STAT_HPP
+#define OPT_HEADER
+#include <core/fitter.hpp>
+#include <math/vector_operation.hpp>
+#include <distributions/kmm_component.hpp>
+#include <iostream>
+#include <cassert>
+
+
+using std::cout;using std::endl;
+namespace opt_utilities
+{
+
+ /**
+ \brief c-statistic, max-likelihood method
+ \tparam Ty the return type of model
+ \tparam Tx the type of the self-var
+ \tparam Tp the type of model parameter
+ \tparam Ts the type of the statistic
+ \tparam Tstr the type of the string used
+ */
+ template<typename T,typename Ts,typename Tstr>
+ class kmm_stat
+ :public statistic<optvec<T>,optvec<T>,optvec<T>,Ts,Tstr>
+ {
+ private:
+ bool verb;
+ int n;
+ public:
+ kmm_stat()
+ :verb(true)
+ {}
+
+ void verbose(bool v)
+ {
+ verb=v;
+ }
+
+ const char* do_get_type_name()const
+ {
+ return "kmm likelihood";
+ }
+
+ public:
+
+ statistic<optvec<T>,optvec<T>,optvec<T>,Ts,Tstr>* do_clone()const
+ {
+ // return const_cast<statistic<Ty,Tx,Tp>*>(this);
+ return new kmm_stat<T,Ts,Tstr>(*this);
+ }
+
+ Ts do_eval(const optvec<T>& p)
+ {
+ Ts result(0);
+
+ kmm_component<T>* kmm=
+ dynamic_cast<kmm_component<T>* >(&const_cast<model<optvec<T>,optvec<T>,optvec<T>,std::string>&>(this->get_fitter().get_model()));
+
+ if(kmm==0)
+ {
+ throw opt_exception("model is not a kmm component");
+ }
+
+ for(int i=(this->get_data_set()).size()-1;i>=0;--i)
+ {
+ optvec<T> unsummed_possibility(kmm->eval_unsumed(this->get_data_set().get_data(i).get_x(),p));
+ T maxp=*max_element(unsummed_possibility.begin(),unsummed_possibility.end());
+ T logp=std::log(maxp);
+ result-=this->get_data_set().get_data(i).get_y()[0]*logp;
+ }
+
+ if(verb)
+ {
+ n++;
+ if(n%10==0)
+ {
+ cout<<"a:"<<result<<"\t";
+ for(size_t i=0;i<get_size(p);++i)
+ {
+ cout<<get_element(p,i)<<",";
+ }
+ cout<<endl;
+ }
+
+ }
+ return result;
+ }
+ };
+}
+
+#endif
+//EOF
+
+
Copyright © 2017-2018 Aaron LI. All Rights Reserved.