#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <windows.h>
#define PI 3.1415926
double guassrnd (double mean, double sigma)
{
int i;
double var;
double random;
random = 0;
for (i = 0; i < 12; i++)
random += (rand() / (RAND_MAX + 1.0));
var = mean + sigma * (random - 6);
return var;
}
double guassrnd2()
{
static double V1, V2, S;
static int phase = 0;
double X;
if ( phase == 0 ) {
do {
double U1 = (double)rand() / RAND_MAX;
double U2 = (double)rand() / RAND_MAX;
V1 = 2 * U1 - 1;
V2 = 2 * U2 - 1;
S = V1 * V1 + V2 * V2;
} while(S >= 1 || S == 0);
X = V1 * sqrt(-2 * log(S) / S);
} else
X = V2 * sqrt(-2 * log(S) / S);
phase = 1 - phase;
return X;
}
double uniform(double a,double b,long int *seed)
{
// double a,b;
//long int *seed;
double t;
*seed=2045*(*seed)+1;
*seed=*seed-(*seed/1048576)*1048576;
t=(*seed)/1048576.0;
t=a+(b-a)*t;
return(t);
}
double guassrnd3(double mean, double sigma, long int *s)
{
//long int *s;
int i;
double x,y;
for(x= 0,i=0;i<12;i++)
x += uniform(0.0,1.0,s);
x = x - 6;
y = mean + x*sigma;
return(y);
}
double distribute (double sigma, double mean, double x)
{
double s, m;
double s1, s2;
double fun;
s = sigma;
m = mean;
s1 = sigma * sqrt (2 * PI);
s2 = sigma * sigma;
fun = 1 / s1 * exp (-(x - mean) * (x - mean) / 2 * s2);
return fun;
}
void main ()
{
FILE *fpt;
int i, test_num;
double var, fun, mean, sigma; //mean-均值，sigma-方差
long int s;
test_num = 10000;
mean = 0.5;
sigma = 0.1;
srand(time(NULL));
fpt = fopen ("test.txt", "w");
for (i = 0; i < test_num; i++)
{
//var = guassrnd (mean, sigma);
//var = guassrnd2();
s = rand();
var = guassrnd3(mean, sigma, &s);
fun = distribute (sigma, mean, var);
fprintf (fpt, "%lf\t%lf\n", var, fun);
}
}
- 1
- 2
- 3
前往页