What's missing in a lot of these answers is that you may not want to fit a single cubic Bézier curve to your data. More generally, you would like to fit a sequence of cubic Bézier curves, i.e., a piecewise cubic Bézier fit, to an arbitrary set of data.
There's a nice thesis dating from 1995, complete with MATLAB code, that does this:
% Lane, Edward J. Fitting Data Using Piecewise G1 Cubic Bezier Curves.
% Thesis, NAVAL POSTGRADUATE SCHOOL MONTEREY CA, 1995
http://www.dtic.mil/dtic/tr/fulltext/u2/a298091.pdf
To use this, you must, at minimum, specify the number of knot points, i.e., the number data points that will be used by the optimization routines to make this fit. Optionally, you can specify the knot points themselves, which increases the reliability of a fit. The thesis shows some pretty tough examples. Note that Lane's approach guarantees G1 continuity (directions of adjacent tangent vectors are identical) between the cubic Bézier segments, i.e., smooth joints. However, there can be discontinuities in curvature (changes in direction of second derivative).
I have reimplemented the code, updating it to modern MATLAB (R2015b). Contact me if you would like it.
Here's an example of using just three knot points (chosen automatically by the code) the fits two cubic Bézier segments to a Lissajous figure.