注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Mr.Right

不顾一切的去想,于是我们有了梦想。脚踏实地的去做,于是梦想成了现实。

 
 
 

日志

 
 
关于我

人生一年又一年,只要每年都有所积累,有所成长,都有那么一次自己认为满意的花开时刻就好。即使一时不顺,也要敞开胸怀。生命的荣枯并不是简单的重复,一时的得失不是成败的尺度。花开不是荣耀,而是一个美丽的结束,花谢也不是耻辱,而是一个低调的开始。

网易考拉推荐

C语言正态分布随机数,高斯白噪声demo  

2014-09-12 01:50:03|  分类: 编程 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

// Generating random numbers from Normal distribution in C
// Before plotting make sure to adjust the MAX_WIDTH macro so that the histogram bars does not wrap lines.

// demo :  mu = 0, sigma = 1, samples = 50000, bins = 70 and the min and max were not provided
// in cmd line: randn 0 1 50000 70
//                : randn 10 4 50000 70
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <limits.h>
 
#define MAX_WIDTH 80
 
 // https://en.wikipedia.org/wiki/Marsaglia_polar_method
double randn (double mu, double sigma);
 
/* Arguments: ./a.out mu sigma samples bins [min] [max] */
 
int main (int argc, char *argv[])
{
  int samples, bins;
  double *range, mu = 0.0, std = 1.0, min, max, current_random_no;
  int i, j, *bin_count, count = 0, bin, max_bin_count = INT_MIN, flag = 0;
 
  if (argc < 5)
    {
      printf ("Usage: %s mu sigma samples bins [min] [max]\n", argv[0]);
      return 0;
    }
 
  srand (time (NULL));
 
  mu = atof (argv[1]);
  std = atof (argv[2]);
  samples = atoi (argv[3]);
  bins = atoi (argv[4]);
 
  /* Allocate memory */
  bin_count = calloc (sizeof (int), bins);
  range = malloc (sizeof (double) * (bins + 1));
 
  /* Automatically set a min and max range */
  min = -(abs (mu) + 5 * std);
  max = abs (mu) + 5 * std;
 
  /* If min and max is specified, use them */
  if (argc >= 6)
    {
      min = atof (argv[5]);
    }
  if (argc >= 7)
    {
      max = atof (argv[6]);
    }
 
  /* Generate the bin ranges */
  range[0] = min;
  for (i = 1; i <= bins; i++)
    {
      range[i] = range[0] + (max - min) * (1.0 / bins) * i;
    }
 
  for (count = 0; count < samples; count++)
    {
      /* Generate random numbers from a distribution, Normal for this */
      current_random_no = randn (mu, std);
 
      /* Check which range the current sample falls */
      for (i = 0, bin = -1; i < bins; i++)
        {
          if ((current_random_no >= range[i])
              && (current_random_no < range[i + 1]))
            {
              bin = i;
              break;
            }
        }
      /* In case we have the number exactly equalto range[bins] */
      if (current_random_no == range[i])
        {
          bin = i - 1;
        }
      if ((bin <= bins) && (bin >= 0))
        {
          bin_count[bin]++;
        }
      /* Not all the random numbers were generated within the [min,max] range. */
      else
        {
          flag = 1;
        }
    }
 
  /* Find the max value of the bin counters. This is used to scale the histogram */
  for (i = 0; i < bins; i++)
    {
      if (bin_count[i] > max_bin_count)
        {
          max_bin_count = bin_count[i];
        }
    }
 
 
  /* Print histogram and ranges */
  printf ("[bin_low, bin_high), count\n");
  for (i = 0; i < bins; i++)
    {
      printf ("[%+7.1f,%+7.1f), %6d:", range[i], range[i + 1], bin_count[i]);
      for (j = 0;
           j < (bin_count[i] / (double) max_bin_count) * (MAX_WIDTH / std);
           j++)
        {
          printf ("*");
        }
      printf ("\n");
 
    }
 
  /* Not all random numbers were in the range of the histogram */
  if (flag == 1)
    {
      printf
        ("\nWARNING: Random numbers generated outside range (%f, %f). These will not be shown in the above histogram.\n",
         min, max);
    }
 
 
  free (bin_count);
  free (range);
 
  printf ("\n");
  return 0;
}
 
 
double
randn (double mu, double sigma)
{
  double U1, U2, W, mult;
  static double X1, X2;
  static int call = 0;
 
  if (call)
    {
      call = !call;
      return (mu + sigma * (double) X2);
    }
 
  do
    {
      U1 = -1 + ((double) rand () / RAND_MAX) * 2;
      U2 = -1 + ((double) rand () / RAND_MAX) * 2;
      W = pow (U1, 2) + pow (U2, 2);
    }
  while (W >= 1 || W == 0);
 
  mult = sqrt ((-2 * log (W)) / W);
  X1 = U1 * mult;
  X2 = U2 * mult;
 
  call = !call;
 
  return (mu + sigma * (double) X1);
}

  评论这张
 
阅读(586)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016