Anonymous user
Talk:Total circles area: Difference between revisions
→Monte Carlo
Line 3:
A semi-important point and semi-nitpicking for doing Monte Carlo: you need to estimate the error in a credible way given the number of samples, and should ''not'' output more significant digits than the confidence allows. The current C code sample size of 100M can't possibly allow 8 decimal point precision. --[[User:Ledrug|Ledrug]] 04:29, 16 September 2012 (UTC)
: Right, in the C/Python outputs only the first few digits are correct. Feel free to add a bit of error estimation code.
:: This is basically what I used on reddit, but I don't feel like changing the example on the task page because MC really isn't the right tool for this task.
:: <lang c>#include <stdio.h>
#include <stdlib.h>
#include <tgmath.h>
#include <time.h>
typedef float flt;
flt data[][3] = {
{ 1.6417233788, 1.6121789534, 0.0848270516},
...
};
#define N sizeof(data)/sizeof(data[0])
int main(void)
{
flt x0 = 1/0., x1 = -1/0., y0 = 1/0., y1 = -1/0.;
inline flt min(flt x, flt y) { return x < y ? x : y; }
inline flt max(flt x, flt y) { return x > y ? x : y; }
inline flt sq(flt x) { return x * x; }
inline flt rndf(flt x, flt d) { return x + d * rand() / RAND_MAX; }
inline int covered(flt x, flt y) {
for (int i = 0; i < N; i++)
if (sq(x - data[i][0]) + sq(y - data[i][1]) < data[i][2])
return 1;
return 0;
}
for (int i = 0; i < N; i++) {
flt *p = data[i];
x0 = min(x0, p[0] - p[2]);
x1 = max(x1, p[0] + p[2]);
y0 = min(y0, p[1] - p[2]);
y1 = max(y1, p[1] + p[2]);
p[2] *= p[2];
}
x1 -= x0, y1 -= y0;
flt total = x1 * y1;
unsigned long to_try = 0x10000, tries = 0, hit = 0;
srand(time(0));
while (1) {
hit += covered(rndf(x0, x1), rndf(y0, y1));
if (++tries == to_try) {
flt area = total * hit / tries;
flt r = (flt) hit / tries;
flt s = area * sqrt(r * (1 - r) / tries);
printf("%.4f +/- %.4f (%lu samples)\n", area, s, tries);
// stop at 3 sigmas
if (s * 3 <= 1e-3) break;
to_try *= 2;
}
}
return 0;
}</lang> --[[User:Ledrug|Ledrug]] 22:09, 16 September 2012 (UTC)
| |||