Why this determinant is conformally invariant? While I was reading a paper about random analytic function I found a statement that I was not able to prove and after try brute force and search for some references I decided to ask for a help here. The statement is the following: Denote by $\mathbb{D}$ the unit dic on the complex plane, and consider the function $p:\mathbb{D}^n\to\mathbb{C}$ given by 
$$
p(z_1,\dots,z_n)\, =\, 4^{-n}\, \det\left(\frac{(1-|z_i|^2)(1-|z_j|^2)}{|1-z_i\overline{z}_j|^2}\right)_{i,j=1}^{n}\,
$$ 
then for any fixed $u\in\mathbb{D}$ the function $p$ is invariant for the Möbius transformation
$$
\phi(u,z)=\frac{u+z}{1+\overline{u}z},
$$
in the sense that
$$
p(z_1,\dots,z_n)\, = p(\phi(u,z_1),\ldots,\phi(u,z_n)).
$$
I appreciate any reference or hint to prove this fact.
Edition: The denominator was modified as point out for David and Anon. 
 A: I suspect you have a typo, and the denominator should be $|1-z_i \overline{z_j}|^2$. At least, I will show that the quantity is Möbius invariant with this modification. It is possible, of course, that both are invariant.
The determinant is a red herring; each individual entry in the matrix is invariant under Möbius transformations. I begin by quoting a fact from hyperbolic geometry: Define
$$\delta(u,v) = \frac{|u-v|^2}{(1-|u|^2)(1-|v|^2)}$$
for $u$ and $v$ in $\mathbb{D}$. Then this quantity is invariant under Möbius transformations. I can't give an intuition for this fact, but its easy to give a reference.
So
$$1+\delta(u,v) = 1+\frac{(u-v)(\overline{u} - \overline{v})}{(1-u \overline{u}) (1-v \overline{v})} = \frac{\left( 1-u \overline{u} - v \overline{v} + u \overline{u} v \overline{v} \right) + \left( u \overline{u} - u \overline{v} - \overline{u} v + v \overline{v} \right)}{(1-u \overline{u}) (1-v \overline{v})}$$
$$=\frac{1- u \overline{v} - \overline{u} v + u \overline{u} v \overline{v}}{(1-u \overline{u}) (1-v \overline{v})} = \frac{(1-u \overline{v})(1-\overline{u} v)}{(1-|u|^2)(1-|v|^2)} = \frac{|1-u \overline{v}|^2}{(1-|u|^2)(1-|v|^2)}$$
is invariant under Möbius tranformations. 
Replacing $u$ and $v$ by $z_i$ and $z_j$, this is the reciprocal of the $(i,j)$ entry in your matrix. So every element of your matrix is Möbius invariant, as claimed.
UPDATE: I just realized a nice way to express anon's solution below. Let $\sigma$ denote Schwarz reflection in the boundary of $\mathbb{D}$. Since Schwarz reflection is a conformally invariant operation, if $\phi$ is a Möbius 
transformation of $\mathbb{P}^1$ preserving $\mathbb{D}$, then $\phi(\sigma(z)) = \sigma(\phi(z))$. Explicitly, $\sigma(z) = \overline{z}^{-1}$.
Then the $(i,j)$ entry in your matrix is the cross ratio $(z_i, z_j; \sigma(z_i), \sigma(z_j))$, by a computation very similar to anon's.
A: More pointedly, each entry in the matrix is a cross-ratio:
$$\frac{(1-|z|^2)(1-|w|^2)}{|1-z\bar{w}|^2}=\frac{(1-z\bar{z})(1-w\bar{w})}{(1-\bar{z}w)(1-z\bar{w})}=(z^{-1},w^{-1};\bar{z},\bar{w}).$$
This can be seen by dividing both numerator and denominator by $zw$ appropriately. Cross-ratios are invariant under Möbius transformations, in that sense that if $f$ is one, then
$$(a,b;c,d)=(f(a),f(b);f(c),f(d)).$$
EDIT: Ah, the last ingredient is provided by Theo Buehler. Since $\phi(u,z)^{-1}=\phi(\bar{u},z)$ and similarly we have $\overline{\phi(u,z)}=\phi(\bar{u},z)$, we can say that $\phi(u,\cdot)$ acting on both $z$ and $w$ is equivalent to the Möbius transformation $\phi(\bar{u},\cdot)$ acting on each individual component the $4$-tuple $(z^{-1},w^{-1},\bar{z},\bar{w})$, hence the cross-ratio is preserved under it.
