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,

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