The two topics are related. Using gradient descent with DiffSharp worked fine,
but wasn’t ideal. For my purposes, it was too slow, and the gradient approach
was a little overly complex. This led me to investigate if perhaps a
simpler maximization technique like Nelder-Mead would do the job, which in turn
led me to develop Quipu.
Fast forward to today: while Quipu is still in pre-release, its core is fairly
solid now, so I figured I would revisit the problem, and demonstrate how you
could go about using Quipu on a Maximum Likelihood Estimation (or MLE in short)
problem.
In this post, we will begin with a simple problem first, to set the stage. In
the next installment, we will dive into a more complex case, to illustrate why
MLE can be such a powerful technique.
The setup
Imagine that you have a dataset, recording when a piece of equipment
experienced failures. You are interested perhaps in simulating that piece of
equipment, and therefore want to model the time elapsed between failures. As a
starting point, you plot the data as a histogram, and observe something like
this:
It looks like observations fall in between 0 and 8, with a peak around 3.
What we would like to do is estimate a distribution that fits the data. Given
the shape we are observing, a LogNormal distribution is a plausible
candidate. It takes only positive values, which we would expect for durations,
and its density climbs to a peak, and then decreases slowly, which is what we
observe here.
In our last installment, I hit a roadblock. I attempted to
implement Delaunay triangulations using the Bowyer-Watson algorithm,
followed this pseudo-code from Wikipedia,
and ended up with a mostly working F# implementation.
Given a list of points, the code produces a triangulation, but
occasionally the outer boundary of the triangulation is not convex,
displaying bends towards the inside, something that should never
supposed happen for a proper Delaunay triangulation.
While I could not figure out the exact issue, by elimination I
narrowed it down a bit. My guess was that the issue was
probably a missing unstated condition, probably related to the initial
super-triangle. As it turns out, my guess was correct.
The reason I know is, a kind stranger on the internet reached out
with a couple of helpful links (thank you!):
The second link in particular mentions that the Wikipedia page is
indeed missing conditions, and suggests that the initial
super triangle should verify the following property to be valid:
it seems that one should rather demand that the vertices of the
super triangle have to be outside all circumcircles of any three
given points to begin with (which is hard when any three points are almost collinear)
That doesn’t look overly complicated, let’s modify our code
accordingly, and check if this fixes our problem!
During the recent weeks, I have been making slow but steady progress
implementing Delaunay triangulation with the Bowyer-Watson algorithm.
However, as I mentioned in the conclusion of my previous post, I spotted a bug,
which I hoped would be an easy fix, but so far no such luck: it has me stumped.
In this post, I will go over how I approached figuring out the problem, which
is interesting in its own right. My hope is that by talking it through I
might either get an idea, or perhaps someone else will spot what I am missing!
a Delaunay triangulation […] of a set of points in the plane subdivides their
convex hull into triangles whose circumcircles do not contain any of the points;
that is, each circumcircle has its generating points on its circumference, but
all other points in the set are outside of it. This maximizes the size of the
smallest angle in any of the triangles, and tends to avoid sliver triangles.
The Bowyer-Watson algorithm seemed straightforward to implement, so that’s
what I did. For those interested, my current code is here, and it mostly
works. Starting from a list of 20 random points, I can generate a triangulation
like so:
In my last two posts, I did a bit of prep work leading to an implementation of
the Bowyer-Watson algorithm. Now that we have the geometry building blocks
we need, we can attack the core of the algorithm, and perform a Delaunay
triangulation on a list of points.
Rather than attempt to explain what a Delaunay triangulation is, I will
leave that out, and simply illustrate on an example. Starting from a collection
of 20 random points, like so:
… their Delaunay triangulation produces something like this: a mesh of
triangles connecting all the points, without any edges crossing each other, and
where the triangles have “reasonably even” angles:
The Bowyer-Watson algorithm is one way to generate such a
triangulation. In this post, I’ll go over implementing it in F#.
In my last installment, I started revisiting some old code of mine around
Delaunay triangulation, the dual of a Voronoi diagram. My goal is to
implement the Bowyer–Watson algorithm, and perhaps use it for procedural
map generation at some point.
I will follow the pseudo-code outlined on the Wikipedia page, and work my way
through it. Last time I took care of the initialization step, computing an
initial super-triangle large enough to completely contain all the points:
function BowyerWatson (pointList)
triangulation := empty triangle mesh data structure
add super-triangle to triangulation
Today I will tackle the next step:
for each point in pointList do
badTriangles := empty set
for each triangle in triangulation do
if point is inside circumcircle of triangle
add triangle to badTriangles
