octave-bug-tracker
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Octave-bug-tracker] [bug #60682] betainc is inaccurate


From: Nicholas Jankowski
Subject: [Octave-bug-tracker] [bug #60682] betainc is inaccurate
Date: Wed, 11 Aug 2021 19:47:07 -0400 (EDT)
User-agent: Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/92.0.4515.131 Safari/537.36

Follow-up Comment #4, bug #60682 (project octave):

curious how matlab compares, here are the outputs of the two commands below in
Matlab 2021a:


>> x=2.^(-(0:20))';[x x.^2 betainc(x,2,1)./x.^2-1]

ans =

     1.000000000000000e+00     1.000000000000000e+00                        
0
     5.000000000000000e-01     2.500000000000000e-01    
2.220446049250313e-16
     2.500000000000000e-01     6.250000000000000e-02    
2.220446049250313e-16
     1.250000000000000e-01     1.562500000000000e-02    
2.220446049250313e-16
     6.250000000000000e-02     3.906250000000000e-03    
2.220446049250313e-16
     3.125000000000000e-02     9.765625000000000e-04   
-4.440892098500626e-16
     1.562500000000000e-02     2.441406250000000e-04    
8.881784197001252e-16
     7.812500000000000e-03     6.103515625000000e-05   
-1.110223024625157e-16
     3.906250000000000e-03     1.525878906250000e-05    
1.110223024625157e-15
     1.953125000000000e-03     3.814697265625000e-06   
-7.771561172376096e-16
     9.765625000000000e-04     9.536743164062500e-07   
-4.440892098500626e-16
     4.882812500000000e-04     2.384185791015625e-07    
4.440892098500626e-16
     2.441406250000000e-04     5.960464477539063e-08                        
0
     1.220703125000000e-04     1.490116119384766e-08                        
0
     6.103515625000000e-05     3.725290298461914e-09    
2.220446049250313e-15
     3.051757812500000e-05     9.313225746154785e-10   
-2.220446049250313e-16
     1.525878906250000e-05     2.328306436538696e-10    
1.776356839400250e-15
     7.629394531250000e-06     5.820766091346741e-11    
6.661338147750939e-16
     3.814697265625000e-06     1.455191522836685e-11    
4.440892098500626e-16
     1.907348632812500e-06     3.637978807091713e-12    
2.220446049250313e-16
     9.536743164062500e-07     9.094947017729282e-13                        
0


>> x=2.^-18*(1+eps(1)*(0:9)');[x x.^2 betainc(x,2,1) betainc(x,2,1)./x.^2-1]

ans =

     3.814697265625000e-06     1.455191522836685e-11     1.455191522836686e-11
    4.440892098500626e-16
     3.814697265625001e-06     1.455191522836686e-11     1.455191522836686e-11
                        0
     3.814697265625002e-06     1.455191522836686e-11     1.455191522836686e-11
   -4.440892098500626e-16
     3.814697265625003e-06     1.455191522836687e-11     1.455191522836686e-11
   -8.881784197001252e-16
     3.814697265625003e-06     1.455191522836688e-11     1.455191522836686e-11
   -1.332267629550188e-15
     3.814697265625004e-06     1.455191522836688e-11     1.455191522836691e-11
    1.776356839400250e-15
     3.814697265625005e-06     1.455191522836689e-11     1.455191522836691e-11
    1.332267629550188e-15
     3.814697265625006e-06     1.455191522836690e-11     1.455191522836691e-11
    8.881784197001252e-16
     3.814697265625007e-06     1.455191522836690e-11     1.455191522836691e-11
    4.440892098500626e-16
     3.814697265625008e-06     1.455191522836691e-11     1.455191522836691e-11
                        0



using Octave 6.3.0:


>> x=2.^(-(0:20))';[x x.^2 betainc(x,2,1)./x.^2-1]
ans =

   1.000000000000000e+00   1.000000000000000e+00                       0
   5.000000000000000e-01   2.500000000000000e-01   2.220446049250313e-16
   2.500000000000000e-01   6.250000000000000e-02                       0
   1.250000000000000e-01   1.562500000000000e-02   4.440892098500626e-16
   6.250000000000000e-02   3.906250000000000e-03   2.220446049250313e-16
   3.125000000000000e-02   9.765625000000000e-04                       0
   1.562500000000000e-02   2.441406250000000e-04   6.661338147750939e-16
   7.812500000000000e-03   6.103515625000000e-05   4.440892098500626e-16
   3.906250000000000e-03   1.525878906250000e-05   4.440892098500626e-16
   1.953125000000000e-03   3.814697265625000e-06   2.220446049250313e-16
   9.765625000000000e-04   9.536743164062500e-07                       0
   4.882812500000000e-04   2.384185791015625e-07  -1.110223024625157e-16
   2.441406250000000e-04   5.960464477539062e-08   1.554312234475219e-15
   1.220703125000000e-04   1.490116119384766e-08  -2.331468351712829e-15
   6.103515625000000e-05   3.725290298461914e-09   1.110223024625157e-15
   3.051757812500000e-05   9.313225746154785e-10  -2.664535259100376e-15
   1.525878906250000e-05   2.328306436538696e-10   6.661338147750939e-16
   7.629394531250000e-06   5.820766091346741e-11   6.661338147750939e-16
   3.814697265625000e-06   1.455191522836685e-11  -3.108624468950438e-15
   1.907348632812500e-06   3.637978807091713e-12   2.220446049250313e-16
   9.536743164062500e-07   9.094947017729282e-13                       0

>> x=2.^-18*(1+eps(1)*(0:9)');[x x.^2 betainc(x,2,1) betainc(x,2,1)./x.^2-1]
ans =

   3.814697265625000e-06   1.455191522836685e-11   1.455191522836681e-11 
-3.108624468950438e-15
   3.814697265625001e-06   1.455191522836686e-11   1.455191522836681e-11 
-3.552713678800501e-15
   3.814697265625002e-06   1.455191522836686e-11   1.455191522836681e-11 
-3.996802888650564e-15
   3.814697265625003e-06   1.455191522836687e-11   1.455191522836681e-11 
-4.440892098500626e-15
   3.814697265625003e-06   1.455191522836688e-11   1.455191522836681e-11 
-4.884981308350689e-15
   3.814697265625004e-06   1.455191522836688e-11   1.455191522836691e-11  
1.776356839400250e-15
   3.814697265625005e-06   1.455191522836689e-11   1.455191522836691e-11  
1.332267629550188e-15
   3.814697265625006e-06   1.455191522836690e-11   1.455191522836691e-11  
8.881784197001252e-16
   3.814697265625007e-06   1.455191522836690e-11   1.455191522836691e-11  
4.440892098500626e-16
   3.814697265625008e-06   1.455191522836691e-11   1.455191522836691e-11      
                0


I thought maybe the patch was already pushed but 

Octave 6.3.0:

>> betainc (0.5, 1, Inf)
ans = NaN



after applying the patch i get:


>> x=2.^(-(0:20))';[x x.^2 betainc(x,2,1)./x.^2-1]
ans =

   1.000000000000000e+00   1.000000000000000e+00                       0
   5.000000000000000e-01   2.500000000000000e-01                       0
   2.500000000000000e-01   6.250000000000000e-02                       0
   1.250000000000000e-01   1.562500000000000e-02                       0
   6.250000000000000e-02   3.906250000000000e-03                       0
   3.125000000000000e-02   9.765625000000000e-04                       0
   1.562500000000000e-02   2.441406250000000e-04                       0
   7.812500000000000e-03   6.103515625000000e-05                       0
   3.906250000000000e-03   1.525878906250000e-05                       0
   1.953125000000000e-03   3.814697265625000e-06                       0
   9.765625000000000e-04   9.536743164062500e-07                       0
   4.882812500000000e-04   2.384185791015625e-07                       0
   2.441406250000000e-04   5.960464477539062e-08                       0
   1.220703125000000e-04   1.490116119384766e-08                       0
   6.103515625000000e-05   3.725290298461914e-09                       0
   3.051757812500000e-05   9.313225746154785e-10                       0
   1.525878906250000e-05   2.328306436538696e-10                       0
   7.629394531250000e-06   5.820766091346741e-11                       0
   3.814697265625000e-06   1.455191522836685e-11                       0
   1.907348632812500e-06   3.637978807091713e-12                       0
   9.536743164062500e-07   9.094947017729282e-13                       0


>> x=2.^-18*(1+eps(1)*(0:9)');[x x.^2 betainc(x,2,1) betainc(x,2,1)./x.^2-1]
ans =

   3.814697265625000e-06   1.455191522836685e-11   1.455191522836685e-11      
                0
   3.814697265625001e-06   1.455191522836686e-11   1.455191522836686e-11      
                0
   3.814697265625002e-06   1.455191522836686e-11   1.455191522836686e-11      
                0
   3.814697265625003e-06   1.455191522836687e-11   1.455191522836687e-11      
                0
   3.814697265625003e-06   1.455191522836688e-11   1.455191522836688e-11      
                0
   3.814697265625004e-06   1.455191522836688e-11   1.455191522836688e-11      
                0
   3.814697265625005e-06   1.455191522836689e-11   1.455191522836689e-11      
                0
   3.814697265625006e-06   1.455191522836690e-11   1.455191522836690e-11      
                0
   3.814697265625007e-06   1.455191522836690e-11   1.455191522836690e-11      
                0
   3.814697265625008e-06   1.455191522836691e-11   1.455191522836691e-11      
                0


>> betainc (0.5, 1, Inf)
ans = 1.000000000000000e+00


the trivial case fix is a good fix.  but I'm not seeing the accuracy issue 

    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?60682>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/




reply via email to

[Prev in Thread] Current Thread [Next in Thread]