R Under development (unstable) (2023-11-25 r85635 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(HyperbolicDist) > > ### Create vector of nu values and of x values > nus <- c(0:5, 10, 20) > x <- seq(1, 4, length.out = 11) > ### Specify the difference between the numerator and denominator orders > k <- 3 > > ### Matrix for results of calculations > results <- matrix(nrow = length(nus), ncol = length(x)) > for (i in 1:length(nus)){ + for (j in 1:length(x)) { + raw <- besselK(x[j], nus[i] + k)/besselK(x[j], nus[i]) + scaled <- besselK(x[j], nus[i] + k, expon.scaled = TRUE)/ + besselK(x[j], nus[i], expon.scaled = TRUE) + results[i,j] <- raw/scaled + } + } > > ### Should yield 1 in every case > results [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 1 1 1 1 1 1 1 1 1 1 1 [2,] 1 1 1 1 1 1 1 1 1 1 1 [3,] 1 1 1 1 1 1 1 1 1 1 1 [4,] 1 1 1 1 1 1 1 1 1 1 1 [5,] 1 1 1 1 1 1 1 1 1 1 1 [6,] 1 1 1 1 1 1 1 1 1 1 1 [7,] 1 1 1 1 1 1 1 1 1 1 1 [8,] 1 1 1 1 1 1 1 1 1 1 1 > max(abs(results - 1)) [1] 4.440892e-16 > > ### Try a negative value for the difference of the orders > k <- -3 > > results <- matrix(nrow = length(nus), ncol = length(x)) > for (i in 1:length(nus)){ + for (j in 1:length(x)) { + raw <- besselK(x[j], nus[i] + k)/besselK(x[j], nus[i]) + scaled <- besselK(x[j], nus[i] + k, expon.scaled = TRUE)/ + besselK(x[j], nus[i], expon.scaled = TRUE) + results[i,j] <- raw/scaled + } + } > > results [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 1 1 1 1 1 1 1 1 1 1 1 [2,] 1 1 1 1 1 1 1 1 1 1 1 [3,] 1 1 1 1 1 1 1 1 1 1 1 [4,] 1 1 1 1 1 1 1 1 1 1 1 [5,] 1 1 1 1 1 1 1 1 1 1 1 [6,] 1 1 1 1 1 1 1 1 1 1 1 [7,] 1 1 1 1 1 1 1 1 1 1 1 [8,] 1 1 1 1 1 1 1 1 1 1 1 > max(abs(results - 1)) [1] 4.440892e-16 > > ### Now use besselRatio function > nus <- c(0:5, 10, 20) > x <- seq(1, 4, length.out = 11) > k <- 3 > > raw <- matrix(nrow = length(nus), ncol = length(x)) > scaled <- matrix(nrow = length(nus), ncol = length(x)) > compare <- matrix(nrow = length(nus), ncol = length(x)) > > for (i in 1:length(nus)){ + for (j in 1:length(x)) { + raw[i,j] <- besselRatio(x[j], nus[i], + orderDiff = k) + scaled[i,j] <- besselRatio(x[j], nus[i], + orderDiff = k, useExpScaled = 1) + compare[i,j] <- raw[i,j]/scaled[i,j] + } + } > raw [,1] [,2] [,3] [,4] [,5] [,6] [1,] 16.86663 10.75385 7.781137 6.090468 5.024655 4.302128 [2,] 73.48710 39.35525 24.822473 17.380739 13.074196 10.355893 [3,] 222.15162 109.48831 64.098882 41.995961 29.773284 22.368514 [4,] 514.53360 244.54113 137.964988 87.165108 59.673423 43.368998 [5,] 999.42586 467.03943 258.597043 160.171259 107.450155 76.520698 [6,] 1724.70940 798.93325 437.859195 268.144833 177.715991 124.972350 [7,] 10633.28765 4862.90579 2623.886703 1578.095601 1024.973207 705.066531 [8,] 74052.62238 33747.89311 18129.755109 10846.870001 7002.415905 4783.910099 [,7] [,8] [,9] [,10] [,11] [1,] 3.785224 3.399798 3.102843 2.867909 2.677938 [2,] 8.525505 7.230002 6.276125 5.550945 4.984885 [3,] 17.564024 14.275849 11.927487 10.191042 8.869535 [4,] 33.006241 26.053230 21.180509 17.641581 14.993686 [5,] 57.081406 44.185346 35.249374 28.831590 24.082062 [6,] 92.029179 70.317805 55.375514 44.718344 36.886643 [7,] 507.122749 377.998297 290.107405 228.152598 183.181231 [8,] 3414.690278 2524.043094 1919.703872 1495.154973 1188.123170 > scaled [,1] [,2] [,3] [,4] [,5] [,6] [1,] 16.86663 10.75385 7.781137 6.090468 5.024655 4.302128 [2,] 73.48710 39.35525 24.822473 17.380739 13.074196 10.355893 [3,] 222.15162 109.48831 64.098882 41.995961 29.773284 22.368514 [4,] 514.53360 244.54113 137.964988 87.165108 59.673423 43.368998 [5,] 999.42586 467.03943 258.597043 160.171259 107.450155 76.520698 [6,] 1724.70940 798.93325 437.859195 268.144833 177.715991 124.972350 [7,] 10633.28765 4862.90579 2623.886703 1578.095601 1024.973207 705.066531 [8,] 74052.62238 33747.89311 18129.755109 10846.870001 7002.415905 4783.910099 [,7] [,8] [,9] [,10] [,11] [1,] 3.785224 3.399798 3.102843 2.867909 2.677938 [2,] 8.525505 7.230002 6.276125 5.550945 4.984885 [3,] 17.564024 14.275849 11.927487 10.191042 8.869535 [4,] 33.006241 26.053230 21.180509 17.641581 14.993686 [5,] 57.081406 44.185346 35.249374 28.831590 24.082062 [6,] 92.029179 70.317805 55.375514 44.718344 36.886643 [7,] 507.122749 377.998297 290.107405 228.152598 183.181231 [8,] 3414.690278 2524.043094 1919.703872 1495.154973 1188.123170 > compare [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 1 1 1 1 1 1 1 1 1 1 1 [2,] 1 1 1 1 1 1 1 1 1 1 1 [3,] 1 1 1 1 1 1 1 1 1 1 1 [4,] 1 1 1 1 1 1 1 1 1 1 1 [5,] 1 1 1 1 1 1 1 1 1 1 1 [6,] 1 1 1 1 1 1 1 1 1 1 1 [7,] 1 1 1 1 1 1 1 1 1 1 1 [8,] 1 1 1 1 1 1 1 1 1 1 1 > > proc.time() user system elapsed 0.14 0.06 0.18