Merge pull request #714 from karlnapf/master

john platt's sigmoid fitting for binary svm
shogun-toolbox · Aug 14, 2012 · 419f224 · 419f224
2 parents d16ae4b + f4c19cb
commit 419f224
Show file tree

Hide file tree

Showing 4 changed files with 230 additions and 3 deletions.
diff --git a/examples/undocumented/libshogun/Makefile b/examples/undocumented/libshogun/Makefile
@@ -54,6 +54,7 @@ TARGETS = basic_minimal \
 		  features_copy_subset_string_features \
 		  features_copy_subset_sparse_features \
 		  features_create_merged_copy \
+		  labels_binary_fit_sigmoid \
 		  mathematics_confidence_intervals \
 		  clustering_kmeans base_parameter_map \
 		  base_load_file_parameters \

diff --git a/examples/undocumented/libshogun/labels_binary_fit_sigmoid.cpp b/examples/undocumented/libshogun/labels_binary_fit_sigmoid.cpp
@@ -0,0 +1,43 @@
+/*
+ * This program 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.
+ *
+ * Written (W) 2012 Heiko Strathmann
+ */
+#include <shogun/labels/BinaryLabels.h>
+
+using namespace shogun;
+
+void test_sigmoid_fitting()
+{
+	CBinaryLabels* labels=new CBinaryLabels(10);
+	labels->set_confidences(SGVector<float64_t>(labels->get_num_labels()));
+
+	for (index_t i=0; i<labels->get_num_labels(); ++i)
+		labels->set_confidence(i%2==0 ? 1 : -1, i);
+
+	labels->get_confidences().display_vector("scores");
+	labels->scores_to_probabilities();
+
+	/* only two probabilities will be the result, repeatedly,
+	 * assert against reference implementation */
+	ASSERT(CMath::abs(labels->get_confidence(0)-0.8571428439385661)<10E-15);
+	ASSERT(CMath::abs(labels->get_confidence(1)-0.14285715606143384)<10E-15);
+
+	SG_UNREF(labels);
+}
+
+int main()
+{
+	init_shogun_with_defaults();
+
+//	sg_io->set_loglevel(MSG_DEBUG);
+
+	test_sigmoid_fitting();
+
+	exit_shogun();
+	return 0;
+}
+
diff --git a/src/shogun/labels/BinaryLabels.cpp b/src/shogun/labels/BinaryLabels.cpp
@@ -73,3 +73,172 @@ ELabelType CBinaryLabels::get_label_type()
 {
 	return LT_BINARY;
 }
