środa, 27 marca 2013

Rzutki, Pi i PostGIS.

Przybliżenie liczby Pi metodą Monte Carlo. Czyli gra w rzutki. Wieszamy tarczę do gry w rzutki na kwadratowej desce o boku równym średnicy tarczy. Po kilku kolejkach gry liczymy stosunek rzutek które trafiły w tarczę do liczby rzutek które trafiły w tarczę lub deskę i mnożymy razy cztery by uzyskać liczbę Pi. Im więcej rzutek tym dokładniejsze przybliżenie.

Na początek musimy się wygenerować losowe punkty w PostGISie:

SELECT st_makepoint((2.0*random()-1),(2.0*random()-1))
FROM generate_series(1,100);
To zapytanie wygeneruje nam sto punktów których x,y leży w przedziale <-1.0,1.0), następnie obliczamy stosunek punktów leżących w kole do ogółu:
SELECT sum(st_dwithin(point,st_makepoint(0,0), 1.0)::int)::float/count(*)*4
FROM
(SELECT st_makepoint((2.0*random()-1),(2.0*random()-1)) AS point
FROM generate_series(1,100)) AS q;
Zmierzymy przybliżenie dla stu punktów a następnie dla tysiąca, dziesięciu tysięcy i stu tysięcy (wizualizacja w OpenJUMP).

100 punktów Pi=3,36
 1000 punktów Pi=3,156
 10 tysięcy punktów Pi=3.1336
100 tysięcy punktów Pi=3.141



I na sam koniec, zapytanie do generowania danych do wykresu:
WITH RECURSIVE t(n,point, circ,pi) AS (
    SELECT 1 as n,point, ST_DWithin(point,ST_MAKEPOINT(0,0), 1.0)::int AS circ, ST_DWithin(point,ST_MAKEPOINT(0,0), 1.0)::int::float/1*4 as pi
    FROM (SELECT ST_MAKEPOINT((2.0*random()-1),(2.0*random()-1)) as point) as q
  UNION ALL
    SELECT t.n+1 as n,q.point, circ + ST_DWithin(q.point,ST_MAKEPOINT(0,0), 1.0)::int,
    ((circ + ST_DWithin(q.point,ST_MAKEPOINT(0,0), 1.0)::int::float)/(n+1))*4 as pi
    FROM (SELECT ST_MAKEPOINT((2.0*random()-1),(2.0*random()-1)) as point) as q,t
    WHERE n<10000  
)
SELECT * FROM t;