git-svn-id: file:///home/svn/opt_utilities@197 ed2142bd-67ad-457f-ba7c-d818d4011675

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

+

+