\(\phi^{-1}(z) = \Phi^{-1}(F(z))\), eq (5.60): "matching the percentiles of the data to those of a standard Gaussian distribution". The \((i-1/2)/n\) plotting position is used so no observation maps to an infinite score, which plain \(i/n\) would do to the largest value.
Examples
V <- c(1, 2, 3, 4, 5, 6, 7, 8)
morie_normal_scores(V)
#> [1] -1.5341205 -0.8871466 -0.4887764 -0.1573107 0.1573107 0.4887764 0.8871466
#> [8] 1.5341205