用C语言实现瑞利分布，莱斯分布，高斯分布的分布函数

# include "stdio.h"

# include "math.h"
# include "stdlib.h"
# include "math.h"
# include "dos.h"
# define MAX_N 3000  /*这个值为N可以定义的最大长度*/
# define N 100       /*产生随机序列的点数，注意不要大于MAX_N*/

/*产生均匀分布的随机变量*/
void randa(float *x,int num); 

/*产生瑞利分布的随机变量*/
void randr(float *x,int num); 

/*产生标准高斯分布的随机变量*/
void randn(float *x,int num);

/*产生莱斯分布的随机变量*/
void randl(float *x, float a, float b, int num);

void fshow(char *name,float *x,int num);


main()
{
   float x[N];
   int i;
/*
   randa(&x,N);
   randr(&x,N);
   randl(&x,10,10,N);

*/
   randn(&x,N);
/*此时x[N]就是所需要的高斯分布的序列*/

/*显示该序列*/
   fshow("x",&x,N);

   getch();

}

void randa(float *x,int num)
{
  int i;
  struct time stime;
  unsigned seed;
  gettime(&stime);
  seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
  srand(seed);
  for(i=0;i<num;i++)
  {
     x[i]=rand();
     x[i]=x[i]/32768;
  }
}

void randr(float *x,int num)
{
  float x1[MAX_N];
  int i;
  struct time stime;
  unsigned seed;
  gettime(&stime);
  seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
  srand(seed);
  for(i=0;i<num;i++)
  {
     x1[i]=rand();
     x[i]=x1[i]/32768;
     x[i]=sqrt(-2*log(x[i]));
  }

}
void randn(float *x,int num)
{
  float x1[MAX_N],x2[MAX_N];
  int i;
  struct time stime;
  unsigned seed;
  gettime(&stime);
  seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
  srand(seed);
  for(i=0;i<num;i++)
  {
     x1[i]=rand();
     x2[i]=rand();
     x1[i]=x1[i]/32768;
     x2[i]=x2[i]/32768;
     x[i]=sqrt(-2*log(x1[i]))*cos(x2[i]*M_PI);
  }

}
void randl(float *x, float a, float b, int num)
{
  float x1[MAX_N],x2[MAX_N];
  float temp[MAX_N];
  int i;
  struct time stime;
  unsigned seed;
  gettime(&stime);
  seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
  srand(seed);
  for(i=0;i<num;i++)
  {
     x1[i]=rand();
     x2[i]=rand();
     x1[i]=x1[i]/32768;
     x2[i]=x2[i]/32768;
     temp[i]=sqrt(-2*log(x1[i]))*cos(x2[i]*M_PI);
     x2[i]=sqrt(-2*log(x1[i]))*sin(x2[i]*M_PI);
     x1[i]=temp[i];
     x[i]=sqrt((a+x1[i])*(a+x1[i])+(b+x2[i])*(b+x2[i]));
  }

}


void fshow(char *name,float *x,int num)
{
  int i,sign,L;
  float temp;
  printf("\n");
  printf(name);
  printf(" = ");
  L=6;
  /*按照每行6个数据的格式显示*/
  for(i=0;i<num;i++)
  {
    temp=i/L;
    sign=temp;
    if((i-sign*L)==0) printf("\n");
    if(x[i]>0) printf(" %f   ",x[i]);
    else       printf("%f   ",x[i]);
  }
}


程序 2
以下程序加入了图形显示的效果，因此更加直观，你可以参考一下。


/* 作者 Leo_nanjing
   时间 2008.5.10
   功能 生成各种分布的随机变量,并显示
*/


# include "stdio.h"
# include "math.h"
# include "graphics.h"
# include "math.h"
# include "dos.h"
# define MAX_N 3000
# define N 1000


void randa(float *x,int num);
void randr(float *x,int num);
void randn(float *x,int num);
void randl(float *x, float a, float b, int num);
void fshow(char *name,float *x,int num);


