I am trying to clarify the following:

Suppose P is the $N\times N$ transition matrix of a finite-dimensional Markov Chain, with invariant distribution given by N-vector $\mu$ (i.e. $\mu^T=\mu^T P$). Define the time-reverse chain with transition matrix $\hat P$ as:
$\hat p\_{ij}=\frac{\mu\_j}{\mu\_i} p\_{ji}$

Now define the reversible transition matrix $R=(P+\hat P)/2$. It can be shown that this reversible chain has the same invariant measure $\mu$.

Now if we define an inner product space with $<x,y>\_\mu=\Sigma x\_i y\_i \mu\_i$, then $R$ is self-adjoint in this space, with real eigenvalues and orthogonal eigenvectors (w.r.t the inner product defined above).

Somehow, numerically, when I try the above computation on a randomly generated Transition matrix, I do not get a self-adjoint R, and the eigenvectors are not orthogonal in the sense described above. Can anyone tell me if there is an error in above statements ? 

Edit (in response to Bryon): Let me elaborate on my statements to see if there is any error

Let us define $D\_\mu=Diag(\mu)$, a diagonal matrix with $\mu$ on the diagonal. Now define $A=D\_\mu^{1/2}RD\_\mu^{-1/2}$. It can be checked that A is symmetric matrix.

Then let eigenvectors of A be $\phi\_j$, and eigenvectors of R be $f\_j$. 
The relation between the two is:
$f\_j=D\_\mu^{-1/2}\phi\_j$

Now consider $\delta\_{ij}=\langle \phi\_i,\phi\_j\rangle=\langle D\_\mu^{1/2} f\_i,D\_\mu^{1/2} f\_j\rangle=\langle f\_i,f\_j\rangle\_\mu$

Edit #2 (Response to Did):

The reversibility of R w.r.t $\mu$ implies $R\_{ij}\mu\_i=R\_{ji}\mu\_j \:\: \forall i,j$

Now, the self-adjointness condition is $\langle Rx,y\rangle\_\mu=\langle x,Ry\rangle\_\mu$

L.H.S: $\Sigma\_i \Sigma\_j R\_{ij}x\_j y\_i \mu\_i=\Sigma\_i \Sigma\_j R\_{ji}x\_j y\_i \mu\_j$ (by the reversibility property above).

Then interchanging $i$ and $j$, we get , L.H.S=$\Sigma\_i \Sigma\_j R\_{ij}x\_i y\_j \mu\_i$, which is the same as R.H.S. 

So the definition of inner-product I wrote seems to work here. I am still confused.

Final Edit:
I feel stupid writing this, but i have tracked down my error to a faulty piece of Matlab code: apparently the command Eigs (based on Arnoldi algorithm to find top k eigenvalues) gives wrong results (such as complex eigenvalues/vectors). I found the whole spectrum using the other command "Eig" and it turns out R has real spectrum, and orthogonal eigenvectors w.r.t to the inner-product I described.