simplenormaldistribution.h

Go to the documentation of this file.

/*

Copyright (c) 2003, Cornell University
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

   - Redistributions of source code must retain the above copyright notice,
       this list of conditions and the following disclaimer.
   - Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
   - Neither the name of Cornell University nor the names of its
       contributors may be used to endorse or promote products derived from
       this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.

*/

// -*- C++ -*-

#if !defined _CLUS_SIMPLENORMALDISTRIBUTION_H_
#define _CLUS_SIMPLENORMALDISTRIBUTION_H_

#include "distribution.h"
#include "multiclasscontinuousdistribution.h"
#include "statfct.h"

namespace CLUS
{

/** Implements a unidimensional normal distribution but the "active" dimension
    can be specified.
*/
class SimpleNormalDistribution : public Distribution
{
protected:
    /// expected value
    double mu;
    
    /// the variance of the distribution = sigma^2
    double variance;
    
    /// the active dimention
    int dimension;
    
    /// number of datapoints
    int N;
    
    /// sum of probabilities
    double S;
    
    /// sum prob*x
    double S_x;
    
    /// sum prob*x^2
    double S_x2;
    
    /// alpha/(sqrt(2PI)*sigma)
    double alpha_over_sigma;

public:
    SimpleNormalDistribution(int InDim = 0, int Dimension = 0):Distribution(InDim), mu(0.0), variance(1.0),
            dimension(Dimension), N(0), S(0.0), S_x(0.0), S_x2(0.0), alpha_over_sigma(1/sqrt(2*M_PI))
    { }

    SimpleNormalDistribution(int InDim, int Dimension, double Alpha, double Mu, double Sigma):Distribution(InDim),
            mu(Mu), dimension(Dimension), N(0), S(0.0), S_x(0.0), S_x2(0.0)
    {
        variance=pow2(Sigma);
        alpha_over_sigma=Alpha/(Sigma*sqrt(2*M_PI));
    }

    double LearnProbability(const double* DataCache)
    {
        if (weight==0.0)
            return 0.0;
        dataCache = DataCache;
        double X=dataCache[dimension]-mu;
        probabilityLearn = alpha_over_sigma*exp(-pow2(X)/(2*variance));
        return probabilityLearn;
    }

    /// InferProbability is not defined for this distribution
    double NormalizeLearnProbability(double Coef, int nrClus=1)
    {
        double pLearn=Distribution::NormalizeLearnProbability(Coef,nrClus);
        UpdateStatistics(dataCache, pLearn);

        return 0.0;
    }

    /// @todo to redefine these two since they are not virtual so the wrong methods are called
    double Probability(const double* DataCache)
    {
        return LearnProbability(DataCache);
    }

    double NormalizeProbability(double Coef, int nrClus=1)
    {
        return Distribution::NormalizeLearnProbability(Coef,nrClus);
    }

    void UpdateStatistics(const double* DataCache, double prob)
    {
        N++;
        S+=prob;
        S_x+=prob*DataCache[dimension];
        S_x2+=prob*pow2(DataCache[dimension]);
    }

    double UpdateParameters(void)
    {
        double oldMu, distP=0.0;

        // if the cluster isdead do nothing
        if (weight==0.0)
            return 0.0;

        if (S<=1.0)
        {
            cout << "Cluster got too small. We have to throw it out since is not anymore reliable." << endl;
            weight=0.0;
            distP=1.0;
            goto cleanup;
        }

        weight=S/N; // S is the number of points on which this clusterdepends

        oldMu = mu;
        mu = S_x/S;
        S_x2 = (S_x2-mu*S_x)/S;

        variance = S_x2;

        if (variance==0.0)
            weight=0;

        distP=pow2(oldMu-mu);

        alpha_over_sigma = weight / sqrt(2*M_PI*S_x2 );

cleanup:

        N=0;
        S=0.0;
        S_x=0.0;
        S_x2=0.0;

        return sqrt(distP)/(inDim+1); //how much the center has moved
    }

    void UpdateANOVAElements(int& numNonzero, double& SSR, double& n, double& sumS_x_overS, double& sumS_x)
    {
        if (weight!=0)
        {
            numNonzero++;
            SSR+=variance*S;
            n+=S;
            sumS_x_overS+=mu*S_x;
            sumS_x+=S_x;
        }
    }

    static string TypeName(void)
    {
        return string("SimpleNormalDistribution");
    }

    // @todo write RandomDistribution with reference to parent if needed
    void RandomDistribution(int NrClusters)
    {
        // reset weithts to 1.
        weight = 1.0;
        // choose mu random
        mu = RANDOM01FLOAT;
        variance = 1.0;
        alpha_over_sigma = 1 / sqrt(2*M_PI);
    }
    

#ifdef CLUS_USE_XML
    void PrintToXmlStream(ostream& output)
    {
        output << "<SimpleNormalDistribution";
        PrintAttribute(output, "mu", mu);
        PrintAttribute(output, "variance", variance);
        PrintAttribute(output, "dimension", dimension);
        PrintAttribute(output, "alpha_over_sigma", alpha_over_sigma);
        output << ">" << endl;

        Distribution::PrintToXmlStream(output);

        output << "</SimpleNormalDistribution>" << endl;
    }
#endif

};

double MulticlassContinuousDistribution<SimpleNormalDistribution>::PValueStatisticalTest(void)
{
    int numNonzero = 0;
    double SSR=0.0;
    double n=0.0;
    double sumS_x_overS=0.0;
    double sumS_x=0.0;
    for (int i=0; i<noClasses; i++)
    {
        distributions[i].UpdateANOVAElements(numNonzero, SSR, n, sumS_x_overS, sumS_x);
    }

    if (numNonzero>=2)
    {
        /*double SSA=sumS_x_overS-pow2(sumS_x)/n;
        double W=(SSA/(numNonzero-1))/(SSR/(n-numNonzero));*/

        /// @todo fix this: return F(W, noClasses-1, n-noClasses, 3);
        return 0.0;        
    }
    else
    {
        return 0.0;
    }
}

}

#endif // _CLUS_SIMPLENORMALDISTRIBUTION_H_

---