/*用于图形显示的部分*/
void init_graphic(unsigned color);
void plotxy(float *x, float *y, int num,int mode);
void plot(float *y,int num, int mode);
float max(float *x, int num);
float min(float *x, int num);
/*画出该随机序列的分布函数曲线*/
void plotpdf(float *x,int num,int part,int mode);


main()
{
   float x[N];
   int i;
   randn(&x,N);
   fshow("x",&x,N);
   getch();


   /*以下为图形显示部分*/
   init_graphic(0);
   /*显示随机序列*/
   plot(&x,N,1);
   getch();
   /*显示其分布函数*/
   plotpdf(&x,N,20,0);
   getch();
}

void randa(float *x,int num)
{
  int i;
  struct time stime;
  unsigned seed;
  gettime(&stime);
  seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
  srand(seed);
  for(i=0;i<num;i++)
  {
     x[i]=rand();
     x[i]=x[i]/32768;
  }
}

void randr(float *x,int num)
{
  float x1[MAX_N];
  int i;
  struct time stime;
  unsigned seed;
  gettime(&stime);
  seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
  srand(seed);
  for(i=0;i<num;i++)
  {
     x1[i]=rand();
     x[i]=x1[i]/32768;
     x[i]=sqrt(-2*log(x[i]));
  }

}
void randn(float *x,int num)
{
  float x1[MAX_N],x2[MAX_N];
  int i;
  struct time stime;
  unsigned seed;
  gettime(&stime);
  seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
  srand(seed);
  for(i=0;i<num;i++)
  {
     x1[i]=rand();
     x2[i]=rand();
     x1[i]=x1[i]/32768;
     x2[i]=x2[i]/32768;
     x[i]=sqrt(-2*log(x1[i]))*cos(x2[i]*M_PI);
  }

}
void randl(float *x, float a, float b, int num)
{
  float x1[MAX_N],x2[MAX_N];
  float temp[MAX_N];
  int i;
  struct time stime;
  unsigned seed;
  gettime(&stime);
  seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
  srand(seed);
  for(i=0;i<num;i++)
  {
     x1[i]=rand();
     x2[i]=rand();
     x1[i]=x1[i]/32768;
     x2[i]=x2[i]/32768;
     temp[i]=sqrt(-2*log(x1[i]))*cos(x2[i]*M_PI);
     x2[i]=sqrt(-2*log(x1[i]))*sin(x2[i]*M_PI);
     x1[i]=temp[i];
     x[i]=sqrt((a+x1[i])*(a+x1[i])+(b+x2[i])*(b+x2[i]));
  }

}


void fshow(char *name,float *x,int num)
{
  int i,sign,L;
  float temp;
  printf("\n");
  printf(name);
  printf(" = ");
  L=6;
  for(i=0;i<num;i++)
  {
    temp=i/L;
    sign=temp;
    if((i-sign*L)==0) printf("\n");
    if(x[i]>0) printf(" %f   ",x[i]);
    else       printf("%f   ",x[i]);
  }
}




/*以下为图形显示的函数*/
void init_graphic(unsigned color)
{
   int graphicdriver,graphicmode;
   graphicdriver=DETECT;
   graphicmode=1;
   initgraph(&graphicdriver,&graphicmode,"E:\\turboc2\\");
   setbkcolor(color);
}

