B4A Library Probability Library

Hi
After viewing Erel's video I managed to produce a simple library.
I used the C# code of the B4PPC Probability library and very easily converted it to Java.

The library caculates the Probability density and commulative functions of:
Normal, ChiSquare, Student, Fisher F, Rayleigh, Erlang, Binomial, Geometric, Poisson and Exponential.
It also includes Random Number Generators with some of the distributions (as B4a has only uniform distribution random numbers), and some mathematical functions that are related to distributions: Beta, Gamma, Factorial and Double Factorial, NChooseK

Please check and let me know if it is understandable.
Edit: updated to correct documentation.
Edit: This version of documentation must be better...
Edit: attached an example program.
 

Attachments

Last edited:
Click to expand...
C# -> Java ?

Hi,
Please, can you tell me how you convert easily C# in JAVA ?
Thank's
 
The basic languages are similar, I'm not talking about the .net.
Here is an example of the same method in both:
In C# :
Code: 
public  double [] StdNormal(double Z, string Mode)
      {
      if ((Mode != "L") & (Mode != "C") & (Mode != "R") & (Mode != "T"))
         return new double[] { -999, -999 } ;   
      double z2 = Z*Z ;
      double z3 = z2*Z ;
      int SignZ = 1 ;
      if (Z < 0) 
      {
         Z = Math.Abs(Z) ;         
         SignZ = -1 ;
      }
      double Fn = Math.Pow(Math.E,(-z2/2))/Math.Sqrt(2*Math.PI) ;
      double Pn = 1-0.5*Math.Pow((1 + 0.0498673470 * Z 
                                    + 0.0211410061 * z2 
                                    + 0.0032776263 * z3
                                    + 0.0000380036 * z2 * z2
                                    + 0.0000488906 * z2 * z3
                                    + 0.0000053830 * z3 * z3),-16) ;
      switch(Mode)
         {
         case "T": 
               Pn = (1-Pn)*2;
               break;
         case "R":
               Pn = (1-Pn) ;
               break ;
         case "L":
               if (SignZ == -1) 
                  Pn = 1-Pn ;
               break ;
         case "C":
               Pn = Pn - (1-Pn) ;
               break;
         default:
               if (SignZ == -1) 
                     Pn = 1-Pn ;
               break ;   
         }
      if (Pn > 1)
         Pn = 1 ;
      return new double[] { Pn, Fn } ;
      }

and in Java :
Code: 
public  double [] StdNormal(double Z, String Mode)
   {
      /*
       * Returns an array of two Double variables, the first is the Standard Normal Probability P(z) and the second is the Density F(z).
       */
   if ((Mode != "L") & (Mode != "C") & (Mode != "R") & (Mode != "T"))
      return new double[] { -999, -999 } ;   
   double z2 = Z*Z ;
   double z3 = z2*Z ;
   int SignZ = 1 ;
   if (Z < 0) 
   {
      Z = Math.abs(Z) ;         
      SignZ = -1 ;
   }
   double Fn = Math.pow(Math.E,(-z2/2))/Math.pow(2*Math.PI, 0.5) ;
   double Pn = 1-0.5*Math.pow((1 + 0.0498673470 * Z 
                                 + 0.0211410061 * z2 
                                 + 0.0032776263 * z3
                                 + 0.0000380036 * z2 * z2
                                 + 0.0000488906 * z2 * z3
                                 + 0.0000053830 * z3 * z3),-16) ;
   
   if (Mode == "C") 
      Pn = 2*Pn -1 ;
   else if (Mode == "R")
      Pn = (1-Pn) ;
   else if (Mode == "L")
      {if (SignZ == -1) 
         Pn = 1-Pn ;}
   else if (Mode == "T")
      Pn = (1-Pn)*2;
   else
      if (SignZ == -1) 
         Pn = 1-Pn ;
   
   if (Pn > 1)
      Pn = 1 ;
   return new double[] { Pn, Fn } ;
   }

The differences are string - String, Math.Pow - Math.pow , Math.sqrt - Math.pow , Math.Abs - Math.abs, the switch doesn't work for strings so I changed it to If. The Eclipse IDE shows the errors and suggests the correction. Isn't it easy ?
 
Last edited:
David,
thanks for sharing your lib.
Many people are struggling with different library wrappers and if you could post your Eclipse library wrapper project so we could pick it apart and learn from it, that would be very helpful. Here is the link to the Library Wrapper Tutorial request, many people voted and the vote count keeps growing:
http://www.b4x.com/forum/bugs-wishl...ndroid-java-class-b4a-library.html#post151416

