25 #include <MathAlgo/MonteCarloTransform.hh>
26 #include <MathAlgo/SVD.hh>
27 #include <MathAlgo/Lapack.hh>
28 #include <Base/Common/CompareFloatingPoint.hh>
29 #include <Base/Debug/Exception.hh>
30 #include <Base/Math/Matrix.hh>
40 UseQuasiRandomNumbers_(true), uniformDistribution_(false)
62 vector<Vector<double> > src_samples;
63 const int num =
GetSamples_(src_mean, src_cov, src_samples);
64 dst_samples.resize(num);
68 of.open(
"dst_samples.txt");
71 #ifdef WRITE_DST_SAMPLES
72 ofstream off(
"dst_samples_p4.txt");
75 for (
int i=0; i<num; i++){
76 if ((res =
Transform_(src_samples[i], dst_samples[i-num_errs]))!=0){
83 of << setw(5) << i <<
"\t";
84 for (
unsigned d=0; d<dst_samples[i-num_errs].Size(); d++){
85 of << dst_samples[i-num_errs][d] <<
"\t";
90 #ifdef WRITE_DST_SAMPLES
92 off << setw(5) << i <<
"\t";
93 for (
unsigned d=0; d<dst_samples[i-num_errs].Size(); d++){
94 off << dst_samples[i-num_errs][d] <<
"\t";
102 BCDOUT(D_MCT_NUM_SAMPLES, num_errs<<
" errors in transform functions, "
103 "using less samples ("<<num-num_errs<<
")\n");
104 dst_samples.resize(num-num_errs);
111 #ifdef WRITE_DST_SAMPLES
114 return (num_errs==0)?(0):(-num_errs);
122 if (src_cov.
IsZero(numeric_limits<double>::epsilon())){
124 const int dim = dst_mean.
Size();
128 vector<Vector<double> > dst_samples;
129 int res =
Transform(src_mean, src_cov, dst_samples);
131 BIASERR(
"error transforming samples. error code: "<<res);
134 const int num = (int)dst_samples.size();
136 const int dim = (int)dst_samples[0].Size();
141 for (
int i=0; i<num; i++){
142 dst_mean += dst_samples[i];
145 dst_mean /= (double)num;
148 for (
int d=0; d<dim; d++) {
149 INFNANCHECK(dst_mean[d]);
155 for (
int i=0; i<num; i++){
156 diff = dst_samples[i] - dst_mean;
159 dst_cov /= (double)(num-1);
162 for (
int row=0; row<dst_cov.
num_rows(); row++) {
163 for (
int col=0; col<dst_cov.
num_cols(); col++) {
164 INFNANCHECK(dst_cov[row][col]);
180 vector<Vector<double> > eigenvec;
183 double factor = src_cov.
Trace();
190 BIASERR(
"error computing eigenvalues: "<<res);
193 BCDOUT(D_MCT_EIGENVALVEC,
"eignevalues: "<<eigenval<<endl);
195 const int dim = (int)src_mean.
Size();
199 for (
int d=0; d<dim; d++){
200 BCDOUT(D_MCT_EIGENVALVEC, d<<
": "<<eigenvec[d]<<endl);
203 if (eigenval[d]<0.0 && eigenval[d]>-1e-22){
204 BCDOUT(D_MCT_EIGENVALVEC,
"["<<d<<
"] numerical inaccurate eigenvalue: "
205 <<eigenval[d]<<endl<<
"corresp eigenvec.NormL2(): "
206 <<eigenvec[d].NormL2()<<endl);
210 if (eigenval[d]<0.0) {
211 cout <<
"["<<d<<
"] illegal eigenvalue: "<<eigenval[d]<<endl;
212 cout <<
"corresp eigenvec: "<<eigenvec[d]<<endl;
215 BIASASSERT(eigenval[d]>=0.0);
217 eigenval[d] *= factor;
219 eigenval[d] = sqrt(eigenval[d]);
221 if (
Equal(eigenvec[d].NormL2(), 0.0)) {
222 cout <<
"illegal eigenvector: "<<eigenvec[d]<<endl;}
224 BIASASSERT(!
Equal(eigenvec[d].NormL2(), 0.0));
225 eigenvec[d] /= eigenvec[d].NormL2();
232 of.open(
"src_samples.txt");
237 cout <<
"using: "<<num<<
" samples\n";
238 vector<Vector<double> > quasi_std_normal_distr;
239 vector<double> normal_distr;
244 BIASWARNONCE(
"uniform-distribution sampling is experimental"
245 ", multi-dim. case needs verification");
248 themax(dim, sqrt(3.0));
250 quasi_std_normal_distr);
262 BCDOUT(D_MCT_NUM_SAMPLES,
"QuasiRandom with mean zero and cov=id computed "
263 <<testmean<<
" and "<<testcov);
268 num*dim, normal_distr);
277 BCDOUT(D_MCT_NUM_SAMPLES,
"Random with mean zero and cov=id computed "
278 <<testmean<<
" and "<<testcov);
284 src_samples.resize(num, zero);
285 for (
int i=0; i<num; i++){
286 for (
int d=0; d<dim; d++){
289 r = quasi_std_normal_distr[i][d] * eigenval[d];
292 r = normal_distr[jj++]*eigenval[d];
303 src_samples[i] += r * eigenvec[d];
315 src_samples[i] += src_mean;
327 for (
int d=0; d<dim; d++){
328 of << src_samples[i][d] <<
"\t";
333 BCDOUT(D_MCT_NUM_SAMPLES,
"x"<<flush);
340 BCDOUT(D_MCT_NUM_SAMPLES,
"\nusing "<<num<<
" src_samples\n");
373 double skewnness)
const
375 BIASERR(
"unfinished");
bool IsZero(double eps=0.0) const
Checks if the matrix is a null matrix.
Subscript num_cols() const
int GetQuasiUniformDistributed(const int dim, const Vector< double > &min, const Vector< double > &max, Vector< double > &res)
On subsequent calls returns a Sobol' sequence of quasi random numbers 'uniformly distributed' in a N-...
Matrix< T > & newsize(Subscript M, Subscript N)
unsigned int Size() const
length of the vector
double GetUniformDistributed(const double min, const double max)
on succesive calls return uniform distributed random variable between min and max ...
void SetZero()
equivalent to matrix call
bool DebugLevelIsSet(const long int lv) const
Vector< T > & newsize(Subscript N)
Matrix< T > OuterProduct(const Vector< T > &v) const
outer product, constructs a matrix.
void SetZero()
Sets all values to zero.
long int NewDebugLevel(const std::string &name)
creates a new debuglevel
bool Equal(const T left, const T right, const T eps)
comparison function for floating point values See http://www.boost.org/libs/test/doc/components/test_...
Subscript num_rows() const
double GetNormalDistributed(const double mean, const double sigma)
on succesive calls return normal distributed random variable with mean and standard deviation sigma ...
BIAS::Vector< double > Upper_symmetric_eigenvalue_solve(const BIAS::Matrix< double > &A)
Solve symmetric eigenvalue problem (eigenvalues only)
class for producing random numbers from different distributions
void GetMeanAndCovariance(const std::vector< double > &samples, double &mean, double &covariance) const
check-function to compute mean and cov from samples
int GetQuasiNormalDistributed(const int dim, Vector< double > &res)
On subsequent calls returns a Sobol' sequence of quasi random numbers 'normal distributed' with mean ...