diff --git a/Cargo.toml b/Cargo.toml index e291172..44fa244 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -6,8 +6,8 @@ edition = "2018" license = "MIT OR Apache-2.0" [dependencies] -lapacke = "0.2" -ndarray = { version = "0.14", features = ["approx"] } +lapacke = "0.5" +ndarray = { version = "0.15", features = ["approx"] } [dev-dependencies] approx = "0.4" diff --git a/src/kernel/affine.rs b/src/kernel/affine.rs index c904932..94f1c6c 100644 --- a/src/kernel/affine.rs +++ b/src/kernel/affine.rs @@ -1,5 +1,3 @@ -use std::iter::FromIterator; - use ndarray::prelude::*; use super::Kernel; @@ -36,7 +34,7 @@ impl Kernel for Affine { } fn k_diag(&self, ts: &[f64]) -> Array1 { - let r = Array1::from_iter(ts.iter().map(|v| (v - self.t0).powi(2))); + let r: Array1<_> = ts.iter().map(|v| (v - self.t0).powi(2)).collect(); (r * self.var_slope) + self.var_offset } diff --git a/src/kernel/wiener.rs b/src/kernel/wiener.rs index ada3721..2d84abc 100644 --- a/src/kernel/wiener.rs +++ b/src/kernel/wiener.rs @@ -1,5 +1,3 @@ -use std::iter::FromIterator; - use ndarray::prelude::*; use super::Kernel; @@ -32,7 +30,7 @@ impl Kernel for Wiener { } fn k_diag(&self, ts: &[f64]) -> Array1 { - (Array1::from_iter(ts.iter().map(|v| v - self.t0)) * self.var) + self.var_t0 + (ts.iter().map(|v| v - self.t0).collect::>() * self.var) + self.var_t0 } fn order(&self) -> usize { diff --git a/src/math.rs b/src/math.rs index 65105d2..9abaaec 100644 --- a/src/math.rs +++ b/src/math.rs @@ -13,310 +13,310 @@ pub fn erfc(x: f64) -> f64 { /// Polynomial coefficients for a numerator of `erf_impl` /// in the interval [1e-10, 0.5]. const ERF_IMPL_AN: &[f64] = &[ - 0.00337916709551257388990745, - -0.00073695653048167948530905, - -0.374732337392919607868241, - 0.0817442448733587196071743, - -0.0421089319936548595203468, - 0.0070165709512095756344528, - -0.00495091255982435110337458, - 0.000871646599037922480317225, + 0.003_379_167_095_512_573_7, + -0.000_736_956_530_481_679_5, + -0.374_732_337_392_919_6, + 0.081_744_244_873_358_73, + -0.042_108_931_993_654_86, + 0.007_016_570_951_209_575, + -0.004_950_912_559_824_351, + 0.000_871_646_599_037_922_5, ]; /// Polynomial coefficients for a denominator of `erf_impl` /// in the interval [1e-10, 0.5] const ERF_IMPL_AD: &[f64] = &[ 1.0, - -0.218088218087924645390535, - 0.412542972725442099083918, - -0.0841891147873106755410271, - 0.0655338856400241519690695, - -0.0120019604454941768171266, - 0.00408165558926174048329689, - -0.000615900721557769691924509, + -0.218_088_218_087_924_64, + 0.412_542_972_725_442_1, + -0.084_189_114_787_310_67, + 0.065_533_885_640_024_16, + -0.012_001_960_445_494_177, + 0.004_081_655_589_261_74, + -0.000_615_900_721_557_769_7, ]; /// Polynomial coefficients for a numerator in `erf_impl` /// in the interval [0.5, 0.75]. const ERF_IMPL_BN: &[f64] = &[ - -0.0361790390718262471360258, - 0.292251883444882683221149, - 0.281447041797604512774415, - 0.125610208862766947294894, - 0.0274135028268930549240776, - 0.00250839672168065762786937, + -0.036_179_039_071_826_25, + 0.292_251_883_444_882_7, + 0.281_447_041_797_604_5, + 0.125_610_208_862_766_94, + 0.027_413_502_826_893_053, + 0.002_508_396_721_680_657_5, ]; /// Polynomial coefficients for a denominator in `erf_impl` /// in the interval [0.5, 0.75]. const ERF_IMPL_BD: &[f64] = &[ 1.0, - 1.8545005897903486499845, - 1.43575803037831418074962, - 0.582827658753036572454135, - 0.124810476932949746447682, - 0.0113724176546353285778481, + 1.854_500_589_790_348_6, + 1.435_758_030_378_314_2, + 0.582_827_658_753_036_5, + 0.124_810_476_932_949_75, + 0.011_372_417_654_635_328, ]; /// Polynomial coefficients for a numerator in `erf_impl` /// in the interval [0.75, 1.25]. const ERF_IMPL_CN: &[f64] = &[ - -0.0397876892611136856954425, - 0.153165212467878293257683, - 0.191260295600936245503129, - 0.10276327061989304213645, - 0.029637090615738836726027, - 0.0046093486780275489468812, - 0.000307607820348680180548455, + -0.039_787_689_261_113_69, + 0.153_165_212_467_878_3, + 0.191_260_295_600_936_24, + 0.102_763_270_619_893_04, + 0.029_637_090_615_738_836, + 0.004_609_348_678_027_549, + 0.000_307_607_820_348_680_2, ]; /// Polynomial coefficients for a denominator in `erf_impl` /// in the interval [0.75, 1.25]. const ERF_IMPL_CD: &[f64] = &[ 1.0, - 1.95520072987627704987886, - 1.64762317199384860109595, - 0.768238607022126250082483, - 0.209793185936509782784315, - 0.0319569316899913392596356, - 0.00213363160895785378615014, + 1.955_200_729_876_277, + 1.647_623_171_993_848_6, + 0.768_238_607_022_126_2, + 0.209_793_185_936_509_78, + 0.031_956_931_689_991_336, + 0.002_133_631_608_957_853_7, ]; /// Polynomial coefficients for a numerator in `erf_impl` /// in the interval [1.25, 2.25]. const ERF_IMPL_DN: &[f64] = &[ - -0.0300838560557949717328341, - 0.0538578829844454508530552, - 0.0726211541651914182692959, - 0.0367628469888049348429018, - 0.00964629015572527529605267, - 0.00133453480075291076745275, - 0.778087599782504251917881e-4, + -0.030_083_856_055_794_972, + 0.053_857_882_984_445_45, + 0.072_621_154_165_191_42, + 0.036_762_846_988_804_936, + 0.009_646_290_155_725_275, + 0.001_334_534_800_752_910_7, + 7.780_875_997_825_043e-5, ]; /// Polynomial coefficients for a denominator in `erf_impl` /// in the interval [1.25, 2.25]. const ERF_IMPL_DD: &[f64] = &[ 1.0, - 1.75967098147167528287343, - 1.32883571437961120556307, - 0.552528596508757581287907, - 0.133793056941332861912279, - 0.0179509645176280768640766, - 0.00104712440019937356634038, - -0.106640381820357337177643e-7, + 1.759_670_981_471_675_3, + 1.328_835_714_379_611_2, + 0.552_528_596_508_757_6, + 0.133_793_056_941_332_87, + 0.017_950_964_517_628_076, + 0.001_047_124_400_199_373_6, + -1.066_403_818_203_573_4e-8, ]; /// Polynomial coefficients for a numerator in `erf_impl` /// in the interval [2.25, 3.5]. const ERF_IMPL_EN: &[f64] = &[ - -0.0117907570137227847827732, - 0.014262132090538809896674, - 0.0202234435902960820020765, - 0.00930668299990432009042239, - 0.00213357802422065994322516, - 0.00025022987386460102395382, - 0.120534912219588189822126e-4, + -0.011_790_757_013_722_784, + 0.014_262_132_090_538_81, + 0.020_223_443_590_296_084, + 0.009_306_682_999_904_321, + 0.002_133_578_024_220_66, + 0.000_250_229_873_864_601_05, + 1.205_349_122_195_881_9e-5, ]; /// Polynomial coefficients for a denominator in `erf_impl` /// in the interval [2.25, 3.5]. const ERF_IMPL_ED: &[f64] = &[ 1.0, - 1.50376225203620482047419, - 0.965397786204462896346934, - 0.339265230476796681555511, - 0.0689740649541569716897427, - 0.00771060262491768307365526, - 0.000371421101531069302990367, + 1.503_762_252_036_204_8, + 0.965_397_786_204_462_9, + 0.339_265_230_476_796_7, + 0.068_974_064_954_156_98, + 0.007_710_602_624_917_683, + 0.000_371_421_101_531_069_3, ]; /// Polynomial coefficients for a numerator in `erf_impl` /// in the interval [3.5, 5.25]. const ERF_IMPL_FN: &[f64] = &[ - -0.00546954795538729307482955, - 0.00404190278731707110245394, - 0.0054963369553161170521356, - 0.00212616472603945399437862, - 0.000394984014495083900689956, - 0.365565477064442377259271e-4, - 0.135485897109932323253786e-5, + -0.005_469_547_955_387_293, + 0.004_041_902_787_317_071, + 0.005_496_336_955_316_117, + 0.002_126_164_726_039_454, + 0.000_394_984_014_495_083_9, + 3.655_654_770_644_424e-5, + 1.354_858_971_099_323_2e-6, ]; /// Polynomial coefficients for a denominator in `erf_impl` /// in the interval [3.5, 5.25]. const ERF_IMPL_FD: &[f64] = &[ 1.0, - 1.21019697773630784832251, - 0.620914668221143886601045, - 0.173038430661142762569515, - 0.0276550813773432047594539, - 0.00240625974424309709745382, - 0.891811817251336577241006e-4, - -0.465528836283382684461025e-11, + 1.210_196_977_736_307_7, + 0.620_914_668_221_143_9, + 0.173_038_430_661_142_77, + 0.027_655_081_377_343_203, + 0.002_406_259_744_243_097_3, + 8.918_118_172_513_365e-5, + -4.655_288_362_833_827e-12, ]; /// Polynomial coefficients for a numerator in `erf_impl` /// in the interval [5.25, 8]. const ERF_IMPL_GN: &[f64] = &[ - -0.00270722535905778347999196, - 0.0013187563425029400461378, - 0.00119925933261002333923989, - 0.00027849619811344664248235, - 0.267822988218331849989363e-4, - 0.923043672315028197865066e-6, + -0.002_707_225_359_057_783_7, + 0.001_318_756_342_502_94, + 0.001_199_259_332_610_023_3, + 0.000_278_496_198_113_446_64, + 2.678_229_882_183_318_6e-5, + 9.230_436_723_150_282e-7, ]; /// Polynomial coefficients for a denominator in `erf_impl` /// in the interval [5.25, 8]. const ERF_IMPL_GD: &[f64] = &[ 1.0, - 0.814632808543141591118279, - 0.268901665856299542168425, - 0.0449877216103041118694989, - 0.00381759663320248459168994, - 0.000131571897888596914350697, - 0.404815359675764138445257e-11, + 0.814_632_808_543_141_6, + 0.268_901_665_856_299_54, + 0.044_987_721_610_304_114, + 0.003_817_596_633_202_484_7, + 0.000_131_571_897_888_596_92, + 4.048_153_596_757_641e-12, ]; /// Polynomial coefficients for a numerator in `erf_impl` /// in the interval [8, 11.5]. const ERF_IMPL_HN: &[f64] = &[ - -0.00109946720691742196814323, - 0.000406425442750422675169153, - 0.000274499489416900707787024, - 0.465293770646659383436343e-4, - 0.320955425395767463401993e-5, - 0.778286018145020892261936e-7, + -0.001_099_467_206_917_422, + 0.000_406_425_442_750_422_67, + 0.000_274_499_489_416_900_7, + 4.652_937_706_466_594e-5, + 3.209_554_253_957_674_6e-6, + 7.782_860_181_450_209e-8, ]; /// Polynomial coefficients for a denominator in `erf_impl` /// in the interval [8, 11.5]. const ERF_IMPL_HD: &[f64] = &[ 1.0, - 0.588173710611846046373373, - 0.139363331289409746077541, - 0.0166329340417083678763028, - 0.00100023921310234908642639, - 0.24254837521587225125068e-4, + 0.588_173_710_611_846_1, + 0.139_363_331_289_409_75, + 0.016_632_934_041_708_368, + 0.001_000_239_213_102_349_1, + 2.425_483_752_158_722_4e-5, ]; /// Polynomial coefficients for a numerator in `erf_impl` /// in the interval [11.5, 17]. const ERF_IMPL_IN: &[f64] = &[ - -0.00056907993601094962855594, - 0.000169498540373762264416984, - 0.518472354581100890120501e-4, - 0.382819312231928859704678e-5, - 0.824989931281894431781794e-7, + -0.000_569_079_936_010_949_6, + 0.000_169_498_540_373_762_25, + 5.184_723_545_811_009e-5, + 3.828_193_122_319_288_5e-6, + 8.249_899_312_818_944e-8, ]; /// Polynomial coefficients for a denominator in `erf_impl` /// in the interval [11.5, 17]. const ERF_IMPL_ID: &[f64] = &[ 1.0, - 0.339637250051139347430323, - 0.043472647870310663055044, - 0.00248549335224637114641629, - 0.535633305337152900549536e-4, - -0.117490944405459578783846e-12, + 0.339_637_250_051_139_37, + 0.043_472_647_870_310_66, + 0.002_485_493_352_246_371, + 5.356_333_053_371_529e-5, + -1.174_909_444_054_595_8e-13, ]; /// Polynomial coefficients for a numerator in `erf_impl` /// in the interval [17, 24]. const ERF_IMPL_JN: &[f64] = &[ - -0.000241313599483991337479091, - 0.574224975202501512365975e-4, - 0.115998962927383778460557e-4, - 0.581762134402593739370875e-6, - 0.853971555085673614607418e-8, + -0.000_241_313_599_483_991_34, + 5.742_249_752_025_015e-5, + 1.159_989_629_273_837_7e-5, + 5.817_621_344_025_938e-7, + 8.539_715_550_856_736e-9, ]; /// Polynomial coefficients for a denominator in `erf_impl` /// in the interval [17, 24]. const ERF_IMPL_JD: &[f64] = &[ 1.0, - 0.233044138299687841018015, - 0.0204186940546440312625597, - 0.000797185647564398289151125, - 0.117019281670172327758019e-4, + 0.233_044_138_299_687_84, + 0.020_418_694_054_644_03, + 0.000_797_185_647_564_398_3, + 1.170_192_816_701_723_2e-5, ]; /// Polynomial coefficients for a numerator in `erf_impl` /// in the interval [24, 38]. const ERF_IMPL_KN: &[f64] = &[ - -0.000146674699277760365803642, - 0.162666552112280519955647e-4, - 0.269116248509165239294897e-5, - 0.979584479468091935086972e-7, - 0.101994647625723465722285e-8, + -0.000_146_674_699_277_760_36, + 1.626_665_521_122_805_3e-5, + 2.691_162_485_091_652_3e-6, + 9.795_844_794_680_92e-8, + 1.019_946_476_257_234_6e-9, ]; /// Polynomial coefficients for a denominator in `erf_impl` /// in the interval [24, 38]. const ERF_IMPL_KD: &[f64] = &[ 1.0, - 0.165907812944847226546036, - 0.0103361716191505884359634, - 0.000286593026373868366935721, - 0.298401570840900340874568e-5, + 0.165_907_812_944_847_22, + 0.010_336_171_619_150_588, + 0.000_286_593_026_373_868_4, + 2.984_015_708_409_003_4e-6, ]; /// Polynomial coefficients for a numerator in `erf_impl` /// in the interval [38, 60]. const ERF_IMPL_LN: &[f64] = &[ - -0.583905797629771786720406e-4, - 0.412510325105496173512992e-5, - 0.431790922420250949096906e-6, - 0.993365155590013193345569e-8, - 0.653480510020104699270084e-10, + -5.839_057_976_297_718e-5, + 4.125_103_251_054_962e-6, + 4.317_909_224_202_509_4e-7, + 9.933_651_555_900_132e-9, + 6.534_805_100_201_047e-11, ]; /// Polynomial coefficients for a denominator in `erf_impl` /// in the interval [38, 60]. const ERF_IMPL_LD: &[f64] = &[ 1.0, - 0.105077086072039915406159, - 0.00414278428675475620830226, - 0.726338754644523769144108e-4, - 0.477818471047398785369849e-6, + 0.105_077_086_072_039_92, + 0.004_142_784_286_754_756, + 7.263_387_546_445_238e-5, + 4.778_184_710_473_988e-7, ]; /// Polynomial coefficients for a numerator in `erf_impl` /// in the interval [60, 85]. const ERF_IMPL_MN: &[f64] = &[ - -0.196457797609229579459841e-4, - 0.157243887666800692441195e-5, - 0.543902511192700878690335e-7, - 0.317472492369117710852685e-9, + -1.964_577_976_092_295_8e-5, + 1.572_438_876_668_007e-6, + 5.439_025_111_927_009e-8, + 3.174_724_923_691_177e-10, ]; /// Polynomial coefficients for a denominator in `erf_impl` /// in the interval [60, 85]. const ERF_IMPL_MD: &[f64] = &[ 1.0, - 0.052803989240957632204885, - 0.000926876069151753290378112, - 0.541011723226630257077328e-5, - 0.535093845803642394908747e-15, + 0.052_803_989_240_957_63, + 0.000_926_876_069_151_753_3, + 5.410_117_232_266_303e-6, + 5.350_938_458_036_424e-16, ]; /// Polynomial coefficients for a numerator in `erf_impl` /// in the interval [85, 110]. const ERF_IMPL_NN: &[f64] = &[ - -0.789224703978722689089794e-5, - 0.622088451660986955124162e-6, - 0.145728445676882396797184e-7, - 0.603715505542715364529243e-10, + -7.892_247_039_787_227e-6, + 6.220_884_516_609_87e-7, + 1.457_284_456_768_824e-8, + 6.037_155_055_427_153e-11, ]; /// Polynomial coefficients for a denominator in `erf_impl` /// in the interval [85, 110]. const ERF_IMPL_ND: &[f64] = &[ 1.0, - 0.0375328846356293715248719, - 0.000467919535974625308126054, - 0.193847039275845656900547e-5, + 0.037_532_884_635_629_37, + 0.000_467_919_535_974_625_3, + 1.938_470_392_758_456_5e-6, ]; /// `erf_impl` computes the error function at `z`. @@ -334,7 +334,7 @@ fn erf_impl(z: f64, inv: bool) -> f64 { let result = if z < 0.5 { if z < 1e-10 { - z * 1.125 + z * 0.003379167095512573896158903121545171688 + z * 1.125 + z * 0.003_379_167_095_512_573_7 } else { z * 1.125 + z * polynomial(z, ERF_IMPL_AN) / polynomial(z, ERF_IMPL_AD) } diff --git a/src/model/binary.rs b/src/model/binary.rs index 48e8502..97f0041 100644 --- a/src/model/binary.rs +++ b/src/model/binary.rs @@ -6,6 +6,12 @@ use crate::kernel::Kernel; use crate::observation::*; use crate::storage::Storage; +/* +trait BinaryObeservation: Observation { + fn new(storage: &mut Storage, elems: &[(usize, f64)], t: f64, margin: f64) -> Self; +} +*/ + pub enum BinaryModelObservation { Probit, Logit, diff --git a/src/utils.rs b/src/utils.rs index 8868fb3..08c3964 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -20,20 +20,20 @@ const CS: [f64; 14] = [ ]; const RS: [f64; 5] = [ - 1.2753666447299659525, - 5.019049726784267463450, - 6.1602098531096305441, - 7.409740605964741794425, - 2.9788656263939928886, + 1.275_366_644_729_965_9, + 5.019_049_726_784_267, + 6.160_209_853_109_631, + 7.409_740_605_964_742, + 2.978_865_626_393_993, ]; const QS: [f64; 6] = [ - 2.260528520767326969592, - 9.3960340162350541504, - 12.048951927855129036034, - 17.081440747466004316, - 9.608965327192787870698, - 3.3690752069827527677, + 2.260_528_520_767_327, + 9.396_034_016_235_054, + 12.048_951_927_855_128, + 17.081_440_747_466_004, + 9.608_965_327_192_788, + 3.369_075_206_982_752_8, ]; /// Normal cumulative density function. @@ -59,7 +59,7 @@ pub fn logphi(z: f64) -> (f64, f64) { (res, dres) } else if z < -11.3137 { // Second case: z very small. - let mut num = 0.5641895835477550741; + let mut num = 0.564_189_583_547_755_1; for r in &RS { num = -z * num / SQRT_2 + r;