Stream: helpdesk (published)

Topic: Inverse normal transformation

Kevin Bonham (Aug 03 2021 at 15:32):

Are there any packages that implement INT? This is defined as

image.png

where $\Phi^{-1}$ is the probit function, $r_{i}$ is the rank of observation $i$, and $N$ is the number of observations.

Kevin Bonham (Aug 03 2021 at 15:34):

I did a bit of googling for this or for probit, but all of the results are about probit link functions in the context of linear modeling for example

Júlio Hoffimann (Aug 03 2021 at 15:46):

StatsFuns.jl has normcdf and norminvcdf

Júlio Hoffimann (Aug 03 2021 at 15:52):

Alternatively, take a look at the erf functions from SpecialFunctions.jl:

https://specialfunctions.juliamath.org/stable/functions_overview/#Error-Functions,-Dawson%E2%80%99s-and-Fresnel-Integrals-1

Júlio Hoffimann (Aug 03 2021 at 15:52):

The Gaussian cdf can be expressed and properly implemented with the erf function

Chad Scherrer (Aug 03 2021 at 16:16):

There might be a better way to get the ranks (is there?), but this works:

julia> x = [3,5,2,7]
4-element Vector{Int64}:
 3
 5
 2
 7

julia> (1:4)[invperm(sortperm(x))]
4-element Vector{Int64}:
 2
 3
 1
 4

Kevin Bonham (Aug 03 2021 at 16:32):

@Chad Scherrer That's the way I'd do it :laughing:... well, skipping using the range

julia> sortperm(x) |> invperm
4-element Vector{Int64}:
 2
 3
 1
 4

Kevin Bonham (Aug 03 2021 at 16:52):

Cool.

using Distributions
using CairoMakie
using StatsFuns
using StatsBase

function invnormaltransform(v; c=3/8)
    (μ, σ) = mean_and_var(v)
    rank = invperm(sortperm(v))
    N = length(v)
    return [norminvcdf(μ, σ, (x - c) / (N - 2c + 1)) for x in rank]
end

x = rand(Beta(3, 1), 1000); y = invnormaltransform(x);

hist(x)
hist!(y)
current_figure()

image.png

Last updated: Oct 02 2023 at 04:34 UTC