VineCopulaFactory: Select native

Author: jschueller

lib/src/VineCopula.cxx

@@ -20,6 +20,9 @@
  */
 #include "otvine/VineCopula.hxx"
 #include <openturns/PersistentObjectFactory.hxx>
+#include <openturns/ClaytonCopula.hxx>
+#include <openturns/GumbelCopula.hxx>
+#include <openturns/FrankCopula.hxx>
 
 #include <vinecopulib.hpp>
 
@@ -137,6 +140,28 @@ String VineCopula::__repr__() const
   return oss;
 }
 
+Distribution VineCopula::asNative() const
+{
+  if (p_vinecop_->get_dim() != 2)
+    throw NotYetImplementedException(HERE) << "dim is not 2";
+  std::vector<std::vector<vinecopulib::Bicop> > pairs = p_vinecop_->get_all_pair_copulas();
+  if (pairs.size() != 1 && pairs[0].size() != 1)
+    throw NotYetImplementedException(HERE) << "tree contains more than one node";
+  vinecopulib::Bicop bicop = pairs[0][0];
+  int rotation = bicop.get_rotation();
+  if (rotation != 0)
+    throw NotYetImplementedException(HERE) << "rotation is not 0";
+  const Scalar theta = bicop.get_parameters()(0, 0);
+  if (bicop.get_family_name() == "Clayton")
+    return ClaytonCopula(theta);
+  else if (bicop.get_family_name() == "Frank")
+    return FrankCopula(theta);
+  else if (bicop.get_family_name() == "Gumbel")
+    return GumbelCopula(theta);
+  else
+    throw NotYetImplementedException(HERE) << "unknown family: " << bicop.get_family_name();
+}
+
 /* Method save() stores the object through the StorageManager */
 void VineCopula::save(Advocate & adv) const
 {

lib/src/VineCopulaFactory.cxx

@@ -47,7 +47,7 @@ VineCopulaFactory * VineCopulaFactory::clone() const
 }
 
 /* example of a func that return a point squared. */
-Distribution VineCopulaFactory::build(const Sample & sample) const
+VineCopula VineCopulaFactory::buildAsNative(const Sample & sample) const
 {
   const UnsignedInteger dimension = sample.getDimension();
   const UnsignedInteger size = sample.getSize();
@@ -56,10 +56,23 @@ Distribution VineCopulaFactory::build(const Sample & sample) const
     for (UnsignedInteger j = 0; j < dimension; ++ j)
       data(i, j) = sample(i, j);
   Pointer <vinecopulib::Vinecop> p_vinecop = new vinecopulib::Vinecop(dimension);
-  p_vinecop->select(data);
+  std::vector<vinecopulib::BicopFamily> family_set;
+#if 1
+  family_set = {
+    vinecopulib::BicopFamily::clayton,
+    vinecopulib::BicopFamily::gumbel,
+    vinecopulib::BicopFamily::frank};
+#endif
+  vinecopulib::FitControlsVinecop controls(family_set);
+  p_vinecop->select(data, controls);
   return VineCopula(p_vinecop);
 }
 
+Distribution VineCopulaFactory::build(const Sample & sample) const
+{
+  return buildAsNative(sample);
+}
+
 /* String converter */
 String VineCopulaFactory::__repr__() const
 {

lib/src/otvine/VineCopula.hxx

@@ -55,7 +55,10 @@ public:
 protected:
   OT::Bool equals(const DistributionImplementation & other) const override;
 public:
-  
+
+  /** Convert to native OT distribution if possible */
+  OT::Distribution asNative() const;
+
   /** String converter */
   OT::String __repr__() const override;
 

lib/src/otvine/VineCopulaFactory.hxx

@@ -48,6 +48,7 @@ public:
 
   /** example of a func that return a point squared. **/
   OT::Distribution build(const OT::Sample & sample) const;
+  VineCopula buildAsNative(const OT::Sample & sample) const;
 
   /** String converter */
   OT::String __repr__() const override;

python/src/VineCopula.i

@@ -1,11 +1,11 @@
-// SWIG file VineCopulaFactory.i
+// SWIG file VineCopula.i
 
 %{
-#include "otvine/VineCopulaFactory.hxx"
+#include "otvine/VineCopula.hxx"
 %}
 
-%include VineCopulaFactory_doc.i
+%include VineCopula_doc.i
 
-%copyctor OTVINE::VineCopulaFactory;
+%copyctor OTVINE::VineCopula;
 
-%include otvine/VineCopulaFactory.hxx
+%include otvine/VineCopula.hxx

python/test/t_VineCopulaFactory_std.py

@@ -16,5 +16,13 @@
 realization = distribution.getRealization()
 print(realization)
 
-sample = distribution.getSample(1000)
-print(sample.computeMean())
+bic = ot.FittingTest.BIC(sample, distribution)
+print("bic=", bic)
+
+sample2 = distribution.getSample(1000)
+print(sample2.computeMean())
+
+factory2 = ot.ClaytonCopulaFactory()
+print(factory2.build(sample))
+
+print(factory.buildAsNative(sample).asNative())