# v_besratinv0

## PURPOSE

V_BESRATINV0 Inverse function of the Modified Bessel Ratio I1(k)/I0(k)

## SYNOPSIS

function k=v_besratinv0(r)

## DESCRIPTION

     V_BESRATINV0 Inverse function of the Modified Bessel Ratio I1(k)/I0(k)

      Inputs: r    Input argument in range [0,1] (scalar or matrix)

     Outputs: k    Ouput satisfying r=v_besselratio(k,0)=I1(k)/I0(k) (same shape as r)

     This routine is a translation of VKAPPA(R) from [2] which is described in [1] and is the
     inverse of r=v_besselratio(k,0). The precision is 8 siginficant digits.
     If the angle q has a von Mises distribution  (a.k.a. Tikhonov distributon) with mean m and
     concentration parameter k, then E(exp(jq)) = I1(k)/I0(k) exp(jm) so k=besratinv0(abs(E(exp(jq)))).

     Refs:
      [1]    G. W. Hill. Evaluation and inversion of the ratios of Modified Bessel functions, I1(x)/I0(x)
           and I1.5(x)/I0.5(x). ACM Transactions on Mathematical Software, 7(2): 199–208, June 1981.
      [2]    G. W. Hill. Algorithm 571: Statistics for von Mises’ and Fisher’s distributions of directions:
           I1(x)/I0(x) and I1.5(x)/I0.5(x) and their inverses [s14].
           ACM Transactions on Mathematical Software, 7(2): 233–238, June 1981. doi: 10.1145/355945.355953.

     Revision History:
     2022/04/21    Original Version

## CROSS-REFERENCE INFORMATION

This function calls:

- v_besselratio V_BESSELRATIO calculate the Bessel function ratio besseli(v+1,x)./besseli(v,x)

This function is called by:

## SOURCE CODE

    0001 function k=v_besratinv0(r)
    0002 % V_BESRATINV0 Inverse function of the Modified Bessel Ratio I1(k)/I0(k)
    0003 %
    0004 %  Inputs: r    Input argument in range [0,1] (scalar or matrix)
    0005 %
    0006 % Outputs: k    Ouput satisfying r=v_besselratio(k,0)=I1(k)/I0(k) (same shape as r)
    0007 %
    0008 % This routine is a translation of VKAPPA(R) from [2] which is described in [1] and is the
    0009 % inverse of r=v_besselratio(k,0). The precision is 8 siginficant digits.
    0010 % If the angle q has a von Mises distribution  (a.k.a. Tikhonov distributon) with mean m and
    0011 % concentration parameter k, then E(exp(jq)) = I1(k)/I0(k) exp(jm) so k=besratinv0(abs(E(exp(jq)))).
    0012 %
    0013 % Refs:
    0014 %  [1]    G. W. Hill. Evaluation and inversion of the ratios of Modified Bessel functions, I1(x)/I0(x)
    0015 %       and I1.5(x)/I0.5(x). ACM Transactions on Mathematical Software, 7(2): 199–208, June 1981.
    0016 %  [2]    G. W. Hill. Algorithm 571: Statistics for von Mises’ and Fisher’s distributions of directions:
    0017 %       I1(x)/I0(x) and I1.5(x)/I0.5(x) and their inverses [s14].
    0018 %       ACM Transactions on Mathematical Software, 7(2): 233–238, June 1981. doi: 10.1145/355945.355953.
    0019 %
    0020 % Revision History:
    0021 % 2022/04/21    Original Version
    0022
    0023 k=zeros(size(r));       % space for outputs
    0024 m1=r<=0.85;             % Use different algorithms for r<=0.85 and r>0.85
    0025 a=r(m1);
    0026 if ~isempty(a)          % r<=0.85 => Use adjusted inverse Taylor Series
    0027     y=2./(1-a);
    0028     x=a.*a;
    0029     s=(((a-5.6076).*a+5.0797).*a-4.6494).*y.*x-1;           % Eqn (7) from [1]
    0030     s=((((s.*x+15).*x+60).*x/360+1).*x-2).*a./(x-1);        % Eqn before (7) from [1]
    0031     m2=a>=0.642;
    0032     b=a(m2);
    0033     if ~isempty(b)
    0034         z=y(m2);
    0035         y=((0.00048.*z-0.1589).*z+0.744).*z-4.2932;         % Two Eqns after (9) in [1]; -dk/dr
    0036         t=s(m2);
    0037         t=(v_besselratio(t,0,9)-b).*y+t;                    % 1st Newton-Raphson Iteration for 0.642<r<=0.85
    0038         m3=b>=0.75;
    0039         c=b(m3);
    0040         if ~isempty(c)
    0041             u=t(m3);
    0042             t(m3)=(v_besselratio(u,0,9)-c).*y(m3)+u;        % 2nd Newton Iteration-Raphson for 0.75<=r<=0.85
    0043         end
    0044         s(m2)=t;        % update the outputs for values needing Newton iterations
    0045     end
    0046     k(m1)=s;            % update the outputs for values 0<=r<=0.85
    0047 end
    0048 a=r(~m1);
    0049 if ~isempty(a) % r>0.85 => Use continued fraction approximation
    0050     y=2./(1-a);
    0051     x=zeros(size(a));
    0052     m2=a>0.95;
    0053     x(m2)=32./(120*a(m2)-131.5+y(m2));                      % Eqn (8b) from [1] for 0.95<=r<=1
    0054     x(~m2)=(-2326*a(~m2)+4317.5526).*a(~m2)-2001.035224;    % Eqn (8c) from [1] for 0.85<r<0.95 (sign-corrected)
    0055     s=(y+1+3./(y-5-12./(y-10-x)))*0.25;
    0056     m2=a<0.95;
    0057     b=a(m2);
    0058     if ~isempty(b)
    0059         z=y(m2);
    0060         y=((0.00048.*z-0.1589).*z+0.744).*z-4.2932;         % Two Eqns after (9) in [1]; -dk/dr
    0061         t=s(m2);
    0062         t=(v_besselratio(t,0,9)-b).*y+t;                    % 1st Newton-Raphson Iteration for 0.85<r<=0.95
    0063         m3=b<=0.875;
    0064         c=b(m3);
    0065         if ~isempty(c)
    0066             u=t(m3);
    0067             t(m3)=(v_besselratio(u,0,9)-c).*y(m3)+u;        % 2nd Newton-Raphson Iteration for 0.85<r<=0.875
    0068         end
    0069         s(m2)=t;        % update the outputs for values needing Newton iterations
    0070     end
    0071     k(~m1)=s;           % update the outputs for values 0.85<r<=1
    0072 end
    0073 if ~nargout && sum(isfinite(k))>1
    0074     plot(r,k);
    0075     ylabel('k');
    0076     xlabel('r = I_1(k)/I_0(k)');
    0077 end
    0078
    0079