Thanks!:sign0188:
 
No problem, only this is not a wrapper, just a piece of code built by Erel's examples.
Click to expand...
Thanks David, that's still useful - one can play with your library code or create a different one using yours as a template and try to compile a library. I sure will.:sign0188::sign0098:
 
I actually was not aware of such distribution (and now read about in wikipedia). The application does not include it but like you wrote above - it should be easy to add it.
You can calculate ERF by this series:
07e87fcd0fa7ce0f9622d1fe8c332b37.png

and get the distribution:
PDF
9e2b928f871663bc2ab9e7478735f4e2.png

CDF
2049e219bac8f47a956631cd9eed2b7d.png

Mean
091d1c00bddc99356e54a72830581b8f.png
 
Last edited:
I have coded these functions but not added them to the library:

Code: 
' Error function:
Sub erf(z As Double) As Double
Dim  res, fct, sum As Double
Dim k As Int
sum = 0
For i = 0 To 30
If i = 0 Then fct = 1    Else fct = fct * i
    k = 2*i+1
    sum = sum + Power(-1,i) * Power(z,k)/(fct*k)
Next
res = 1.1283791670955126 *sum  ' = (2/sqrt(cpi))  * sum
Return NumberFormat(res,0,7)
End Sub

'LogNormal density:
Sub DpLogNorm(x As Double, mu As Double,sigma As Double) As Double
Dim res As Double
res = Power(cE,-Power((Logarithm(x,cE)-mu),2)/(2*sigma*sigma)) / (x * sigma * (Sqrt((2*cPI))))
Return NumberFormat(res,0,7)
End Sub

'LogNormal cumulative:
Sub Cplognorm(x As Double, mu As Double,sigma As Double) As Double
Dim res As Double
res = 0.5 + 0.5*erf((Logarithm(x,cE)-mu)/(Sqrt(2)*sigma))
Return NumberFormat(res,0,7)
End Sub

The last function can be implemented by the following as well but I have found it to be less accurate :
Code: 
Sub Cplognorm(x As Double, mu As Double,sigma As Double) As Double()
Return prob.StdNormal((Logarithm(x,cE)-mu)/sigma, "L")
End Sub
 
Last edited:
To add to the collection here is another ERFC function and an Inverse Normal cumulative distribution:

Code: 
#Region Normal Distribution Formula
Sub ERFC(x As Double) As Double
  ' Uses Chebyshev curve fit to approximate the complimentary error function
  ' with fractional error every where less than 1.2 x 10-7.
  ' Note that ERFC(x) = 1-ERF(x)and x is typically between -4 and +4
 
  Dim t, z, ans As Double
 
  z = Abs(x)
  t = 1 / (1 + 0.5 * z)
  ans = t * Power(cE,(-z * z - 1.26551223 + t * (1.00002368 + t * (0.37409196 + t * (0.09678418 + t * (-0.18628806 + t * (0.27886807 + t * (-1.13520398 + t * (1.48851587 + t * (-0.82215223 + t * 0.17087277))))))))))
  If x >= 0 Then
    Return ans
  Else
    Return 2 - ans
  End If
End Sub

Sub InvNormCDF(p As Double) As Double
'  Calculate inverse Normal cummulative distribution
'  using rational approximation
'  "An algorithm for computing the inverse normal cumulative distribution function"
'  (http://home.online.no/~pjacklam/notes/invnorm/)
'  by John Herrero (3-Jan-03)

'Define coefficients in rational approximations

  Dim a1  As Double = -39.6968302866538
  Dim a2  As Double = 220.946098424521
  Dim a3  As Double = -275.928510446969
  Dim a4  As Double = 138.357751867269
  Dim a5  As Double = -30.6647980661472
  Dim a6  As Double = 2.50662827745924

  Dim b1  As Double = -54.4760987982241
  Dim b2  As Double = 161.585836858041
  Dim b3  As Double = -155.698979859887
  Dim b4  As Double = 66.8013118877197
  Dim b5  As Double = -13.2806815528857

  Dim c1  As Double = -7.78489400243029E-03
  Dim c2  As Double = -0.322396458041136
  Dim c3  As Double = -2.40075827716184
  Dim c4  As Double = -2.54973253934373
  Dim c5  As Double = 4.37466414146497
  Dim c6  As Double = 2.93816398269878

  Dim d1  As Double = 7.78469570904146E-03
  Dim d2  As Double = 0.32246712907004
  Dim d3  As Double = 2.445134137143
  Dim d4  As Double = 3.75440866190742

