Let's put together a simple mathematical model of an accelerometer - from this, we can work out some calibration options.
Ignoring non-linearity and other nasty effects, the output measurement of an accelerometer is given by:
$$\hat{\mathbf{f}} = \mathbf{M} \mathbf{f} + \mathbf{b}_a + \mathbf{n}_a$$
where \$\hat{\mathbf{f}}\$ is the actual measurement, \$\mathbf{b}_a\$ is the accelerometer bias, \$\mathbf{n}_a\$ is a random noise vector, \$\mathbf{f}\$ is the true specific force (i.e. acceleration) and \$\mathbf{M}\$ is the Scale Factor/Misalignment Matrix.
The individual elements of the SFA matrix are:
$$ \mathbf{M} = \begin{bmatrix} S_x && \gamma_{xy} && \gamma_{xz} \\ \gamma_{yx} && S_{yy} && \gamma_{yz} \\ S_x && \gamma_{zy} && S_{zz} \\ \end{bmatrix} $$
So, each scale factor is represented by an \$S\$ and each cross-axis sensitivity is represented by a \$\gamma\$.
Ideally, if the scale factor is 1 and there is no cross-axis sensitivity, then then the resulting matrix is \$\mathbf{M} = \mathbf{I}\$.
Representing it like this allows us to develop a compensation model. If we happen to know \$\mathbf{M}\$ and \$\mathbf{b}_a\$ and assume \$\mathbf{n}_a\$ to be small (i.e. close to zero), we can make a good estimate of the "true" acceleration from the measurements:
$$
\mathbf{f} = \mathbf{M}^{-1}\left(\hat{\mathbf{f}} - \mathbf{b}_a\right)
$$
The trick is, of course, working out \$\mathbf{M}\$ and \$\mathbf{b}_a\$.
I'll describe a procedure called the six position test, which is an easy and cheap way to calibrate an accelerometer. Step 1 is to mount the accelerometer in a rectangular box with perfectly \$90^\circ\$ sides (or as close as you can get). Place this on a perfectly level surface (or, again, as close as you can get) - you'd be surprised how good you can do this.
At this point, we know what the value should be: gravity on the z-accelerometer:
$$
\mathbf{f}_1 = \begin{bmatrix} 0 \\ 0 \\ g\end{bmatrix}
$$
So, this becomes:
$$
\hat{\mathbf{f}}_1 = \mathbf{M} \mathbf{f}_1 + \mathbf{b}_a + \mathbf{n}_a
$$
Noting that \$\hat{\mathbf{f}}_1\$ will close, but not the same as \$\mathbf{f}_1\$
If we put the box on it's head, the force acting is \$-g\$:
$$
\mathbf{f}_2 = \begin{bmatrix} 0 \\ 0 \\ -g\end{bmatrix}
$$
And when placed on one side:
$$
\mathbf{f}_3 = \begin{bmatrix} 0 \\ -g \\ 0\end{bmatrix}
$$
And so on for the remaining three sides.
Now, let's write out one of the equations longhand:
$$
\hat{\mathbf{f}}_1 = \mathbf{M} \mathbf{f}_1 + \mathbf{b}_a + \mathbf{n}_a
= \begin{bmatrix} S_{xx} f_x + \gamma_{xy} f_y + \gamma_{xz} f_z + b_x \\
\gamma_{yx} f_x + S_{yy} f_y + \gamma_{yz} f_z + b_y \\
\gamma_{xz} f_x + \gamma_{yz} f_y + S_{zz} f_z + b_z \end{bmatrix}
$$
And even longer hand (for the first one):
$$
\hat{\mathbf{f}}_1 = \begin{bmatrix}
f_x S_{xx} + f_y \gamma_{xy} + f_z \gamma_{xz} +
0 \gamma_{yx} + 0 S_{yy} + 0 \gamma_{yz} +
0 \gamma_{xz} + 0 \gamma_{yz} + 0 S_{zz} +
1 b_x + 0 b_y + 0 b_z \\
0 S_{xx} + 0 \gamma_{xy} + 0 \gamma_{xz} +
f_x \gamma_{yx} + f_y S_{yy} + f_z \gamma_{yz} +
0 \gamma_{xz} + 0 \gamma_{yz} + 0 S_{zz} +
0 b_x + 1 b_y + 0 b_z \\
0 S_{xx} + 0 \gamma_{xy} + 0 \gamma_{xz} +
0 \gamma_{yx} + 0 S_{yy} + 0 \gamma_{yz} +
f_x \gamma_{xz} + f_y \gamma_{yz} + f_z S_{zz} +
0 b_x + 0 b_y + 1 b_z
\end{bmatrix}
$$
So we can create a stacked vector of the unknowns
$$
\mathbf{z} = \mathbf{A} \mathbf{\beta}
$$
Where
$$
\mathbf{z} =
\begin{bmatrix}
\hat{\mathbf{f}}_1 \\
\hat{\mathbf{f}}_2 \\
\vdots \\
\hat{\mathbf{f}}_6
\end{bmatrix}
$$
And
$$
\mathbf{\beta} = \begin{bmatrix}
S_{xx} \\
\gamma_{xy} \\
\gamma_{xz} \\
\gamma_{yx} \\
S_{yy} \\
\gamma_{yz} \\
\gamma_{xz} \\
\gamma_{yz} \\
S_{zz} \\
b_x \\
b_y \\
b_z \\
\end{bmatrix}
$$
The design matrix is (for one set of measurements):
$$
\hat{A}_1 = \begin{bmatrix}
f_x && f_y && f_z &&
0 && 0 && 0 &&
0 && 0 && 0 &&
1 && 0 && 0 \\
0 && 0 && 0 &&
f_x && f_y && f_z &&
0 && 0 && 0 &&
0 && 1 && 0 \\
0 && 0 && 0 &&
0 && 0 && 0 &&
f_x && f_y && f_z &&
0 && 0 && 1
\end{bmatrix}
$$
Now, once this is setup, one may solve for \$\mathbf{\beta}\$ (and hence sensitivity and bias) via least squares.
A similar procedure may be performed with a robotic arm if you can precisely control the angles - it simply replies knowing the precise gravity at that angle which, if you know the angle, is easy to calculate.
