blob: 5c9370b021f9cb9938b431bf601310b1997399e9 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
|
/*
* Normaldistribution.cpp
*
* Created on: 20.05.2011
* Author: sven
*/
#include "Normaldistribution.h"
#include <cmath>
#include <cstdlib>
#include <iostream>
Normaldistribution::Normaldistribution() {
}
Normaldistribution::~Normaldistribution() {
}
void Normaldistribution::calcRandomVector(const CholeskyMatrix& D, const std::vector<double>& mu , std::vector<double>& y)
{
std::vector<double> u(nenEndOfCompanies);
calcStandardNormalDistributedRandomVector(nenEndOfCompanies,u);
for (int i=0;i<nenEndOfCompanies;i++)
{
double sum = 0.0;
for (int j=0; j< nenEndOfCompanies ; j++)
{
sum += (D.at(i).at(j) * u.at(j));
}
y.at(i) = sum + mu.at(i);
// std::cout << y.at(i) << std::endl;
}
// std::cout << std::endl;
}
void Normaldistribution::calcStandardNormalDistributedRandomVector(const int size, std::vector<double>& u)
{
double x1,x2;
double v1,v2;
double u1,u2;
double s;
int counter = 0;
while (counter < size)
{
do
{
x1 = ((double) rand()) / (RAND_MAX + 1.0);
x2 = ((double) rand()) / (RAND_MAX + 1.0);
v1 = (2.0 * x1) - 1.0;
v2 = (2.0 * x2) - 1.0;
s = (v1*v1) + (v2*v2);
} while (s >= 1.0);
u1 = v1 * std::sqrt( - (2.0 / s) * std::log(s) );
u2 = v2 * std::sqrt( - (2.0 / s) * std::log(s) );
if (counter < size) {
u.at(counter) = u1;
counter++;
}
if (counter < size) {
u.at(counter) = u2;
counter++;
}
}
}
|
