rand.c

数据挖掘经典的hierarchial clustering algorithm

    count++;
    sum += value;
    sumOfSqu += value*value;
    if (value > sampleMax) sampleMax = value;
    if (value < sampleMin) sampleMin = value;
    };

void FloatRandomStream::CollectStat(float value) {
    count++;
    sum += value;
    sumOfSqu += value*value;
    if (value > sampleMax) sampleMax = value;
    if (value < sampleMin) sampleMin = value;
    };

DiscreteUniform::DiscreteUniform(int n, int stream) {
    this->stream = stream;
    this->n   = n;
    statFlag  = 2;
    mean      = (n+1)*0.5;
    if (n == 1){
        std = 0;
    } else {
        std = sqrt((n*n-1)/(12.0));
    };
    };

int  DiscreteUniform::Draw() {
    int result;
    float sample = rint(n)*Uniform(this->stream);
#if defined(linux ) || (defined(i386) && defined(sun))
    result = int(rint(floor(sample))) + 1;  // I hope that was the right idea
#else
    result = irint(floor(sample)) + 1;
#endif
    if (statFlag  == 2)  CollectStat(result);
    return result;
    };

void DiscreteUniform::PrintProfile() {
        printf("\n DiscreteUniform [1-%d] generator (%f, %f)",n,mean,std);
    };

Choose::Choose(float prob, int stream) {
    this->stream = stream;
    this->prob   = prob;
    statFlag     = 2;
    mean         = prob;
    std          = sqrt(prob*(1-prob));
    };

int Choose::Draw() {
    int sample;
    sample = (Uniform(this->stream) <= prob) ? 1 : 0;
    if (statFlag  & 2) CollectStat(sample);
    return sample;
    };

void Choose::PrintProfile() {
    printf("\n Choose[%f] generator (%f, %f)",prob,mean,std);
    }

Uniform::Uniform(float low, float high, int stream) {
    this->stream = stream;
    this->low  = low;
    this->high = high;
    statFlag   = 2;
    mean       = (low+high)*0.5;
    std        = (high-low)/sqrt(12.0);
    };

float   Uniform::Draw() {
    double sample = low+(high-low)*Uniform(this->stream);
    if (statFlag  == 2) CollectStat(sample);
    return sample;
    };

void Uniform::PrintProfile() {
        printf("\n Uniform [%f,%f] generator (%f, %f)",low,high,mean,std);
    };

Uniform01::Uniform01(int stream) {
    this->stream = stream;
    statFlag     = 2;
    mean         = 0.5;
    std          = sqrt(1.0/12.0);
    };

float   Uniform01::Draw() {
    double sample = Uniform(this->stream);
    if (statFlag  == 2) CollectStat(sample);
    return sample;
    };

void Uniform01::PrintProfile() {
        printf("\n Uniform [0,1] generator (%f, %f)",mean,std);
    };

Exponential::Exponential(float mean, int stream) {
        this->stream = stream;
        this->mean = mean;
        this->std = mean;
        statFlag  = 2;
}

float   Exponential::Draw()
{
    float sample = -mean * ( log ( Uniform(this->stream)));
    if (statFlag  == 2) CollectStat(sample);
    return sample;
    };

void Exponential::PrintProfile() {
        printf("\n Exponential(%f) generator (%f, %f)",mean,mean,std);
    };

HyperExp::HyperExp(float mean, float variance, int stream)
{
        float   cx2;

        this->stream = stream;
        this->mean = mean;
        this->std = sqrt(variance);
        cx2  = variance/ (mean * mean);
        alpha = 1/(1 + cx2);
        mu1  =  mean * (1/( 1 + sqrt( 1/(alpha - 1.0)*(cx2 -1)/2)));
        mu2  =  mean * (1/( 1 - sqrt((1/(1 - alpha) - 1) * (cx2 -1)/2)));
}

HyperExp::HyperExp(float alpha, float mu1, float mu2, int stream)
{
        this->stream = stream;
        this->alpha = alpha;
        this->mu1 = mu1;
        this->mu2 = mu2;

        mean = alpha*mu1 + (1 - alpha)*mu2;
        std = sqrt(alpha*mu1*mu1 + (1 - alpha)*mu2*mu2);
        std = sqrt(2*(alpha*mu1*mu1 + (1 - alpha)*mu2*mu2)
                - (alpha*mu1+(1-alpha)*mu2)*(alpha*mu1+(1-alpha)*mu2));
}

float   HyperExp::Draw()
{
        float   sample;

        if (Uniform(stream) <= alpha)
            sample = -log(Uniform(stream))*mu1;
        else
            sample = -log(Uniform(stream))*mu2;
        CollectStat(sample);
        return sample;
}

