sábado, 21 de mayo de 2016

Métodos de generación de números aleatorios en lenguaje C

Aquí muestro la implementación, en lenguaje C, de varios tipos de generadores aleatorios mostrados en el siguiente paper

// Method-1: Sum of Uniform Random Variables
double metodo1_Sum_Of_Uniform_Random_Variables() {
    double X[12];
    double z= 0.0;

    /* Generar 12 numeros con distribucion U(0,1)  */
    for (int i = 0; i < 12; i++ ){
        X[i] = ( (double)rand() / (double)RAND_MAX );
        z+= X[i];
    }

    return z-6;
}

// Box-Muller Method
double metodo2_Box_Muller()
{
    //Generar numero aleatorio uniforme U(0,1)
    double u1 = (double)rand() / (double)RAND_MAX;
    double u2 = (double)rand() / (double)RAND_MAX;

    double z1 = sqrt(-2.0*log(u1))*cos(2*PI*u2); // usa la box-muller transformation
    double z2 = sqrt(-2.0*log(u1))*sin(2*PI*u2); // usa la box-muller transformation

    // Aqui promediamos el valor de las 2 Z para evitar tener que reiterar
    // Este paso es opcional si se retorna el valor de z1 o z2
    double zt= (z1+z2)/2;
    return zt;
}

// Polar Box Muller
double metodo3_Polar_Box_Muller() {

    double S, Z1, Z2, U1, U2, V1,V2;

    do
    {
        // #1 Genera dos números entre -1 y 1
        U1 = (double)rand() / RAND_MAX;
        U2 = (double)rand() / RAND_MAX;

        // #2 V1(-1,1)
        V1 = 2. * U1 - 1.; // V1(-1,1)
        V2 = 2. * U2 - 1.; // V2(-1,1)
        S = V1 * V1 + V2 * V2;

        // #3
    } while(S > 1);

    // #4
    Z1 = sqrt (-2. * log(S) / S)*V1;
    Z2 = sqrt (-2. * log(S) / S)*V2;

    return (Z1+Z2)/2;
}

// Inversion Method
double metodo4_Inversion_Method()
{
    // #1
    double U=(double)rand() / RAND_MAX;

    if(U == 0.5)
        return 0.0;

    // #2 3
    double w=sqrt(-2*log(U < 0.5 ? U : 1-U));

    double a[]={ 2.515517, 0.802853, 0.010328 };
    double b[]={ 1.000000, 1.432788, 0.189269, 0.001308 };

    double Z = -w +      ( (a[0] + a[1]*w + a[2] * w*w) /
            (b[0] + b[1]*w + b[2] * w*w + b[3]*w*w*w));

    // #3
    return U > 0.5 ? -Z : Z;
}

// Acceptance Rejection Method
double metodo5_Acceptance_Rejection_Method()
{
    double Y1=0.0;
    double Y2=0.0;

    // #1 Genera dos Ui ~ U(0,1).
    do{
        double U1=(double)rand()/RAND_MAX;
        double U2=(double)rand()/RAND_MAX;

        Y1= -log(U1);
        Y2= -log(U2);

        // #2
    }while(Y2 > ((Y1-1)*(Y1-1))/2.0 );

    double Z=Y1;

    // #3
    double U3=(double)rand()/RAND_MAX;

    return U3<=0.5 ? fabs(Z) : -fabs(Z);
}

// Using Generalized Exponential Distribution
double metodo6_Using_Generalized_Exponential_Distribution() {

    //Generar numero aleatorio uniforme U(0,1)
    double u = (double)rand() / (double)RAND_MAX;
    double x = -log(1-pow(u,0.0775));

    double z=(log(x)-1.0821)/0.3807;

    //retornar valor de z
    return z;
}

// Bol’shev Formula
double metodo7_Bolshev_Formula() {

    // Array de 12 numeros entre 0 y 1
    double U[12];

    // acumulador
    double X=0.0;
    /* Generar 5 numeros con distribucion U(0,1)  */
    for (int j = 0; j < 12; j++ ){
        U[j] = ( (double)rand() / (double)RAND_MAX );

        X+=sqrt(3)*( 2 * U[j] -1 );
    }
    X*=(1/sqrt(5));
    double z= X-0.01*(3*X-pow(X,3));
    return z;
}

// Inversion Method
double metodo8_Inversion_Method() {
    //Generar numero aleatorio uniforme U(0,1)
    double u = (double)rand() / (double)RAND_MAX;

    //retornar valor de z
    return (-log(1)+ log(u)+ 1);
}

// Proposed Method
double metodo9_Proposed_Method(){
    double U  = (double)rand() / (double)RAND_MAX;
    double X1 = tanh(-31.35694+28.77154*U);
    double X2 = tanh(-2.57136-31.16364*U);
    double X3 = tanh(3.94963-1.66888*U);
    double X4 = tanh(2.31229+1.84289*U);
    double Z=0.46615+90.72192*X1-89.36967*X2-96.55499*X3+97.36346*X4;
    return Z;
}

double calcValorReal(double M, double a, double Z)
{
    //// X=Z*a+M
    return Z*a+M;
}

No hay comentarios:

Publicar un comentario