+
+void CBinaryLabels::scores_to_probabilities()
+{
+	SG_DEBUG("entering CBinaryLabels::scores_to_probabilities()\n");
+
+	REQUIRE(m_confidences.vector, "%s::scores_to_probabilities() requires "
+			"confidences vector!\n", get_name());
+
+	/* count prior0 and prior1 if needed */
+	int32_t prior0=0;
+	int32_t prior1=0;
+	SG_DEBUG("counting number of positive and negative labels\n");
+	{
+		for (index_t i=0; i<m_confidences.vlen; ++i)
+		{
+			if (m_confidences[i]>0)
+				prior1++;
+			else
+				prior0++;
+		}
+	}
+	SG_DEBUG("%d pos; %d neg\n", prior1, prior0);
+
+	/* parameter setting */
+	/* maximum number of iterations */
+	index_t maxiter=100;
+
+	/* minimum step taken in line search */
+	float64_t minstep=1E-10;
+
+	/* for numerically strict pd of hessian */
+	float64_t sigma=1E-12;
+	float64_t eps=1E-5;
+
+	/* construct target support */
+	float64_t hiTarget=(prior1+1.0)/(prior1+2.0);
+	float64_t loTarget=1/(prior0+2.0);
+	index_t length=prior1+prior0;
+
+	SGVector<float64_t> t(length);
+	for (index_t i=0; i<length; ++i)
+	{
+		if (m_confidences[i]>0)
+			t[i]=hiTarget;
+		else
+			t[i]=loTarget;
+	}
+
+	/* initial Point and Initial Fun Value */
+	/* result parameters of sigmoid */
+	float64_t a=0;
+	float64_t b=CMath::log((prior0+1.0)/(prior1+1.0));
+	float64_t fval=0.0;
+
+	for (index_t i=0; i<length; ++i)
+	{
+		float64_t fApB=m_confidences[i]*a+b;
+		if (fApB>=0)
+			fval+=t[i]*fApB+CMath::log(1+CMath::exp(-fApB));
+		else
+			fval+=(t[i]-1)*fApB+CMath::log(1+CMath::exp(fApB));
+	}
+
+	index_t it;
+	float64_t g1;
+	float64_t g2;
+	for (it=0; it<maxiter; ++it)
+	{
+		SG_DEBUG("Iteration %d, a=%f, b=%f, fval=%f\n", it, a, b, fval);
+
+		/* Update Gradient and Hessian (use H' = H + sigma I) */
+		float64_t h11=sigma; //Numerically ensures strict PD
+		float64_t h22=h11;
+		float64_t h21=0;
+		g1=0;
+		g2=0;
+
+		for (index_t i=0; i<length; ++i)
+		{
+			float64_t fApB=m_confidences[i]*a+b;
+			float64_t p;
+			float64_t q;
+			if (fApB>=0)
+			{
+				p=CMath::exp(-fApB)/(1.0+CMath::exp(-fApB));
+				q=1.0/(1.0+CMath::exp(-fApB));
+			}
+			else
+			{
+				p=1.0/(1.0+CMath::exp(fApB));
+				q=CMath::exp(fApB)/(1.0+CMath::exp(fApB));
+			}
+
+			float64_t d2=p*q;
+			h11+=m_confidences[i]*m_confidences[i]*d2;
+			h22+=d2;
+			h21+=m_confidences[i]*d2;
+			float64_t d1=t[i]-p;
+			g1+=m_confidences[i]*d1;
+			g2+=d1;
+		}
+
+		/* Stopping Criteria */
+		if (CMath::abs(g1)<eps && CMath::abs(g2)<eps)
+			break;
+
+		/* Finding Newton direction: -inv(H') * g */
+		float64_t det=h11*h22-h21*h21;
+		float64_t dA=-(h22*g1-h21*g2)/det;
+		float64_t dB=-(-h21*g1+h11*g2)/det;
+		float64_t gd=g1*dA+g2*dB;
+
+		/* Line Search */
+		float64_t stepsize=1;
+
+		while (stepsize>=minstep)
+		{
+			float64_t newA=a+stepsize*dA;
+			float64_t newB=b+stepsize*dB;
+
+			/* New function value */
+			float64_t newf=0.0;
+			for (index_t i=0; i<length; ++i)
+			{
+				float64_t fApB=m_confidences[i]*newA+newB;
+				if (fApB>=0)
+					newf+=t[i]*fApB+CMath::log(1+CMath::exp(-fApB));
+				else
+					newf+=(t[i]-1)*fApB+CMath::log(1+CMath::exp(fApB));
+			}
+
+			/* Check sufficient decrease */
+			if (newf<fval+0.0001*stepsize*gd)
+			{
+				a=newA;
+				b=newB;
+				fval=newf;
+				break;
+			}
+			else
+				stepsize=stepsize/2.0;
+		}
+
+		if (stepsize<minstep)
+		{
+			SG_WARNING("%s::scores_to_probabilities(): line search fails, A=%f, "
+					"B=%f, g1=%f, g2=%f, dA=%f, dB=%f, gd=%f\n",
+					get_name(), a, b, g1, g2, dA, dB, gd);
+		}
+	}
+
+	if (it>=maxiter-1)
+	{
+		SG_WARNING("%s::scores_to_probabilities(): reaching maximal iterations,"
+				" g1=%f, g2=%f\n", get_name(), g1, g2);
+	}
+
+	SG_DEBUG("fitted sigmoid: a=%f, b=%f\n", a, b);
+
+	/* now the sigmoid is fitted, convert all confidences to probabilities */
+	for (index_t i=0; i<m_confidences.vlen; ++i)
+	{
+		float64_t fApB=m_confidences[i]*a+b;
+		m_confidences[i]=fApB>=0 ? CMath::exp(-fApB)/(1.0+exp(-fApB)) :
+				1.0/(1+CMath::exp(fApB));
+	}
+
+	SG_DEBUG("leaving CBinaryLabels::scores_to_probabilities()\n");
+}
diff --git a/src/shogun/labels/BinaryLabels.h b/src/shogun/labels/BinaryLabels.h
@@ -6,7 +6,7 @@
  *
  * Written (W) 1999-2009 Soeren Sonnenburg
  * Written (W) 1999-2008 Gunnar Raetsch
- * Written (W) 2011 Heiko Strathmann
+ * Written (W) 2011-2012 Heiko Strathmann
  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
  */
 
@@ -26,6 +26,12 @@ namespace shogun
 /** @brief Binary Labels for binary classification 
  *
  * valid values for labels are +1/-1
+ *
+ * Confidences may be converted into calibrated probabilities using
+ * scores_to_probabilities(), which implements the method described in
+ * Lin, H., Lin, C., and Weng, R. (2007).
+ * A note on Platt's probabilistic outputs for support vector machines.
+ * Should only be used in conjunction with SVM.
  */
 class CBinaryLabels : public CDenseLabels
 {
@@ -74,10 +80,18 @@ class CBinaryLabels : public CDenseLabels
 		 */
 		virtual ELabelType get_label_type();
 
+		/** Converts all confidences to calibrated probabilities by fitting a
+		 * sigmoid function using the method described in
+		 * Lin, H., Lin, C., and Weng, R. (2007).
+		 * A note on Platt's probabilistic outputs for support vector machines.
+		 *
+		 * Should only be used in conjunction with SVM.
+		 * The fitted sigmoid is used to replace all confidence values.
+		 */
+		void scores_to_probabilities();
+
 		/** @return object name */
 		inline virtual const char* get_name() const { return "BinaryLabels"; }
-
-
 };
 }
 #endif