'Define break-points
  Dim p_low As Double  = 0.02425
  Dim p_high As Double = 1 - p_low

'Define work variables
  Dim q As Double, r As Double

'If argument out of bounds, raise error
  If p <= 0 Or p >= 1 Then Return 999

  If p < p_low Then
  'Rational approximation for lower region
    q = Sqrt(-2 * Logarithm(p,cE))
    Return (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / _
           ((((d1 * q + d2) * q + d3) * q + d4) * q + 1)
  Else
' -----------------------------------------------
    If p <= p_high Then
    'Rational approximation for central region
      q = p - 0.5
      r = q * q
      Return (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q / _
             (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1)
    Else
'   -----------------------------------------------  
     If p < 1 Then
     'Rational approximation for upper region
       q = Sqrt(-2 * Logarithm(1 - p,cE))
       Return -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / _
               ((((d1 * q + d2) * q + d3) * q + d4) * q + 1)
     End If
   End If
' ----------------------------------------------
  End If 
  Return 999
End Sub
#End Region
 
To add to the collection here is another ERFC function and an Inverse Normal cumulative distribution:
Click to expand...
Thank you bluejay. The Library already has a InvNormal but it does it by another algorithm. Would be interesting to compare.
The code below is similar to the library's code.

Code: 
Sub InvNormal(Pn As Double, Mean As Double, Variance As Double, NMod As String) As Double
Dim a As Double
a = InvStdNormal(Pn, NMod)
a =  a * Sqrt(Variance) + Mean
Return a
End Sub

Sub InvStdNormal(Pn As Double,  NMod As String) As Double
Dim signz, halfp As Int

signz = 1
halfp = 1
If Pn > 0.5 Then
    halfp = -1
    Pn = 1-Pn
End If
Select NMod
    Case "T"
        If (halfp = -1) Then
            Pn = 0.5-Pn/2
        Else
            Pn = Pn/2
        End If
    Case "R"
            If (halfp = -1) Then signz = -1

    Case "L"
            If (halfp = 1) Then signz = -1

    Case "C"
            If (halfp = -1) Then
                Pn = Pn/2
            Else
                Pn = 0.5-Pn/2
            End If
End Select
Return  signz * Abr(Pn)
End Sub

Sub Abr(Pz As Double) As Double
Dim  t,tr,nom,denom As Double
t = Logarithm(1/(Pz*Pz),cE)
tr = Sqrt(t)
nom = 2.515517 + 0.802853 * tr + 0.010328 * t
denom = 1 + 1.432788 * tr + 0.189269 * t + 0.001308 *t*tr
Return  tr - nom / denom
End Sub
 
Last edited:
This is the inverse Normal cumulative distribution not inverse normal distribution.
 
I don't understand the difference.
I compared the results with this code:
Code: 
Dim q As Double
For i = 1 To 9
    q = i/10
    Log(NumberFormat(i/10,1,1) & " " & InvNormCDF(q) & " " & prob.InvStdNormal(q,"L"))
Next
got these results:
0.1 -1.281551564140072 -1.281728756502709
0.2 -0.8416212327266098 -0.8414567173547837
0.3 -0.5244005132792943 -0.5240018703826799
0.4 -0.2533471028599986 -0.25293326782658254
0.5 0 1.0100667546808495E-7
0.6 0.2533471028599986 0.25293326782658254
0.7 0.5244005132792942 0.5240018703826799
0.8 0.8416212327266092 0.841456717354784
0.9 1.281551564140072 1.281728756502709
Click to expand...

It seems that this is the same thing although there are differences in accuracy.
 
You must log in or register to reply here.

Similar Threads

J
  • Question Question
Android Question error QR scanner with R8
Replies
1
Views
26
Erel
Erel
Erel
Server software is being upgraded
Replies
3
Views
97
LucaMs
LucaMs
A
eMedNY
Replies
0
Views
95
Alex_197
A
Daestrum
When tech works for you. (minimum effort)
Replies
0
Views
123
Daestrum
Daestrum
Robert Valentino
  • Question Question
Android Question #R8 and B4XLib crashing method not found in: okhttp3.Request$Builder
Replies
5
Views
174
Erel
Erel
Back
Top