Our aim is to understand the stress-dependent seismic anisotropy of the overburden shale in an oil field in the North West Shelf of Western Australia. We analyze data from measurements of ultrasonic P-wave velocities in 132 directions for confining pressures of 0.1–400 MPa on a spherical shale sample. First, we find the orientation of the symmetry axis, assuming that the sample is transversely isotropic, and then transform the ray velocities to the symmetry axis coordinates. We use two parameterizations of the phase velocity; one, in terms of the Thomsen anisotropy parameters α, β, ɛ, δ as the main approach, and the other in terms of α, β, η, δ. We invert the ray velocities to estimate the anisotropy parameters α, ɛ, δ, and η using a very fast simulated reannealing algorithm. Both approaches result in the same estimation for the anisotropy parameters but with different uncertainties. The main approach is robust but produces higher uncertainties, in particular for η, whereas the alternative approach is unstable but gives lower uncertainties. These approaches are used to find the anisotropy parameters for the different confining pressures. The dependency of P-wave velocity, α, on pressure has exponential and linear components, which can be contributed to the compliant and stiff porosities. The exponential dependence at lower pressures up to 100 MPa corresponds to the closure of compliant pores and microcracks, whereas the linear dependence at higher pressures corresponds to contraction of the stiff pores. The anisotropy parameters ɛ and δ are quite large at lower pressures but decrease exponentially with pressure. For lower pressures up to 10 MPa, δ always is larger than ɛ; this trend is reversed for higher pressures. Despite the hydrostatic pressure, the symmetry axis orientation changes noticeably, in particular at lower pressures.