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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
|
#include <cinttypes>
#include <cmath>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <iostream>
#include <iomanip>
long double lfactorial(unsigned i) { return lgamma((long double)(i + 1)); }
long double lbinom(unsigned N, unsigned n) {
return lfactorial(N) - lfactorial(n) - lfactorial(N-n);
}
long double B(unsigned N, unsigned n, long double p) {
return expl(lbinom(N, n) + log(p) * n + log(1-p) * (N - n));
}
int main(int argc, char *argv[]) {
std::vector<std::string> all_args;
if (argc > 1) {
all_args.assign(argv + 1, argv + argc);
}
for (std::string filename : all_args) {
std::ifstream in_file;
in_file.open(filename);
std::vector<std::vector<double>> dat;
unsigned line_num = 0;
std::string line;
std::getline(in_file, line);
while (std::getline(in_file, line)) {
dat.resize(line_num + 1);
std::istringstream iss(line);
double a;
while (iss >> a) {
dat[line_num].push_back(a);
}
line_num++;
}
in_file.close();
unsigned n_moms = dat.size();
unsigned n_ps = 100;
double ε = 1e-5;
std::vector<std::vector<double>> output(n_moms);
for (unsigned i = 0; i < n_moms; i++) {
output[i].resize(n_ps);
}
for (unsigned i = 0; i < n_moms; i++) {
unsigned N = dat[i].size();
for (unsigned j = 0; j < n_ps; j++) {
long double p = 0.5 - 0.5 * expl((j + 1) * logl(ε) / n_ps);
long double mom_p = pow(1 - p, N);
long double tot = pow(1 - p, N);
for (unsigned n = 1; n <= N; n++) {
mom_p += B(N, n, p) * dat[i][n - 1];
tot += B(N, n, p);
}
output[i][j] = mom_p;
}
}
std::ofstream out_file;
out_file.open(filename + ".ps");
for (unsigned i = 0; i < n_moms; i++) {
for (unsigned j = 0; j < output[i].size(); j++) {
long double p = 0.5 - 0.5 * expl((j + 1) * logl(ε) / n_ps);
out_file << std::fixed << std::setprecision(8) << p << " " << output[i][j] << " ";
}
out_file << "\n";
}
out_file.close();
}
return 0;
}
|
