| | 80 | double normcoef(const epdf *ep, const vec &xb, const vec &yb, |
| | 81 | int xn, int yn) { |
| | 82 | mat Pdf(xn + 1, yn + 1); |
| | 83 | vec rgr(2); |
| | 84 | |
| | 85 | double xstep = (xb(1) - xb(0)) / xn; |
| | 86 | double ystep = (yb(1) - yb(0)) / yn; |
| | 87 | |
| | 88 | int i = 0; |
| | 89 | for (double x = xb(0); x <= xb(1); x += xstep, i++) { |
| | 90 | rgr(0) = x; |
| | 91 | int j = 0; |
| | 92 | for (double y = yb(0); y <= yb(1); y += ystep, j++) { |
| | 93 | rgr(1) = y; |
| | 94 | Pdf(i, j) = exp(ep->evallog(rgr)); |
| | 95 | } |
| | 96 | } |
| | 97 | |
| | 98 | return sumsum(Pdf) * xstep * ystep; |
| | 100 | |
| | 101 | vec num_mean2(const epdf *ep, const vec &xb, const vec &yb, |
| | 102 | int xn, int yn) { |
| | 103 | mat Pdf(xn + 1, yn + 1); |
| | 104 | vec rgr(2); |
| | 105 | |
| | 106 | double xstep = (xb(1) - xb(0)) / xn; |
| | 107 | double ystep = (yb(1) - yb(0)) / yn; |
| | 108 | |
| | 109 | vec Mu(xn + 1); |
| | 110 | vec Si(yn + 1); |
| | 111 | |
| | 112 | int i = 0; |
| | 113 | for (double x = xb(0); x <= xb(1); x += xstep, i++) { |
| | 114 | Mu(i) = x; |
| | 115 | rgr(0) = x; |
| | 116 | int j = 0; |
| | 117 | for (double y = yb(0); y <= yb(1); y += ystep, j++) { |
| | 118 | Si(j) = y; |
| | 119 | rgr(1) = y; |
| | 120 | Pdf(i, j) = exp(ep->evallog(rgr)); |
| | 121 | } |
| | 122 | } |
| | 123 | |
| | 124 | vec fm = sum(Pdf, 2); |
| | 125 | vec fs = sum(Pdf, 1); |
| | 126 | return vec_2(Mu * fm / sum(fm), Si * fs / sum(fs)); |
| | 127 | } |
| | 128 | |
| | 129 | } |