Potential for the "Exponentially scaled modified Bessel function of the first kind"

Developers
1

To implement periodic HSGPs the so called Exponentially scaled modified Bessel function of the first kind is needed.

The above linked implementation by @avehtari et al. appears to have issues for small length scales, as the modified Bessel function of the first kind explodes rather quickly (note the log y-scaling) even though the exponentially scaled modified Bessel function does not. Stan already has log_modified_bessel_first_kind, but that just crashed my notebook.

Edit: Also, this is my current (hopefully C1 and hopefully capturing the asymptotics) workaround for up to 32x2 basis functions:


vector approx_sqrt_ive(int[] idxs, real z){
  real log_threshold = 6.214608098422191;
  real log_z = log(z);
  if(log_z < log_threshold){
    return sqrt(to_vector(modified_bessel_first_kind(idxs, z)))
    * exp(-.5*z);
  }
  real coeffs[32,4] = {{-0.45947204529166896,-0.24999981244346867,14033.07814847007,-2688.6129885357695},{-0.4594831553850529,-0.24999906255493107,2806.210005067361,-537.65317061544},{-0.4595016763612292,-0.24999781244003452,1202.7199215199269,-230.4326271206091},{-0.45952760224842537,-0.24999606253103046,668.1752966301113,-128.0177771572948},{-0.4595609379212817,-0.2499938124750806,425.2008557602532,-81.4656412715292},{-0.4596016798582645,-0.24999106252708206,294.3683082132673,-56.39907948567251},{-0.45964983069099263,-0.24998781249655191,215.86895335995402,-41.3591716285243},{-0.4597053898175463,-0.24998406242705853,165.07451899633904,-31.62733985164608},{-0.45976835434106,-0.24997981252828408,130.3213719642655,-24.968881742142795},{-0.45983872892599287,-0.2499750624626037,105.49731203818817,-20.21277337668723},{-0.4599165103097116,-0.24996981246617475,87.14910105568313,-16.697390146538893},{-0.4600016999748666,-0.24996406243167932,73.20447550707883,-14.025701324099892},{-0.46009429617038755,-0.24995781248586416,62.35873118853225,-11.947733731596811},{-0.4601943009518963,-0.24995106247993837,53.75678048399213,-10.299664918400788},{-0.4603017150187685,-0.2499438123632795,46.8198012894012,-8.970590083186076},{-0.46041653348954004,-0.24993606248921962,41.144075337605784,-7.883162642561775},{-0.4605387618609318,-0.24992781245989287,36.44134922540376,-6.982154281434953},{-0.46066839864072806,-0.24991906238330944,32.501220528640616,-6.227254367224537},{-0.46080544238838517,-0.2499098123637393,29.167265336840575,-5.588492831685634},{-0.4609498917150048,-0.2499000625017143,26.3212272694486,-5.043212156126089},{-0.46110175326367164,-0.24988981231639187,23.872271287360284,-4.5740105266822795},{-0.4612610191770723,-0.24987906237650342,21.749859704481267,-4.167371956344471},{-0.461427692915819,-0.24986781243156142,19.898396027954117,-3.812644973487975},{-0.4616017765841689,-0.2498560623292547,18.273642040999555,-3.5013539291970077},{-0.4617832667046158,-0.2498438123212935,16.84003483870306,-3.2266852385565095},{-0.46197216454590384,-0.24983106231584307,15.568723480491911,-2.983111227476678},{-0.4621684695506203,-0.24981781235325032,14.43610078378073,-2.766108951527705},{-0.46237218356055676,-0.24980406230020188,13.422701522222953,-2.5719490172442057},{-0.4625833047132035,-0.24978981229151073,12.512359936685524,-2.3975341807683073},{-0.46280183448430456,-0.24977506222035895,11.691557397963575,-2.2402744362551323},{-0.4630277692363669,-0.24975981235003702,10.948931976480697,-2.0979927271367704},{-0.46326111571363393,-0.24974406219238043,10.274854131199614,-1.9688442992881874}};
  vector[size(idxs)] rv;
  for(idx in idxs){
    rv[idx] = coeffs[idx, 1] + coeffs[idx, 2] * log_z + 1 / (
      coeffs[idx, 3] + coeffs[idx, 4] * log_z
    );
  }
  return exp(rv);
}
5 Likes