mathkit  1.5
 All Classes Namespaces Functions Variables Typedefs
mathkit.cpp
1 #include <ctime>
2 #include <cstdlib>
3 #include <cstdarg>
4 #include <climits>
5 #include <sstream>
6 #include <numeric>
7 #include <stdexcept>
8 #include <algorithm>
9 #include <functional>
10 
11 #include "mathkit.hpp"
12 
13 using namespace std;
14 
15 namespace mathkit {
16  bool Epsilon::operator()(double val) const {
17  if (_eps == 0) return val == 0;
18  else return abs(val) < _eps;
19  }
20 
21  double pv(double amount, double rate, double period) {
22  return amount / pow(1 + rate, period);
23  }
24 
25  double fv(double amount, double rate, double period) {
26  return amount * pow(1 + rate, period);
27  }
28 
29  double pv_coef(double rate, double period) {
30  return pow(1 + rate, -period);
31  }
32 
33  double fv_coef(double rate, double period) {
34  return pow(1 + rate, period);
35  }
36 
37  double apv(double annuity, double rate, int period, bool prepaid) {
38  return annuity * apv_coef(rate, period, prepaid);
39  }
40 
41  double afv(double annuity, double rate, int period, bool prepaid) {
42  return annuity * afv_coef(rate, period, prepaid);
43  }
44 
45  double apv_coef(double rate, int period, bool prepaid) {
46  if (rate == 0) return period;
47  double coef = (1 - pow(1 + rate, -period)) / rate;
48  if (prepaid) coef *= (1 + rate);
49  return coef;
50  }
51 
52  double afv_coef(double rate, int period, bool prepaid) {
53  if (rate == 0) return period;
54  double coef = (pow(1 + rate, period) - 1) / rate;
55  if (prepaid) coef *= (1 + rate);
56  return coef;
57  }
58 
59  double spv(const Vector & amount, double rate, bool prepaid) {
60  double pvalue = 0;
61  int period = 1;
62  Vector::const_iterator citer = amount.begin();
63  if (prepaid) {
64  pvalue += *citer;
65  ++citer;
66  }
67  while (citer != amount.end()) {
68  pvalue += pv(*citer, rate, period);
69  ++period;
70  ++citer;
71  }
72  return pvalue;
73  }
74 
75  double sfv(const Vector & amount, double rate, bool prepaid) {
76  double fvalue = 0;
77  int period = 1;
78  Vector::const_reverse_iterator criter = amount.rbegin();
79  if (!prepaid) {
80  fvalue += *criter;
81  ++criter;
82  }
83  while (criter != amount.rend()) {
84  fvalue += fv(*criter, rate, period);
85  ++period;
86  ++criter;
87  }
88  return fvalue;
89  }
90 
91  double comp_rate(double pval, double fval, double period) {
92  return pow(fval / pval, 1 / period) - 1;
93  }
94 
95  double mean(const Vector & data) {
96  double sum = 0;
97  Vector::const_iterator citer;
98  for (citer = data.begin(); citer != data.end(); ++citer) sum += *citer;
99  return sum / data.size();
100  }
101 
102  double median(const Vector & data, bool sorted) {
103  if (!sorted) {
104  Vector temp = data;
105  sort(temp.begin(), temp.end());
106  return median(temp, true);
107  }
108  if (data.size() % 2 == 0) return (data[data.size() / 2] + data[data.size() / 2 - 1]) / 2;
109  else return data[data.size() / 2];
110  }
111 
112  double var(const Vector & data, bool sample) {
113  double avg = mean(data);
114  double dev = 0;
115  Vector::const_iterator citer;
116  for (citer = data.begin(); citer != data.end(); ++citer) dev += (*citer - avg) * (*citer - avg);
117  if (sample) return dev / (data.size() - 1);
118  else return dev / data.size();
119  }
120 
121  double sd(const Vector & data, bool sample) {
122  return sqrt(var(data, sample));
123  }
124 
125  double cov(const Vector & data1, const Vector & data2, bool sample) {
126  double avg1 = mean(data1);
127  double avg2 = mean(data2);
128  double dev = 0;
129  Vector::const_iterator citer1, citer2;
130  citer1 = data1.begin();
131  citer2 = data2.begin();
132  while (citer1 != data1.end() && citer2 != data2.end()) {
133  dev += (*citer1 - avg1) * (*citer2 - avg2);
134  ++citer1;
135  ++citer2;
136  }
137  if (sample) return dev / (data1.size() - 1);
138  else return dev / data1.size();
139  }
140 
141  double cor(const Vector & data1, const Vector & data2) {
142  return cov(data1, data2) / (sd(data1) * sd(data2));
143  }
144 
145  double moment(const Vector & data, int k, bool central) {
146  double result = 0;
147  Vector::const_iterator citer;
148  if (central) {
149  double mu = mean(data);
150  for (citer = data.begin(); citer != data.end(); ++citer)
151  result += pow(*citer - mu, k);
152  }
153  else {
154  for (citer = data.begin(); citer != data.end(); ++citer)
155  result += pow(*citer, k);
156  }
157  return result / data.size();
158  }
159 
160  double skew(const Vector & data) {
161  return moment(data, 3) / pow(moment(data, 2), 1.5);
162  }
163 
164  double kurt(const Vector & data, bool excess) {
165  double result = moment(data, 4) / pow(moment(data, 2), 2);
166  if (excess) result -= 3;
167  return result;
168  }
169 
170  double dnorm(double x, double mu, double sigma) {
171  double z = (x - mu) / sigma;
172  return 1 / (sqrt(2 * pi) * sigma) * exp(-0.5 * z * z);
173  }
174 
175  double pnorm(double x, double mu, double sigma, const Epsilon & eps) {
176  double z = (x - mu) / sigma;
177  double d = dnorm(z);
178  double s = z, sn = z;
179  double p = 0.5 + d * s;
180  int i = 3;
181  while (true) {
182  sn *= z * z / i;
183  s += sn;
184  double pn = 0.5 + d * s;
185  if (eps(pn - p)) {
186  p = pn;
187  break;
188  }
189  p = pn;
190  i += 2;
191  }
192  return p;
193  }
194 
195  double qnorm(double p, double mu, double sigma, const Epsilon & eps) {
196  Epsilon peps(1e-15);
197  if (peps(p - 0.5)) return mu;
198  double lx, rx;
199  if (p > 0.5) {
200  rx = sigma;
201  while (pnorm(rx, mu, sigma, peps) < p) rx += sigma;
202  lx = rx - sigma;
203  }
204  else {
205  lx = -sigma;
206  while (pnorm(lx, mu, sigma, peps) > p) lx -= sigma;
207  rx = lx + sigma;
208  }
209  if (peps(p - pnorm(lx, mu, sigma, peps))) return lx;
210  if (peps(p - pnorm(rx, mu, sigma, peps))) return rx;
211  while (true) {
212  double mid = (lx + rx) / 2;
213  double pmid = pnorm(mid, mu, sigma, peps);
214  if (eps(rx - lx) || peps(p - pmid)) return mid;
215  if (pmid > p) rx = mid;
216  else lx = mid;
217  }
218  }
219 
220  double dt(double x, int n) {
221  return gamma((n + 1) / 2.0) / (sqrt(n * pi) * gamma(n / 2.0)) * pow(1 + x * x / n, (n + 1) / -2.0);
222  }
223 
224  double pt(double x, int n, const Epsilon & eps) {
225  if (x == 0) return 0.5;
226  if (x > 0) {
227  Pair region(-x, x);
228  double p = integrate(bind2nd(ptr_fun<double, int, double>(dt), n), region, eps);
229  return p + (1 - p) / 2;
230  }
231  else {
232  Pair region(x, -x);
233  double p = integrate(bind2nd(ptr_fun<double, int, double>(dt), n), region, eps);
234  return (1 - p) / 2;
235  }
236  }
237 
238  double qt(double p, int n, const Epsilon & eps) {
239  Epsilon peps(1e-10);
240  if (peps(p - 0.5)) return 0;
241  double lx, rx;
242  if (p > 0.5) {
243  rx = 1;
244  while (pt(rx, n, peps) < p) rx += 1;
245  lx = rx - 1;
246  }
247  else {
248  lx = -1;
249  while (pt(lx, n, peps) > p) lx -= 1;
250  rx = lx + 1;
251  }
252  if (peps(p - pt(lx, n, peps))) return lx;
253  if (peps(p - pt(rx, n, peps))) return rx;
254  while (true) {
255  double mid = (lx + rx) / 2;
256  double pmid = pt(mid, n, peps);
257  if (eps(rx - lx) || peps(p - pmid)) return mid;
258  if (pmid > p) rx = mid;
259  else lx = mid;
260  }
261  }
262 
263  double pois(double lmd, int k, bool cum) {
264  if (!cum) return pow(lmd, k) * exp(-lmd) / fac(k);
265  double result = 0;
266  for (int i = 0; i <= k; ++i) result += pois(lmd, i, false);
267  return result;
268  }
269 
270  double binom(int n, double p, int k, bool cum) {
271  if (!cum) return comb(n, k) * pow(p, k) * pow(1 - p, n - k);
272  double result = 0;
273  for (int i = 0; i <= k; ++i) result += binom(n, p, i, false);
274  return result;
275  }
276 
277  double expo(double theta, double x) {
278  if (x > 0) return 1 - exp(-1 / theta * x);
279  else return 0;
280  }
281 
282  double gamma(double x) {
283  if (!(x > 0)) throw invalid_argument("Invalid argument.");
284  double n = x - 1;
285  double ipart;
286  if (modf(n, &ipart) == 0 && n <= UINT_MAX) return fac((unsigned int) n);
287  double p[] = {1.000000000190015, 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 1.208650973866179e-3, -5.395239384953e-6};
288  double s = p[0];
289  for (int i = 1; i <= 6; ++i) s += p[i] / (x + i);
290  return sqrt(2 * pi) / x * s * pow(x + 5.5, x + 0.5) * exp(-x - 5.5);
291  }
292 
293  LineParam linregress(const Vector & xdata, const Vector & ydata) {
294  LineParam regline;
295  regline.slope = cov(xdata, ydata) / var(xdata);
296  regline.intercept = mean(ydata) - regline.slope * mean(xdata);
297  return regline;
298  }
299 
300  Vector seq(double from, double to, double step) {
301  Vector vec;
302  if (step > 0) {
303  if (from > to) return vec;
304  do {
305  vec.push_back(from);
306  from += step;
307  } while (!(from > to));
308  }
309  else if (step < 0) {
310  if (from < to) return vec;
311  do {
312  vec.push_back(from);
313  from += step;
314  } while (!(from < to));
315  }
316  return vec;
317  }
318 
319  Vector linspace(double start, double end, int count) {
320  Vector vec;
321  if (count < 2) return vec;
322  if (count == 2) {
323  vec.push_back(start);
324  vec.push_back(end);
325  }
326  else {
327  double delta = (end - start) / (count - 1);
328  vec.push_back(start);
329  for (int i = 0; i < count - 2; ++i) {
330  start += delta;
331  vec.push_back(start);
332  }
333  vec.push_back(end);
334  }
335  return vec;
336  }
337 
338  double max(const Vector & data) {
339  return *(max_element(data.begin(), data.end()));
340  }
341 
342  double min(const Vector & data) {
343  return *(min_element(data.begin(), data.end()));
344  }
345 
346  double sum(const Vector & data) {
347  return accumulate(data.begin(), data.end(), 0.0);
348  }
349 
350  double prod(const Vector & data) {
351  return accumulate(data.begin(), data.end(), 1.0, multiplies<double>());
352  }
353 
354  Vector add(const Vector & vec1, const Vector & vec2) {
355  Vector result;
356  Vector::const_iterator citer1, citer2;
357  citer1 = vec1.begin();
358  citer2 = vec2.begin();
359  while (citer1 != vec1.end() && citer2 != vec2.end()) {
360  result.push_back(*citer1 + *citer2);
361  ++citer1;
362  ++citer2;
363  }
364  return result;
365  }
366 
367  Vector add(const Vector & vec, double scalar) {
368  Vector result;
369  Vector::const_iterator citer;
370  for (citer = vec.begin(); citer != vec.end(); ++citer) result.push_back(*citer + scalar);
371  return result;
372  }
373 
374  Vector operator+(const Vector & vec1, const Vector & vec2) {
375  return add(vec1, vec2);
376  }
377 
378  Vector operator+(const Vector & vec, double scalar) {
379  return add(vec, scalar);
380  }
381 
382  Vector operator+(double scalar, const Vector & vec) {
383  return add(vec, scalar);
384  }
385 
386  Vector sub(const Vector & vec1, const Vector & vec2) {
387  Vector result;
388  Vector::const_iterator citer1, citer2;
389  citer1 = vec1.begin();
390  citer2 = vec2.begin();
391  while (citer1 != vec1.end() && citer2 != vec2.end()) {
392  result.push_back(*citer1 - *citer2);
393  ++citer1;
394  ++citer2;
395  }
396  return result;
397  }
398 
399  Vector sub(const Vector & vec, double scalar, bool dir) {
400  Vector result;
401  Vector::const_iterator citer;
402  for (citer = vec.begin(); citer != vec.end(); ++citer) {
403  double element = *citer - scalar;
404  if (!dir) element = -element;
405  result.push_back(element);
406  }
407  return result;
408  }
409 
410  Vector operator-(const Vector & vec1, const Vector & vec2) {
411  return sub(vec1, vec2);
412  }
413 
414  Vector operator-(const Vector & vec, double scalar) {
415  return sub(vec, scalar);
416  }
417 
418  Vector operator-(double scalar, const Vector & vec) {
419  return sub(vec, scalar, false);
420  }
421 
422  Vector mul(const Vector & vec1, const Vector & vec2) {
423  Vector result;
424  Vector::const_iterator citer1, citer2;
425  citer1 = vec1.begin();
426  citer2 = vec2.begin();
427  while (citer1 != vec1.end() && citer2 != vec2.end()) {
428  result.push_back(*citer1 * *citer2);
429  ++citer1;
430  ++citer2;
431  }
432  return result;
433  }
434 
435  Vector mul(const Vector & vec, double scalar) {
436  Vector result;
437  Vector::const_iterator citer;
438  for (citer = vec.begin(); citer != vec.end(); ++citer) result.push_back(*citer * scalar);
439  return result;
440  }
441 
442  Vector operator*(const Vector & vec1, const Vector & vec2) {
443  return mul(vec1, vec2);
444  }
445 
446  Vector operator*(const Vector & vec, double scalar) {
447  return mul(vec, scalar);
448  }
449 
450  Vector operator*(double scalar, const Vector & vec) {
451  return mul(vec, scalar);
452  }
453 
454  Vector div(const Vector & vec1, const Vector & vec2) {
455  Vector result;
456  Vector::const_iterator citer1, citer2;
457  citer1 = vec1.begin();
458  citer2 = vec2.begin();
459  while (citer1 != vec1.end() && citer2 != vec2.end()) {
460  result.push_back(*citer1 / *citer2);
461  ++citer1;
462  ++citer2;
463  }
464  return result;
465  }
466 
467  Vector div(const Vector & vec, double scalar, bool dir) {
468  Vector result;
469  Vector::const_iterator citer;
470  for (citer = vec.begin(); citer != vec.end(); ++citer) {
471  if (dir) result.push_back(*citer / scalar);
472  else result.push_back(scalar / *citer);
473  }
474  return result;
475  }
476 
477  Vector operator/(const Vector & vec1, const Vector & vec2) {
478  return div(vec1, vec2);
479  }
480 
481  Vector operator/(const Vector & vec, double scalar) {
482  return div(vec, scalar);
483  }
484 
485  Vector operator/(double scalar, const Vector & vec) {
486  return div(vec, scalar, false);
487  }
488 
489  double dot_prod(const Vector & vec1, const Vector & vec2) {
490  double result = 0;
491  Vector::const_iterator citer1, citer2;
492  citer1 = vec1.begin();
493  citer2 = vec2.begin();
494  while (citer1 != vec1.end() && citer2 != vec2.end()) {
495  result += *citer1 * *citer2;
496  ++citer1;
497  ++citer2;
498  }
499  return result;
500  }
501 
502  Vector cross_prod(const Vector & vec1, const Vector & vec2) {
503  Vector vec;
504  vec.push_back(vec1[1] * vec2[2] - vec1[2] * vec2[1]);
505  vec.push_back(vec1[2] * vec2[0] - vec1[0] * vec2[2]);
506  vec.push_back(vec1[0] * vec2[1] - vec1[1] * vec2[0]);
507  return vec;
508  }
509 
510  double norm(const Vector & vec) {
511  return sqrt(dot_prod(vec, vec));
512  }
513 
514  string to_str(const Vector & vec, string sep, string delim) {
515  ostringstream oss;
516  Vector::size_type n = vec.size();
517  Vector::size_type last = n - 1;
518  for (Vector::size_type i = 0; i < n; ++i) {
519  oss << vec[i];
520  if (i != last) oss << sep;
521  }
522  string result = oss.str();
523  if (delim.size() == 2) result = delim[0] + result + delim[1];
524  return result;
525  }
526 
527  Vector make_vec(int n, ...) {
528  Vector vec;
529  va_list argptr;
530  va_start(argptr, n);
531  for (int i = 0; i < n; ++i) vec.push_back(va_arg(argptr, double));
532  va_end(argptr);
533  return vec;
534  }
535 
536  Vector load(istream & ins, string delim) {
537  Vector data;
538  string line;
539  if (getline(ins, line)) {
540  if (delim != " ") {
541  string::size_type i;
542  while ((i = line.find(delim)) != string::npos) line.replace(i, delim.size(), " ");
543  }
544  istringstream iss(line);
545  while (!iss.eof()) {
546  double item;
547  iss >> item;
548  data.push_back(item);
549  }
550  }
551  return data;
552  }
553 
554  Table load(istream & ins, int nrow, int ncol, string delim) {
555  Table data(ncol, Vector());
556  for (int i = 0; i < nrow; ++i) {
557  Vector row = load(ins, delim);
558  for (int j = 0; j < ncol; ++j) data[j].push_back(row[j]);
559  }
560  return data;
561  }
562 
563  void save(ostream & outs, const Vector & data, string delim) {
564  Vector::const_iterator citer, last;
565  last = data.end() - 1;
566  int prec = outs.precision();
567  outs.precision(15);
568  for (citer = data.begin(); citer != data.end(); ++citer) {
569  outs << *citer;
570  if (citer != last) outs << delim;
571  }
572  outs << endl;
573  outs.precision(prec);
574  }
575 
576  void save(ostream & outs, const Table & data, string delim) {
577  int nrow = data[0].size();
578  int ncol = data.size();
579  for (int i = 0; i < nrow; ++i) {
580  Vector row;
581  for (int j = 0; j < ncol; ++j) row.push_back(data[j][i]);
582  save(outs, row, delim);
583  }
584  }
585 
586  ostream & operator<<(ostream & outs, const Vector & vec) {
587  save(outs, vec);
588  return outs;
589  }
590 
591  istream & operator>>(istream & ins, Vector & vec) {
592  vec = load(ins);
593  return ins;
594  }
595 
596  double randf(double low, double high) {
597  return rand() / (double) RAND_MAX * (high - low) + low;
598  }
599 
600  Vector randf(double low, double high, int n) {
601  Vector result;
602  for (int i = 0; i < n; ++i) result.push_back(randf(low, high));
603  return result;
604  }
605 
606  int randi(int low, int high) {
607  return (int) floor(rand() / ((double) RAND_MAX + 1) * ((double) high - low + 1) + low);
608  }
609 
610  Vector randi(int low, int high, int n) {
611  Vector result;
612  for (int i = 0; i < n; ++i) result.push_back(randi(low, high));
613  return result;
614  }
615 
616  bool randp(double p) {
617  return !(randf(0, 1) > p);
618  }
619 
620  void randseed() {
621  srand((unsigned int) time(NULL));
622  }
623 
624  void randseed(unsigned int s) {
625  srand(s);
626  }
627 
628  double rnorm(double mu, double sigma) {
629  double u = randf(0, 1);
630  double v = randf(0, 1);
631  double z = sqrt(-2 * log(u)) * cos(2 * pi * v);
632  return mu + sigma * z;
633  }
634 
635  Vector rnorm(int n, double mu, double sigma) {
636  Vector result;
637  for (int i = 0; i < n; ++i) result.push_back(rnorm(mu, sigma));
638  return result;
639  }
640 
641  int rpois(double lmd) {
642  double L = exp(-lmd);
643  double p = 1;
644  int k = 0;
645  do {
646  k += 1;
647  p *= randf(0, 1);
648  } while (p > L);
649  return k - 1;
650  }
651 
652  Vector rpois(double lmd, int n) {
653  Vector result;
654  for (int i = 0; i < n; ++i) result.push_back(rpois(lmd));
655  return result;
656  }
657 
658  int rbinom(int n, double p) {
659  int k = 0;
660  for (int i = 0; i < n; ++i)
661  if (randp(p)) k += 1;
662  return k;
663  }
664 
665  Vector rbinom(int n, double p, int count) {
666  Vector result;
667  for (int i = 0; i < count; ++i) result.push_back(rbinom(n, p));
668  return result;
669  }
670 
671  double rexp(double theta) {
672  return -theta * log(randf(0, 1));
673  }
674 
675  Vector rexp(double theta, int n) {
676  Vector result;
677  for (int i = 0; i < n; ++i) result.push_back(rexp(theta));
678  return result;
679  }
680 
681  double prec(double num, int ndec) {
682  double intpart;
683  double decpart;
684  decpart = modf(num, &intpart);
685  stringstream ss;
686  ss.precision(ndec + 1);
687  ss << decpart + 1;
688  ss >> decpart;
689  return intpart + decpart - 1;
690  }
691 
692  Vector prec(const Vector & vec, int ndec) {
693  Vector result;
694  Vector::const_iterator citer;
695  for (citer = vec.begin(); citer != vec.end(); ++citer)
696  result.push_back(prec(*citer, ndec));
697  return result;
698  }
699 
700  double fac(unsigned int n) {
701  if (n == 0) return 1;
702  else {
703  double result = 1;
704  while (n > 0) {
705  result *= n;
706  --n;
707  }
708  return prec(result, 0);
709  }
710  }
711 
712  double perm(unsigned int m, unsigned int n) {
713  double result = 1;
714  while (n > 0) {
715  result *= m;
716  --m;
717  --n;
718  }
719  return prec(result, 0);
720  }
721 
722  double comb(unsigned int m, unsigned int n) {
723  if (m - n < n) return comb(m, m - n);
724  return prec(perm(m, n) / fac(n), 0);
725  }
726 
727  unsigned int gcd(unsigned int a, unsigned int b) {
728  while (b != 0) {
729  int temp = b;
730  b = a % b;
731  a = temp;
732  }
733  return a;
734  }
735 
736  unsigned int lcm(unsigned int a, unsigned int b) {
737  return a * b / gcd(a, b);
738  }
739 }