Inverse of this function in Matlab

Question

im trying to implement the following functions in Matlab:

$$\Phi = \text{erf}(\frac{h^{-1}(\rho)}{\sqrt{2}}), \quad \rho > 0$$ where $$h(\tau) = \sqrt{\frac{2}{\pi}} \frac{\exp{-\frac{\tau^2}{2}}}{\tau} + \text{erf}(\frac{\tau}{\sqrt{2}}) - 1.$$erf here is the error function which is also provided by Matlab. The difficulty I'm having with this is implementing the inverse of $h$ which is needed for the evaluation of $\Phi$. What I'm doing so far is the following:

function y = Phi(r)

h = @(tau) sqrt(2/pi) * exp(-(tau^2 / 2))/tau + erf(tau / sqrt(2)) - 1;

x = fzero(@(x)(h(x) - r), x0);

y = erf(x / sqrt(2));

end

What I'm having trouble with now is choosing a right estimate for $x0$. I have tried values close to $0$, but for larger values of $\tau$ the program tells me that

Exiting fzero: aborting search for an interval containing a sign change
because no sign change is detected during search.
Function may not have a root.

Also, $\Phi$ should obey $\lim_{\rho \rightarrow 0} \Phi(\rho) = 1$, but I'm getting values very close to zero, $10^{-16}$, if the argument is close to zero itself. Then I've tried intervalls $[a,b]$ with $a$ very close to zero and $b$ arbitraty. That worked better, but for large values I now get:

Error using fzero (line 290)
The function values at the interval endpoints must differ in sign.

Could anybody tell me how to better implement the function $\Phi$? Thanks in advance!

Answer 1

I'm not sure if this is a great solution, but it seems to work.

The idea is that we want to find the root of a function $h$ which has a domain the extends to infinity. So the idea is just to scale it into a finite interval. In this case, the scaling function I chose is $g(t)=\tan(\pi t/2)$, which maps the interval $[0,1)$ into $[0,\infty)$.

We're looking for $\tau_0$ such that $h(\tau_0)=\rho$ for some $\rho$. So equivalently, we can find $t_0$ such that $h(g(t_0))=\rho$, then $\tau_0$ is given by $\tau_0=g(t_0)$.

Here's the code:

function y = Phi(r)

h = @(tau) sqrt(2/pi) * exp(-(tau.^2)/2) * ( 1./tau ) + erf( tau / sqrt(2) ) - 1;

g = @(t) tan( pi * t / 2 );

t0 = fzero( @(t) ( h( g( t ) ) - r ), 0.5 );

x0 = g( t0 );

y = erf( x0 / sqrt(2) );

end

and here's a plot of the function $\Phi$ on the interval $(0,0.1]$ (I'm able to evaluate $\Phi$ at $10^{-4}$--didn't try any further).