Integrating a function with apache.math3, converting code from python

Viewed 22

I'm having trouble integrating a function using the apache.math3 library - I converted the code from a Python snippet I found online. I believe the first two methods (pdflp,cdfmaxlp) are working but not sure I found the correct method to integrate over a function.

Perhaps SimpsonIntegrator isn't correct integrator, but I couldn't find documentation on how scipy.integrate.quad was implemented

Any insight would be appreciated, Thanks!

Original Python

import scipy
import math 
from scipy import integrate 
from scipy import stats

p0 = float(128.60)
v = float(.2149)
p1 = float(130)
p2 = float(130)
t1 = float(1)
t2 = float(504)

def pdflp(x,p0,v,t1): 
  return scipy.stats.norm(math.log(p0), math.sqrt(t1)*v/math.sqrt(365.2425*24)).pdf(x)

def cdfmaxlp(x,v,t1,t2):
  return (1-scipy.special.erf(x/(v*math.sqrt(t2-t1)/math.sqrt(24*365.2425/math.sqrt(2.0)))))

def upandin(p0,v,p1,p2,t1,t2):
  integrand = lambda x: (pdflp(x,p0,v,t1) * cdfmaxlp(math.log(p1)-x,v,t1,t2))
  return scipy.integrate.quad(integrand,0,math.log(p1))

//upandin(p0,v,p1,p2,t1,t2)
//(0.8027563748048668, 2.5901298362509227e-09)

Converted Scala

import org.apache.commons.math3.analysis.integration.SimpsonIntegrator
import org.apache.commons.math3.distribution.NormalDistribution
import org.apache.commons.math3.special.Erf
import org.apache.commons.math3.analysis.UnivariateFunction

val p0: Double = 128.60
val v: Double  = .2149
val p1: Double = 130.0
val p2: Double = 130.0
val t1: Double = 1.0
val t2: Double = 504.0

def pdflp(x: Double, p0: Double, v: Double, t1: Double) = {
  val norm = new NormalDistribution(math.log(p0), math.sqrt(t1) * v / math.sqrt(365.2425 * 24))
  norm.density(x)
}

def cdfmaxlp(x: Double, v: Double, t1: Double, t2: Double) = {
  val erf = Erf.erf(x / (v * math.sqrt(t2 - t1) / math.sqrt(24 * 365.2425 / math.sqrt(2.0))))
  (1 - erf)
}

//believe the issue is in here, perhaps SimpsonIntegrator isn't correct 
//but I couldn't find documentation on how `scipy.integrate.quad` was implemented 
def upandin(p0: Double, v: Double, p1: Double, p2: Double, t1: Double, t2: Double) = {
  val function: UnivariateFunction = (x: Double) => pdflp(x, p0, v, t1) * cdfmaxlp(math.log(p1) - x, v, t1, t2)
  val integrator                   = new SimpsonIntegrator(1e-4, 1e-8, 1, 64)
  val i                            = integrator.integrate(100, function, 0, math.log(p1))

  i
}

//upandin(p0, v, p1, p2, t1, t2)
//-0.0020752245549497347
0 Answers
Related