# How can I estimate theta for the inverse hyperbolic sine transformation?

I would like help with R code to estimate theta for the Inverse Hyperbolic Sine Transformation. This transformation is useful to transform skewed data that contain negative values or zeros.

There are a couple of related posts that discuss the IHS and suggest there is a maximum likelihood approach to estimating theta but I can't work out how to apply it. These related posts are:

Please see example code below. I've made skewed data by using the inverse of the IHS. Now, given the skewed data, and no prior knowledge of theta, how can I work out what theta should be? I would be most grateful for R code to undertake this analysis.

# Define the IHS transformation and its inverse
IHS <- function(x, theta){  # Inverse IHS transformation
(1/theta)*asinh(theta * x)
}

Inv.IHS <- function(x, theta){  # IHS transformation
(1/theta)*sinh(theta * x)
}

set.seed(1)
# generate some normal data
x <- rnorm(1000)
hist(x, breaks='FD')

# skew it by applying the Inverse of the IHS transformation
xt <- Inv.IHS(x, theta=2)
hist(xt, breaks='FD') # yep this is skewed.  How could we estimate theta?


• If you want to estimate it, maximum likelihood estimation, perhaps? Dec 8, 2013 at 9:52

After a couple more days of thinking about the problem, I have two tentative answers.

1. Select theta so that the transformed data is close to normal as measured by goodness of fit. For example choose theta to maximize the p-value of the Shapiro-Wik test.

set.seed(1)
x <- rnorm(1000)
xt <- Inv.IHS(x, theta=2)

Shapiro.test.pvalue <- function(theta, x){
x <- IHS(x, theta)
shapiro.test(x)\$p.value
}

optimise(Shapiro.test.pvalue, lower=0.001, upper=50, x=xt, maximum=TRUE)  # 2.069838

2. Maximum likelihood estimation of theta.

Looking at the paper by Burbidge et al. I think the likelihood function for a single variable can be expressed as follows.

IHS.loglik <- function(theta,x){

IHS <- function(x, theta){  # function to IHS transform
asinh(theta * x)/theta
}

n <- length(x)
xt <- IHS(x, theta)

log.lik <- -n*log(sum((xt - mean(xt))^2))- sum(log(1+theta^2*x^2))
return(log.lik)
}

# try this on our data
optimise(IHS.loglik, lower=0.001, upper=50, x=xt, maximum=TRUE) # 2.0407


In both cases we get close to the expected theta, which is encouraging. But when I try the maximum likelihood approach on my real data it doesn't seem to give reasonable answers.

• +1 (on both Q and A). ... "when I try the maximum likelihood approach on my real data it doesn't seem to give reasonable answers" - if the other approach does give reasonable answers, it might suggest a convergence problem with the maximization routine. What happens if you start it at the other estimate? Dec 9, 2013 at 22:32
• Unfortunately the shapiro version just works for small samples. It would be interesting to see an stable answer using ML. Nov 6, 2014 at 20:20