extra_newhb.h

 
#define q_factor_val 2.0
 
namespace Pds{
namespace LogisticModel{
 
double CostXqEntropy( const Pds::Vector &W,
                       const Pds::Matrix &X,
                       const Pds::Vector &Y,
                       double q);
 
Pds::Vector GradientCostXqEntropy( const Pds::Vector &W,
                                    const Pds::Matrix &X,
                                    const Pds::Vector &Y,
                                    double q,
                                    double h);
 
Pds::Vector GradientCostXqEntropySoft(  const Pds::Vector &W,
                                        const Pds::Matrix &X,
                                        const Pds::Vector &Y,
                                        double q);
                                    
 
}//LogisticModel
}//Pds
 
 
 
 
 
double Pds::LogisticModel::CostXqEntropy( const Pds::Vector &W,
                                           const Pds::Matrix &X,
                                           const Pds::Vector &Y,
                                           double q)
{
    if(W.IsEmpty())             return Pds::Ra::Nan;
    if(X.IsEmpty())             return Pds::Ra::Nan;
    if(Y.IsEmpty())             return Pds::Ra::Nan;
    
    if(X.Nlin()!=Y.Nlin())      return Pds::Ra::Nan;
    if(X.Ncol()!=(W.Nlin()-1))  return Pds::Ra::Nan;
    
    Pds::Vector Yh=Y.Geq(0.5);
    
    Pds::Vector z=Pds::LogisticModel::Classify(W,X);
    z=z.Geq(0.5);
    
    unsigned int L=Yh.Nlin();
    
    double Np=Yh.Sum();
    double Na=z.Sum();
    double N1=Yh.Dot(z);
    
    // New entropy 1
    double E1;
    if(0==Na)   E1=0;
    else        E1=Pds::XqHb(N1/Na,q);
    
    // New entropy 0
    double E0;
    double P=(Np-N1)/(L-Na); 
    if(P>1.0)   P=1.0;
    if(P<0)     P=0.0; 
    if(L==Na)   E0=0.0;
    else        E0=Pds::XqHb(P,q);
 
    
    //return pow(Na/L,2)*E1+pow((L-Na)/L,2)*E0;
    return (Na/L)*E1+((L-Na)/L)*E0;
}
                                                        
Pds::Vector Pds::LogisticModel::GradientCostXqEntropy( const Pds::Vector &W,
                                                        const Pds::Matrix &X,
                                                        const Pds::Vector &Y,
                                                        double q,
                                                        double h)
{
    if(W.IsEmpty())             return Pds::Vector();
    if(X.IsEmpty())             return Pds::Vector();
    if(Y.IsEmpty())             return Pds::Vector();
    
    if(X.Nlin()!=Y.Nlin())      return Pds::Vector();
    if(X.Ncol()!=(W.Nlin()-1))  return Pds::Vector();
    
    if(h==0.0)  return Pds::Vector();
    
    unsigned int N=W.Nel();
    
    Pds::Vector D(N);
    Pds::Vector Whp(N);
    Pds::Vector Whm(N);
    double dn;
    
    Whp.Copy(W);
    Whm.Copy(W);
 
    for(unsigned int n=0;n<N;n++)
    {
        Whp.SetRaw(n,W.GetRaw(n)+h);
        Whm.SetRaw(n,W.GetRaw(n)-h);
 
        dn=( Pds::LogisticModel::CostXqEntropy(Whp,X,Y,q)
            -Pds::LogisticModel::CostXqEntropy(Whm,X,Y,q) )/(2*h);
                    
        D.SetRaw(n,dn);
        
        Whp.SetRaw(n,W.GetRaw(n));
        Whm.SetRaw(n,W.GetRaw(n));
    }
    return D;
}
 
Pds::Vector Pds::LogisticModel::GradientCostXqEntropySoft(  const Pds::Vector &W,
                                                            const Pds::Matrix &X,
                                                            const Pds::Vector &Y,
                                                            double q)
{
    if(W.IsEmpty())             return Pds::Vector();
    if(X.IsEmpty())             return Pds::Vector();
    if(Y.IsEmpty())             return Pds::Vector();
    
    if(X.Nlin()!=Y.Nlin())      return Pds::Vector();
    if(X.Ncol()!=(W.Nlin()-1))  return Pds::Vector();
    
    Pds::Vector Yt=Y.Geq(0.5);
    unsigned int N=Yt.Sum();// numero total de 1s
    unsigned int L=Yt.Nel();// numero de elementos
    
    Pds::Matrix R=Pds::RegressorMatrix(X);
    Pds::Vector z=Pds::Sigmoid(R*W);
    Pds::Vector Dz=z.Product(1.0-z);
    
    double Na=z.Sum();
    double N1=Yt.Dot(z);
    double p1=N1/Na;
    double p2=(N-N1)/(L-Na);
    
    double h1=Pds::XqHb(p1,q);
    double h2=Pds::XqHb(p2,q);
    double Dh1=Pds::DXqHb(p1,q);
    double Dh2=Pds::DXqHb(p2,q);
    
    double factor1=-(-h1+h2+p1*Dh1-p2*Dh2)/L;
    double factor2=-(-Dh1+Dh2)/L;
    Pds::Vector Der=  factor1*R.TMul(Dz)+factor2*R.TMul(Dz.Product(Yt));
    
    return Der;
}