Monte Carlo approximation of the area of heart

Problem

A popular mathematical equation for 2D heart is the following,

\[f(x,y) = (x^2 + y^2 - 1)^3 - x^2 y^3 = 0\]

Any $(x,y)$ values that result in $f(x,y) < 0$ represent the coordinates of a point that falls inside the heart shape. Write a script that generates a heart shape via a set of random points uniformly distributed inside the heart shape. Estimate the area of the heart by drawing an appropriate box around the heart and sampling points randomly uniformly from this box and counting those that fall inside the heart. Here is a plot of the heart,

Solution

Here is an example Python code that estimates $\pi$,

import matplotlib.pyplot as plt 
import seaborn as sns
import numpy as np
sns.set()


# set the number of simulated points

nsim = 20000

# set the sampling box range and limits

xlim = [-1.5,1.5] # sampling box x limits
ylim = [-1.5,1.5] # sampling box y limits

xrange = xlim[1] - xlim[0]
yrange = ylim[1] - ylim[0]

x = np.random.random(nsim) * xrange + xlim[0]
y = np.random.random(nsim) * yrange + ylim[0]

# define heart

def getHeart(x,y): return (x**2 + y**2 - 1)**3 - x**2 * y**3

# select points that are inside 3D heart

hdistSq = np.zeros(nsim)
hpoints = np.zeros((nsim,2))
npoints = 0
for i in range(nsim):
    distSq = getHeart(x[i],y[i])
    if distSq <= 0:
        npoints += 1
        hpoints[npoints-1,:] = np.array([x[i],y[i]])
        hdistSq[npoints-1] = distSq

# plot heart

import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

fig = plt.figure(figsize = (4.5,4), dpi = 200)
ax = plt.axes()

ax.scatter( hpoints[0:npoints,0]
          , hpoints[0:npoints,1]
          , color = "red"
          , s = 0.5
          )
ax.set_xlim(xlim)
ax.set_ylim(ylim)
plt.tight_layout()
plt.savefig("heart.png")
plt.show()

print("area = {}".format( xrange * yrange * npoints / nsim ))

area = 3.661299