void plotxy(float *x, float*y, int num,int mode)
{
   int i;
   float max_x,max_y,min_x,min_y;
   float x0,y0,x1,y1;
   clrscr(0);
   cleardevice();
   setbkcolor(0);
   max_x=max(x,num);
   max_y=max(y,num);
   min_x=min(x,num);
   min_y=min(y,num);
   setlinestyle(0,2,1);
   line(65,35,65,445);
   line(65,445,575,445);
   setlinestyle(3,0,1);
   line(65,35,575,35);
   line(575,35,575,445);
   setlinestyle(0,2,1);
   if(max_x==min_x)
     x0=320;
   else
     x0=(x[0]-min_x)*500/(max_x-min_x)+70;
   if(max_y==min_y)
     y0=240;
   else
     y0=480-((y[0]-min_y)*400/(max_y-min_y)+40);
   if(mode==0) circle(x0,y0,2);
   for(i=1;i<num;i++)
   {
      if(max_x==min_x)
        x1=320;
      else
        x1=(x[i]-min_x)*500/(max_x-min_x)+70;
      if(max_y==min_y)
        y1=240;
      else
        y1=480-((y[i]-min_y)*400/(max_y-min_y)+40);
      if(mode==0) circle(x1,y1,2);
      line(x0,y0,x1,y1);
      x0=x1;y0=y1;
   }
   printf("\n\n");
   printf("%f",max_y);
   printf("\n\n\n\n\n\n\n\n\n\n");
   printf("\n\n\n");
   printf("%f",(max_y+min_y)/2);
   printf("\n\n\n\n\n\n\n\n\n\n");
   printf("\n\n");
   printf("%f",min_y);
   printf("\n       %f",min_x);
   printf("                       ");
   printf("%f",(max_x+min_x)/2);
   printf("                      ");
   printf("%f",max_x);

}

void plot(float*y, int num,int mode)
{
   int i;
   float max_x,max_y,min_x,min_y;
   float x0,y0,x1,y1;
   float x[MAX_N];
   clrscr(0);
   cleardevice();
   setbkcolor(0);
   for(i=0;i<num;i++) x[i]=i+1;
   max_x=max(x,num);
   max_y=max(y,num);
   min_x=min(x,num);
   min_y=min(y,num);
   setlinestyle(0,2,1);
   line(65,35,65,445);
   line(65,445,575,445);
   setlinestyle(3,0,1);
   line(65,35,575,35);
   line(575,35,575,445);
   setlinestyle(0,2,1);
   if(max_x==min_x)
     x0=320;
   else
     x0=(x[0]-min_x)*500/(max_x-min_x)+70;
   if(max_y==min_y)
     y0=240;
   else
     y0=480-((y[0]-min_y)*400/(max_y-min_y)+40);
   if(mode==0) circle(x0,y0,2);
   for(i=1;i<num;i++)
   {
      if(max_x==min_x)
        x1=320;
      else
        x1=(x[i]-min_x)*500/(max_x-min_x)+70;
      if(max_y==min_y)
        y1=240;
      else
        y1=480-((y[i]-min_y)*400/(max_y-min_y)+40);
      if(mode==0) circle(x1,y1,2);
      line(x0,y0,x1,y1);
      x0=x1;y0=y1;
   }
   printf("\n\n");
   printf("%f",max_y);
   printf("\n\n\n\n\n\n\n\n\n\n");
   printf("\n\n\n");
   printf("%f",(max_y+min_y)/2);
   printf("\n\n\n\n\n\n\n\n\n\n");
   printf("\n\n");
   printf("%f",min_y);
   printf("\n       %f",min_x);
   printf("                       ");
   printf("%f",(max_x+min_x)/2);
   printf("                      ");
   printf("%f",max_x);

}

void plotpdf(float *x,int num,int part,int mode)
{
   int i,j;
   float max_x,min_x,round,deltax,up,down,sum;
   float xl[MAX_N],yl[MAX_N];
   sum=0;
   max_x=max(x,num);
   min_x=min(x,num);
   round=max_x-min_x;
   deltax=round/part;
   xl[0]=min_x;
   for(i=1;i<=part;i++)
   {
      xl[i]=min_x+deltax*i;
      yl[i-1]=0;
      up=xl[i];
      down=xl[i-1];
      for(j=0;j<num;j++)
      {
  if((x[j]<up) && (x[j]>=down)) yl[i-1]=yl[i-1]+1;

      }
      yl[i-1]=yl[i-1]/num/deltax;
   }
   for(i=0;i<part;i++) sum=sum+yl[i];
   plotxy(&xl,&yl,part,mode);

}


float max(float *x, int num)
{
   int i;
   float max;
   max=x[0];
   for(i=1;i<num;i++)
   {
      if(x[i]>max) max=x[i];
   }
   return max;
}
float min(float *x, int num)
{
   int i;
   float min;
   min=x[0];
   for(i=1;i<num;i++)
   {
      if(x[i]<min) min=x[i];
   }
   return min;

}