v_besratinv0

 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

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