As you say, four points (without their associated "$t$" values) are not enough to define a planar cubic Bezier curve.
In fact, 6 points are needed. But, on the other hand, the 6 points have to satisfy certain conditions, or else the required curve does not exist. In your case, you know that the curve exists, because it was used to generate the points in the first place.
The full story is given in this paper: J. Kozak, M. Krajnc, Geometric interpolation by planar cubic polynomial curves, Comput. Aided Geom. Des., 24 (2007), pp. 67-78. You can get a copy of the paper at the authors website.
Another approach (instead of the clever techniques suggested in the paper) is to use brute-force numerical methods. Start with 6 points $X_1, \ldots, X_6$ and 6 parameter values $t_1, \ldots, t_6$. We can assume that $t_1 = 0$ and $t_6 = 1$, but $t_2, t_3, t_4, t_5$ are unknown.We want to find control points $P_1, P_2, P_3, P_4$ for our curve. Obviously we can set $P_1 = X_1$ and $P_4 = X_6$, but $P_2$ and $P_3$ are unknown. We can (numerically) find values of $t_2, t_3, t_4, t_5$ and $P_2$ and $P_3$ that minimize the quantity
$$ \sum_{i=2}^{i=5}\Vert C(t_i) - X_i \Vert^2 $$
where $C(t)$ is the Bezier curve with control points $P_1, P_2, P_3, P_4$. This is a fairly nasty non-linear problem, but a good numerical algorithm will usually converge to the desired solution, given a decent starting point.
Another numerical approach: with the same notation we used above, solve the equations:
$$ C(t_i) = X_i \quad (i = 2,3,4,5)$$
Since the $X_i$ are 2D points, there are 8 equations here. We have 8 unknowns -- $t_2, t_3, t_4, t_5$ and the $x$ and $y$ coordinates of $P_2$ and $P_3$. The equations are non-linear, but a good numerical root-finding package should be able to handle them.
A much better approach is to get the Bezier curves from the CAD package. If the design guys won't give you the curves, you should charge them more money for your services, since they are making your life so much more difficult. Their choice -- curves or cash :-)
Let $(x,y)^R=(-y,x)$ represent rotation by $\pi/2$ counterclockwise and
$$
\gamma(t)=(1-t)^3p_0+3t(1-t)^2p_1+3t^2(1-t)p_2+t^3p_3
$$
define a cubic bezier with control points $\{p_0,p_1,p_2,p_3\}$.
Suppose $p_0=(x_0,y_0)$, $p_3=(x_1,y_1)$, and $c=(c_x,c_y)$ are given so that $|p_0-c|=|p_3-c|=r$ ($p_3$ is counterclockwise from $p_0$). Then $p_1=p_0+\alpha(p_0-c)^R$ and $p_2=p_3-\alpha(p_3-c)^R$ where
$$
\alpha=\frac43\tan\left(\frac14\cos^{-1}\left(\frac{(p_0-c)\cdot(p_3-c)}{r^2}\right)\right)
$$
For a quarter of a circle, $\alpha=\frac43(\sqrt2-1)$, and $\gamma$ is no more than $0.00027$ of the radius of the circle off.
Here is a plot of $\gamma$ in red over the quarter circle in black. We really don't see the circle since it is no more than $0.1$ pixels off from $\gamma$ when the radius is $400$ pixels.
$\hspace{3.5cm}$
Computation of $\boldsymbol{\alpha}$
Looking at an arc with an angle of $\theta=\cos^{-1}\left(\frac{(p_0-c)\cdot(p_3-c)}{r^2}\right)$
$\hspace{1.5cm}$
we see that the distance from $c$ to the middle of the arc is
$$
r\cos(\theta/2)+\frac34\alpha r\sin(\theta/2)
$$
we wish to choose $\alpha$ so that this is equal to $r$. Solving for $\alpha$ gives
$$
\begin{align}
\alpha
&=\frac43\frac{1-\cos(\theta/2)}{\sin(\theta/2)}\\
&=\frac43\tan(\theta/4)
\end{align}
$$
A Slight Improvement
Using a circle of radius $1$, the maximum error in radius produced using $\alpha=\frac43\tan(\theta/4)$ is approximately
$$
0.0741\cos^4(\theta/4)\tan^6(\theta/4)
$$
and the error is always positive; that is, the cubic spline never passes inside the circle. Reducing $\alpha$ reduces the midpoint distance by $\frac34\sin(\theta/2)=\frac32\tan(\theta/4)\cos^2(\theta/4)$ times as much, so to distribute the error evenly between the positive and negative, a first guess, assuming that the amplitude of the radius is unchanged, would be to reduce $\alpha$ by $0.0247\cos^2(\theta/4)\tan^5(\theta/4)$.
A bit of investigation shows that, when equalizing the positive and negative swings of the radius, the amplitude increases and that
$$
\alpha=\frac43\tan(\theta/4)-0.03552442\cos^2(\theta/4)\tan^5(\theta/4)
$$
gives pretty even distribution of the error between positive and negative for $\theta\le\pi/2$. The maximum error, both positive and negative, is approximately
$$
0.0533\cos^4(\theta/4)\tan^6(\theta/4)
$$
When $\theta=\pi/2$, this agrees with the article mentioned by bubba in comments.
Note however, that in minimizing the radial error from the circle, the actual variation in radius is increased. Using the simple formula for $\alpha$, which puts the cubic bezier outside the circle, the radius varies by $0.0741\cos^4(\theta/4)\tan^6(\theta/4)$. However, when we minimize the error, the radial variation increases to $0.1066\cos^4(\theta/4)\tan^6(\theta/4)$.
Best Answer
The simplest heuristic approximation method is the one proposed by Tiller and Hanson. They just offset the legs of the control polygon in perpendicular directions:
The blue curve is the original one, and the three blue lines are the legs of its control polygon. We offset these three lines, and intersect/trim them to get the points A and B. The red curve is the offset. It's only an approximation of the true offset, of course, but it's often adequate.
If the approximation is not good enough for your purposes, you split the curve into two, and approximate the two halves individually. Keep splitting until you're happy. You will certainly have to split the original curves at inflexion points, if any, for example.
Here is an example where the approximation is not very good, and splitting would probably be needed:
There is a long discussion of the Tiller-Hanson algorithm plus possible improvements on this web page.
The Tiller-Hanson approach is compared with several others (most of which are more complex) in this paper.
Another good reference, with more up-to-date materials, is this section from the Patrikalakis-Maekawa-Cho book.
For even more references, you can search for "offset" in this bibliography.