|
Lines 1033-1039
double NormalVariableImpl::GetValue()
|
Link Here
|
|---|
|
|---|
| 1033 |
} |
1033 |
} |
| 1034 |
while(1) |
1034 |
while(1) |
| 1035 |
{ // See Simulation Modeling and Analysis p. 466 (Averill Law) |
1035 |
{ // See Simulation Modeling and Analysis p. 466 (Averill Law) |
| 1036 |
// for algorithm |
1036 |
// for algorithm; basically a Box-Muller transform: |
|
|
1037 |
// http://en.wikipedia.org/wiki/Box-Muller_transform |
| 1037 |
double u1 = m_generator->RandU01(); |
1038 |
double u1 = m_generator->RandU01(); |
| 1038 |
double u2 = m_generator->RandU01();; |
1039 |
double u2 = m_generator->RandU01();; |
| 1039 |
double v1 = 2 * u1 - 1; |
1040 |
double v1 = 2 * u1 - 1; |
|
Lines 1082-1088
double NormalVariableImpl::GetSingleValu
|
Link Here
|
|---|
|
|---|
| 1082 |
} |
1083 |
} |
| 1083 |
while(1) |
1084 |
while(1) |
| 1084 |
{ // See Simulation Modeling and Analysis p. 466 (Averill Law) |
1085 |
{ // See Simulation Modeling and Analysis p. 466 (Averill Law) |
| 1085 |
// for algorithm |
1086 |
// for algorithm; basically a Box-Muller transform: |
|
|
1087 |
// http://en.wikipedia.org/wiki/Box-Muller_transform |
| 1086 |
double u1 = m_static_generator->RandU01(); |
1088 |
double u1 = m_static_generator->RandU01(); |
| 1087 |
double u2 = m_static_generator->RandU01();; |
1089 |
double u2 = m_static_generator->RandU01();; |
| 1088 |
double v1 = 2 * u1 - 1; |
1090 |
double v1 = 2 * u1 - 1; |
|
Lines 1092-1102
double NormalVariableImpl::GetSingleValu
|
Link Here
|
|---|
|
|---|
| 1092 |
{ // Got good pair |
1094 |
{ // Got good pair |
| 1093 |
double y = sqrt((-2 * log(w))/w); |
1095 |
double y = sqrt((-2 * log(w))/w); |
| 1094 |
m_static_next = m + v2 * y * sqrt(v); |
1096 |
m_static_next = m + v2 * y * sqrt(v); |
| 1095 |
if (fabs(m_static_next) > b) m_static_next = b * (m_static_next)/fabs(m_static_next); |
1097 |
//if next is in bounds, it is valid |
| 1096 |
m_static_nextValid = true; |
1098 |
m_static_nextValid = fabs(m_static_next-m) <= b;; |
| 1097 |
double x1 = m + v1 * y * sqrt(v); |
1099 |
double x1 = m + v1 * y * sqrt(v); |
| 1098 |
if (fabs(x1) > b) x1 = b * (x1)/fabs(x1); |
1100 |
//if x1 is in bounds, return it |
| 1099 |
return x1; |
1101 |
if (fabs(x1-m) <= b) |
|
|
1102 |
{ |
| 1103 |
return x1; |
| 1104 |
} |
| 1105 |
//otherwise try and return m_next if it is valid |
| 1106 |
else if (m_static_nextValid) |
| 1107 |
{ |
| 1108 |
m_static_nextValid = false; |
| 1109 |
return m_static_next; |
| 1110 |
} |
| 1111 |
//otherwise, just run this loop again |
| 1100 |
} |
1112 |
} |
| 1101 |
} |
1113 |
} |
| 1102 |
} |
1114 |
} |