For positive definite $A,B$ why does $AB+BA$ tend to be positive definite?

Your question appears to be based on a false premise. In fact $AB+BA$ does not tend to be positive definite as $n$ increases, even within the particular distribution you happen to be using.

To demonstrate this, here is a simple piece of Mathematica code that implements precisely the numerical experiment you described:

randMat[n_] := With[{
    eigval = RandomReal[{0, 1}, n],
    eigvec = Orthogonalize @ RandomVariate[NormalDistribution[], {n, n}]},
  Transpose[eigvec] . DiagonalMatrix[eigval] . eigvec];

prob[n_] := Table[With[{A = randMat[n], B = randMat[n]},
  PositiveDefiniteMatrixQ[A.B + B.A]], {10000}];

The results obtained (from running 10,000 trials for each value of $n$) are as follows:

+-----+-------------------------------+-----------------------------------+
| n   | probability that AB+BA is pd  | probability that AB+BA is NOT pd  |
+-----+-------------------------------+-----------------------------------+
| 2   | 94.67%                        | 05.33%                            |
+-----+-------------------------------+-----------------------------------+
| 3   | 87.09%                        | 12.91%                            |
+-----+-------------------------------+-----------------------------------+
| 4   | 79.57%                        | 20.43%                            |
+-----+-------------------------------+-----------------------------------+
| 5   | 71.26%                        | 28.74%                            |
+-----+-------------------------------+-----------------------------------+
| 10  | 39.94%                        | 60.06%                            |
+-----+-------------------------------+-----------------------------------+
| 15  | 21.78%                        | 78.22%                            |
+-----+-------------------------------+-----------------------------------+
| 20  | 10.31%                        | 89.69%                            |
+-----+-------------------------------+-----------------------------------+
| 50  | 00.16%                        | 99.84%                            |
+-----+-------------------------------+-----------------------------------+
| 100 | 00.00%                        | 100.00%                           |
+-----+-------------------------------+-----------------------------------+

From this it is clear that you must have accidentally switched the two columns. As $n$ increases, it is becomes rarer for $AB + BA$ to be positive definite.

As an aside, I find it surprising that this question and its currently accepted answer have received a combined 17 upvotes, when seemingly nobody has even tried to replicate the OP's trivial numerical experiment.


$\text{tr}(AB+BA) = 2 \text{tr}(A^{1/2} B A^{1/2}) > 0$, so that may produce some bias toward positive eigenvalues. In particular if you generate your "random" matrices in such a way that the eigenvalues of $AB+BA$ will tend to be concentrated very close together, this may produce the results you observed.

But I tried a different experiment: $A = X^T X$ and $B = Y^T Y$ where $X$ and $Y$ are random $n \times n$ matrices with integer entries in $[-100,100]$.
For the case $n=10$, I found that it was very rare (0 occurrences in 3000 trials) for $AB + BA$ to be positive definite.