added john platt's method for computung calibrated posterior outputs for

Binary Labels
shogun-toolbox · Aug 14, 2012 · 47bbe6a · 47bbe6a
1 parent d16ae4b
commit 47bbe6a
Show file tree

Hide file tree

Showing 2 changed files with 182 additions and 3 deletions.
diff --git a/src/shogun/labels/BinaryLabels.cpp b/src/shogun/labels/BinaryLabels.cpp
@@ -73,3 +73,170 @@ 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_DEBUG("line search fails, A=%f, B=%f, g1=%f, g2=%f, dA=%f, "
+			"dB=%f, gd=%f\n", a, b, g1, g2, dA, dB, gd);
+		}
+	}
+
+	if (it>=maxiter-1)
+	{
+		SG_DEBUG("reaching maximal iterations, g1=%f, g2=%f\n", 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,11 @@ 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.
  */
 class CBinaryLabels : public CDenseLabels
 {
@@ -74,10 +79,17 @@ 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.
+		 *
+		 * 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