Okamoto has obtained a sequence of τ-functions for the PVI system expressed as a double Wronskian determinant based on a solution of the Gauss hypergeometric equation. Starting with integral solutions of the Gauss hypergeometric equation, we show that the determinant can be re-expressed as multidimensional integrals, and these in turn can be identified with averages over the eigenvalue probability density function for the Jacobi unitary ensemble (JUE), and the Cauchy unitary ensemble (CyUE) (the latter being equivalent to the circular Jacobi unitary ensemble (cJUE)). Hence these averages, which depend on four continuous parameters and the discrete parameter N, can be characterised as the solution of the second order second degree equation satisfied by the Hamiltonian in the PVI theory. We show that the Hamiltonian also satisfies an equation related to the discrete PV equation, thus providing an alternative characterisation in terms of a difference equation. In the case of the cJUE, the spectrum singularity scaled limit is considered, and the evaluation of a certain four parameter average is given in terms of the general PV transcendent in σ form. Applications are given to the evaluation of the spacing distribution for the circular unitary ensemble (CUE) and its scaled counterpart, giving formulas more succinct than those known previously; to expressions for the hard edge gap probability in the scaled Laguerre orthogonal ensemble (LOE) (parameter a a non-negative integer) and Laguerre symplectic ensemble (LSE) (parameter a an even non-negative integer) as finite dimensional combinatorial integrals over the symplectic and orthogonal groups respectively; to the evaluation of the cumulative distribution function for the last passage time in certain models of directed percolation; to the τ-function evaluation of the largest eigenvalue in the finite LOE and LSE with parameter a = 0; and to the characterisation of the diagonal-diagonal spin-spin correlation in the two-dimensional Ising model.
