Compensating for Accelerometer Cross-Axis Sensitivity
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 &amp;&amp; \gamma_{xy} &amp;&amp; \gamma_{xz} \\ \gamma_{yx} &amp;&amp; S_{yy} &amp;&amp; \gamma_{yz} \\ S_x &amp;&amp; \gamma_{zy} &amp;&amp; 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 &amp;&amp; f_y &amp;&amp;  f_z  &amp;&amp;
0   &amp;&amp; 0   &amp;&amp;  0    &amp;&amp; 
0   &amp;&amp; 0   &amp;&amp;  0    &amp;&amp;
1   &amp;&amp; 0   &amp;&amp;  0    \\
0   &amp;&amp; 0   &amp;&amp;  0    &amp;&amp; 
f_x &amp;&amp; f_y &amp;&amp;  f_z  &amp;&amp;
0   &amp;&amp; 0   &amp;&amp;  0    &amp;&amp;
0   &amp;&amp; 1   &amp;&amp;  0    \\
0   &amp;&amp; 0   &amp;&amp;  0    &amp;&amp; 
0   &amp;&amp; 0   &amp;&amp;  0    &amp;&amp;
f_x &amp;&amp; f_y &amp;&amp;  f_z  &amp;&amp;
0   &amp;&amp; 0   &amp;&amp;  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.