void    HyperExp::PrintProfile() {
        printf("\nHyperexponential[%f, %f,%f] generator (%f, %f)",
            alpha,mu1,mu2,mean,std);
}

Normal::Normal(float mean, float sigma, int stream) {
        this->stream = stream;
        this->mean = mean;
        this->std = sigma;
        statFlag  = 2;
    };

float  Normal::Draw() {
    float       v1, v2;
    float       sample;

    do {
        v1 = - log(Uniform(this->stream));
        v2 = - log(Uniform(this->stream));
        } while( (2.0* v1) < ((v2 - 1.0)*(v2 - 1.0)));

    if (Uniform(this->stream) < 0.5)
        v2 *= -1;
    sample = std * v2 + mean;

    if (statFlag  == 2) CollectStat(sample);
    return sample;
    };
void  Normal::PrintProfile() {
        printf("\n Normal[%f,%f] generator (%f, %f)",mean,std,mean,std);
    };

LogNormal::LogNormal(float mean, float sigma, int stream) {
        this->stream = stream;
        this->mean = mean;
        this->std = sigma;
        statFlag  = 2;
    };

float  LogNormal::Draw() {
    float       v1, v2;
    float       sample;

    do {
        v1 = - log(Uniform(this->stream));
        v2 = - log(Uniform(this->stream));
        } while( (2.0* v1) < ((v2 - 1.0)*(v2 - 1.0)));

    if (Uniform(this->stream) < 0.5)
        v2 *= -1;
    sample = exp(std * v2 + mean);

    if (statFlag  == 2) CollectStat(sample);
    return sample;
    };

void  LogNormal::PrintProfile() {
        printf("\n LogNormal[%f,%f] generator (%f, %f)",mean,std,mean,std);
    };

Gamma::Gamma(float alpha, float beta, int stream) {
        this->alpha = alpha;
        this->beta  = beta;
        mean        = alpha/beta;
        std         = sqrt(alpha)/beta;
        statFlag  = 2;
}

float   GammaDist(int stream, float alpha, float beta) {
        float a, v1, v2;

        a = alpha - 1.0;
        do {
            v1 = -log( Uniform(stream));
            v2 = -log( Uniform(stream));
        } while ( v2 < a* (v1 - log(v1) - 1.0));

        return alpha/beta*v1;
    };

float   Gamma::Draw() {
        float sample = GammaDist(this->stream,alpha,beta);
        if (statFlag  == 2) CollectStat(sample);
        return sample;
    };

void Gamma::PrintProfile()
{
        printf("\n Gamma [%f, %f]  generator (%f, %f)",alpha,beta,mean,std);
}

Beta::Beta(float alpha, float beta, int stream)
{
        this->alpha  = alpha;
        this->beta  = beta;
        statFlag     = 2;
        mean = alpha/(alpha + beta);
        std = sqrt((alpha*beta / (alpha + beta + 1))/((alpha + beta)*(alpha + beta)));
}

float Beta::Draw()
{
        float   y1, y2, sample;

        y1 = GammaDist(stream,alpha,beta);
        y2 = GammaDist(stream,alpha,beta);
        sample  = y1/(y1+y2);
        if (statFlag  == 2) CollectStat(sample);
        return sample;
}

void Beta::PrintProfile() {
    printf("\n Beta[%f, %f] generator (%f, %f)",alpha,beta,mean,std);
    };

GeneralCDF::GeneralCDF(int n, float* v, float* p , int stream) {
    int i;

    this->n = n;
    this->v = new float[n];
    this->p = new float[n];
    q       = new float[n];

    this->v[0] = v[0];
    this->p[0] = p[0];
    this->q[0] = p[0];
    for(i=1; i < n; i++){
        this->v[i] = v[i];
        this->p[i] = p[i];
        this->q[i] = p[i] + p[i-1];
        }
    }

float   GeneralCDF::Draw() {
        float   sample, result;

        sample = Uniform(stream);
        result = v[binary_search(0, n-1, q, sample)];
        CollectStat(result);
        return result;
}

GeneralCDF::~GeneralCDF()
{
        delete this->v;
        delete this->p;
        delete this->q;
}

void    GeneralCDF::PrintProfile()
{
        printf("\nGeneral CDF generator (%f, %f)",mean,std);
}

int     binary_search(int low, int high, float* q, float val)
{
        int     mid;

        while (low <= high) {
                mid = (low + high)/2;
                if (q[mid] == val) return mid;
                else if (q[mid] < val) {
                        low = mid + 1;
                        return binary_search(low, high, q, val);
                }
                else {
                        high = mid - 1;
                        return binary_search(low, high, q, val);
                }
        }
        return low;
}