Free Web Site - Free Web Space and Site Hosting - Web Hosting - Internet Store and Ecommerce Solution Provider - High Speed Internet
Search the Web

Generación de la Distribución Normal Acumulada

Se sabe que la función que genera la conocida tabla de probabilidades es :


obviamente después de hacer una integración numérica entre ciertos valores.

Lo primero que podemos hacer sin mayor dificultad es definir una función que retorne lo anterior

double func_normal(double t)
{
    return (double) exp(-pow(t,2)/ 2) / sqrt(2 * 3.14159260);
}

El valor de PI se define como 3.14159260 , para obtener precision.

Proceso de Integración

Para el proceso de integración, necesitamos un metodo numérico que se pueda implementar facilmente. Para esto podemos usar el metodo de simpson, que calcula el valor aproximado del area bajo la curva, bajo ciertos criterios.

Para empezar vamos a definir una función que aplica la siguiente formula :

Ap = Yo + 4 * F(Xo + h) + F(Xo + 2h)

Para no realizar dos veces la evaluación de la función y obtener un manejo especial de la variable x, se opto por lo siguiente :

El valor de Yo y Xo van a ser conocidos. Despues de evaluar la función, estos parametros se modifican segun :

Yo : sale con el valor Yo = F(Xo + 2h)
Xo : sale con el valor Xo = Xo + 2h

Para mayor flexibilidad, la función a evaluar se pasará como parámetro, a través de u puntero.

double suma_parcial(double *x , double *y , double h , double(*f)(double))
{
   double s = *y;
   s += 4 * f(*x + h);
   *x += 2 * h;
   *y = f(*x);
   s += *y;
   return s;
}

Ahora implementamos una función que evalúe lo anterior, en todos los intervalos.

double integral(double a , double b ,int n , double (*f)(double))
{
   double s = 0.0;
   double y = f(a);
   if (n % 2) n++;
   double h = (b - a) / n;
   int cont = 0;
   do
   {
       s += suma_parcial(&a,&y,h,f);
       cont += 2;
   } while(cont < n);
   s *= h / 3;
   return s;
}

Después de esto, solo nos resta integrar desde t = -3.6 hasta z con incrementos de 0.01.
Luego de varios ajustes, se vio que la probabilidad que z fuera menor que -3.6 era de 0.000159 para que los números generados por estas funciones, fueran idénticos a los de las tablas estadisticas.

double normal_acum(double z)
{
   double lim_inf = -3.60;
   double acum = 0.000159;
   double inc = 0.01;
   do
   {
      acum += integral(lim_inf,lim_inf + inc, 4 , func_normal);
      lim_inf += inc;
   } while(lim_inf < z - inc);
   
   return acum; 
}

Por simetría, se pueden evitar los calculos para todas las probabilidades mayores de 0.50 , las cuales son asociadas para todo z menor que cero, con la siguiente función :

double dist_normal(double z)
{
    return (z <= 0.0)? normal_acum(z) : 1.0 - normal_acum(-z);
}

*nota : Al salir los datos a pantalla, hay que redondear el resultado de la función a sólo cuatro decimales, tal como el formato de las tablas estadísticas ( %.4lf ).

Después de estos calculos, lo unico que falta es la función inversa : "Dada cierta probabilidad, encontrar el z asociado ". Para ello modificamos un poco la funcion normal_acum...

double func_z(double prob)
{
   double lim_inf = -3.60;
   double acum = 0.000159;
   double inc = 0.01;
   do
   {
      acum += integral(lim_inf,lim_inf + inc, 4 , func_normal);
      lim_inf += inc;
   } while(acum < prob + 0.000159);
   
   return lim_inf - inc; 
}

Igualmente, por simetría ...

double obtener_z(double prob)
{
   double z_aux = (prob <= 0.5)? func_z(prob) : -func_z(1.0 - prob);
   return z_aux; 
}

Con respecto a la función obtener_z sucede algo especial. Ya que en las tablas estadisticas se da por ejemplo que :

P(z < -3.18) = P(z < -3.19) = 0.0007;
P(z < -3.20) = P(z < -3.21) = 0.0007;

sobretodo para los valores de z cercanos a los limites (-3.5 y 3.5) , la evaluación de la función obtener_z retorna el menor valor de z.

Para poder comprobrar la exactitud del proceso de calculo, se realizo un programa de prueba. El programa consiste en un numero aleatorio de valores de z, almacenados en un arreglo, para ser evaluados en la función dist_normal y obtener_z.

La salida del programa es la siguiente :

Codigo fuente :

#include<stdio.h>
#include<stdlib.h>
#include "distnormal.h"

double *arr_aleat(int);

int main()
{
    
  int opcion = 1;
  int n;
  double p;
  double *arr;
    
  printf("\nIngrese numero de pruebas (0 .. 36) : ");
  scanf("%d",&n);
  printf("\n");
  arr = arr_aleat(n);

  printf("\tz\tP(t < z)\tf(p)\t***");
  printf("\tz\tP(t < z)\tf(p)");
  printf("\n\n");
        
  for(int i = 0; i < n; i++)
  {
    p = dist_normal(arr[i]);
    
    if (arr[i] > 0.0)
    {
      printf("\t+%4.2lf\t%8.4lf",arr[i],p);
      printf("\t+%4.2lf",obtener_z(p));
    }
    else
    {
      printf("\t%4.2lf\t%8.4lf",arr[i],p);
      printf("\t%4.2lf",obtener_z(p));
    }
            
    if (i % 2 != 0) printf("\n");
    else printf("\t***");
  }
  
  free(arr);        
  printf("\n");
    
  return 0;
}

double *arr_aleat(int n)
{
  double lim = -260;
  int cont = 0;
  int digit_aleat;
    
  double *arr = (double *) calloc(n,sizeof(double));

  while(cont < n)
  {
    digit_aleat = rand() % 10;
    arr[cont] = (double) (lim + digit_aleat) / 100;
    lim += 15;
    cont++;
  }
    
  return arr;
}