AssertionError: [MInfSegment, [-inf, 0.05648840000000001]]==MInfSegment, [-inf, 0.05648840000000001]==Segment, [0.05648840000000001, -0.08108632]==-inf==0.0564884000000000

[michelleowen]

AssertionError: [MInfSegment, [-inf, 0.05648840000000001]]==MInfSegment, [-inf, 0.05648840000000001]==Segment, [0.05648840000000001, -0.08108632]==-inf==0.0564884000000000

[jszymon]

Hi Michelle, thanks for reporting the problem. Two questions: 1. do you use the latest git version? 2. would it be possible for you to send me the actual code? Otherwise it is hard for me to trace the problem. I suspect that during optimization the resulting distribution becomes degenerate in some sense. Nevertheless, this should not produce an error, so I'd like to fix this. All the best Szymon

…

On Sat, Jun 16, 2018 at 11:36 PM, michelleowen ***@***.***> wrote: I try to use pacal to create pdf of summation of three random variables. Then I use Nelder-Mead to minimize the negative loglikelihood. However, after several iterations, I came across this error message raised by segments.py in pacal package. param: [2. 0. 0.2] processed in 7.367134094238281e-05 seconds --- 5289.1225545465795 param: [2.1 0. 0.2] processed in 9.1552734375e-05 seconds --- 5413.69569454499 param: [2.0e+00 2.5e-04 2.0e-01] processed in 0.0001373291015625 seconds --- 5289.2735727455265 param: [2. 0. 0.21] processed in 7.677078247070312e-05 seconds --- 5291.688762705109 param: [1.90000000e+00 1.66666667e-04 2.06666667e-01] processed in 0.00013709068298339844 seconds --- 5162.96097874654 param: [1.8e+00 2.5e-04 2.1e-01] processed in 0.0030710697174072266 seconds --- 5032.661772965377 param: [1.86666667e+00 3.33333333e-04 1.96666667e-01] processed in 9.465217590332031e-05 seconds --- 5116.824460425303 param: [1.77777778e+00 1.38888889e-04 2.04444444e-01] processed in 0.00011873245239257812 seconds --- 5001.204618755317 param: [1.66666667e+00 8.33333333e-05 2.06666667e-01] processed in 0.00010085105895996094 seconds --- 4851.356632260949 param: [1.55555556e+00 4.44444444e-04 2.08888889e-01] processed in 8.20159912109375e-05 seconds --- 4698.843942739831 param: [1.33333333e+00 6.66666667e-04 2.13333333e-01] processed in 0.00016498565673828125 seconds --- 4389.884695448772 param: [1.33333333e+00 3.33333333e-04 2.23333333e-01] processed in 9.131431579589844e-05 seconds --- 4394.404234423141 param: [1.08888889e+00 4.72222222e-04 2.18888889e-01] processed in 7.963180541992188e-05 seconds --- 4069.8662612067 param: [7.33333333e-01 5.83333333e-04 2.23333333e-01] processed in 0.00010180473327636719 seconds --- 3861.2544906000508 param: [0.6 0.00097222 0.23333333] processed in 8.392333984375e-05 seconds --- 4055.4810152942673 param: [0.44444444 0.00114815 0.22333333] processed in 0.0001323223114013672 seconds --- 4867.380442059709 param: [1.11111111e+00 5.37037037e-04 2.23333333e-01] processed in 8.893013000488281e-05 seconds --- 4099.60266263902 param: [0.2962963 0.0007284 0.24 ] processed in 9.965896606445312e-05 seconds --- 6952.753469381866 param: [1.07407407e+00 6.82098765e-04 2.20000000e-01] processed in 9.369850158691406e-05 seconds --- 4053.412098162286 param: [0.49382716 0.00095473 0.22777778] processed in 9.465217590332031e-05 seconds --- 4496.614564486516 param: [9.56790123e-01 6.41460905e-04 2.24444444e-01] processed in 0.0001659393310546875 seconds --- 3933.7490338058774 param: [1.24279835e+00 2.99039781e-04 2.11851852e-01] processed in 0.00011301040649414062 seconds --- 4264.013825788384 param: [0.76069959 0.00080393 0.22796296] processed in 0.00011706352233886719 seconds --- 3854.024649773798 param: [0.55980796 0.00067038 0.23049383] processed in 8.177757263183594e-05 seconds --- 4179.071532320201 param: [9.45507545e-01 6.79169524e-04 2.22623457e-01] processed in 7.414817810058594e-05 seconds --- 3922.888324502457 param: [0.66957019 0.00073616 0.22483539] processed in 9.107589721679688e-05 seconds --- 3920.1594117114482 param: [0.49689453 0.00073644 0.228131 ] processed in 8.392333984375e-05 seconds --- 4476.307250583266 param: [8.33354291e-01 6.93487877e-04 2.24000343e-01] processed in 8.225440979003906e-05 seconds --- 3855.589457865707 param: [8.82021287e-01 6.51006473e-04 2.25362369e-01] processed in 0.00012493133544921875 seconds --- 3878.5581421462252 param: [8.28908512e-01 6.72294540e-04 2.25230624e-01] processed in 8.392333984375e-05 seconds --- 3854.914557579754 param: [8.81974928e-01 8.63139352e-04 2.28129287e-01] processed in 0.00015616416931152344 seconds --- 3880.874955501708 param: [7.70493732e-01 6.53284838e-04 2.24532322e-01] processed in 0.0001010894775390625 seconds --- 3849.8283418083765 param: [7.40046931e-01 7.26182784e-04 2.27816930e-01] processed in 8.344650268554688e-05 seconds --- 3860.5974550157202 param: [8.10027451e-01 7.01661603e-04 2.24954490e-01] processed in 7.557868957519531e-05 seconds --- 3850.5827846618563 param: [0.73190533 0.00076695 0.22640256] processed in 0.00021266937255859375 seconds --- 3864.0475105394494 param: [8.04657718e-01 6.95959446e-04 2.25523608e-01] processed in 0.0002453327178955078 seconds --- 3850.1306397142275 param: [8.29419679e-01 5.63343980e-04 2.22043983e-01] processed in 8.988380432128906e-05 seconds --- 3852.830156121806 param: [8.12239656e-01 6.23489638e-04 2.23523728e-01] processed in 7.796287536621094e-05 seconds --- 3849.8769339170635 param: [7.81566620e-01 6.13494344e-04 2.24098615e-01] processed in 9.965896606445312e-05 seconds --- 3848.289064454587 param: [7.67336205e-01 5.69410715e-04 2.23670678e-01] processed in 7.510185241699219e-05 seconds --- 3849.618741910505 param: [7.71542287e-01 5.64219768e-04 2.22579502e-01] processed in 9.179115295410156e-05 seconds --- 3848.3443620554585 param: [7.36828770e-01 5.97176329e-04 2.23949898e-01] processed in 0.0001723766326904297 seconds --- 3859.835450980743 param: [7.93386935e-01 6.16911311e-04 2.23630271e-01] processed in 0.00013828277587890625 seconds --- 3847.90997596603 param: [7.93836829e-01 5.43132110e-04 2.22339937e-01] processed in 0.00010633468627929688 seconds --- 3846.965369220417 param: [8.05508378e-01 4.88055746e-04 2.21243745e-01] processed in 7.462501525878906e-05 seconds --- 3847.1230922372747 param: [8.07651302e-01 6.18138742e-04 2.24133046e-01] processed in 8.988380432128906e-05 seconds --- 3849.4799136312854 param: [7.80569541e-01 5.77699511e-04 2.22967888e-01] processed in 0.0001895427703857422 seconds --- 3847.624068337901 param: [7.96962250e-01 5.45000944e-04 2.21860115e-01] processed in 0.00012063980102539062 seconds --- 3846.8595405953224 param: [8.04660065e-01 5.10754244e-04 2.20740865e-01] processed in 7.963180541992188e-05 seconds --- 3846.8092937871143 param: [7.92657356e-01 4.70812600e-04 2.20402190e-01] processed in 0.00019478797912597656 seconds --- 3845.616883105677 param: [7.92292566e-01 3.97763244e-04 2.18788150e-01] processed in 0.00010895729064941406 seconds --- 3844.477681428711 param: [8.13290099e-01 3.90066887e-04 2.18278080e-01] processed in 9.059906005859375e-05 seconds --- 3846.2981525343926 param: [8.12991657e-01 3.22590806e-04 2.16198126e-01] processed in 0.00011730194091796875 seconds --- 3844.841948197754 param: [8.07722817e-01 2.29526382e-04 2.14768705e-01] processed in 7.653236389160156e-05 seconds --- 3842.915472472104 param: [8.09254193e-01 8.89124506e-05 2.11782625e-01] processed in 7.677078247070312e-05 seconds --- 3840.999820460356 param: [7.96402178e-01 1.49444113e-04 2.12901187e-01] processed in 7.796287536621094e-05 seconds --- 3840.5714275841674 param: [7.87958218e-01 2.91327265e-05 2.10212741e-01] processed in 8.0108642578125e-05 seconds --- 3838.60465661308 param: [7.80011660e-01 2.12814745e-05 2.10990884e-01] processed in 7.748603820800781e-05 seconds --- 3839.422861748212 param: [ 7.92523481e-01 -3.04878810e-04 2.03202684e-01] processed in 8.559226989746094e-05 seconds --- 3833.724530052107 param: [ 7.92638939e-01 -6.56199837e-04 1.95409951e-01] processed in 9.489059448242188e-05 seconds --- 3828.458931342521 param: [ 7.64485019e-01 -4.92769541e-04 1.99293092e-01] processed in 0.00016760826110839844 seconds --- 3834.26295517022 param: [ 7.83376457e-01 -7.67839242e-04 1.92286305e-01] processed in 0.00013875961303710938 seconds --- 3826.8256800975155 param: [ 0.78505885 -0.0011624 0.18293402] processed in 7.748603820800781e-05 seconds --- 3820.6920915476903 param: [ 0.77349699 -0.00157005 0.17487863] processed in 9.608268737792969e-05 seconds --- 3816.8808179592206 param: [ 0.76626638 -0.00236963 0.15721158] processed in 0.0001201629638671875 seconds --- 3807.7752830896975 param: [ 0.79815776 -0.00229939 0.15774394] processed in 7.915496826171875e-05 seconds --- 3804.5435945323006 param: [ 0.81499413 -0.00320269 0.13696936] processed in 0.00013756752014160156 seconds --- 3793.957732312433 param: [ 0.7849073 -0.00383362 0.12266668] processed in 0.00015473365783691406 seconds --- 3785.3001119282135 param: [ 0.78104149 -0.00542233 0.08629505] processed in 0.0001010894775390625 seconds --- 3768.599582643904 param: [ 0.78980914 -0.00616737 0.07071664] processed in 7.772445678710938e-05 seconds --- 3761.2207949550857 param: [ 0.79218429 -0.00866986 0.01460795] processed in 0.0001475811004638672 seconds --- 3744.789138756524 param: [ 0.82588023 -0.00916029 0.00137 ] processed in 7.534027099609375e-05 seconds --- 3746.9002257836974 param: [ 0.78440987 -0.01229896 -0.06878736] ------------------------------ AssertionError Traceback (most recent call last) in () ----> 1 res2 = minimize(loglikelihood, [2.0,0,0.2], method='Nelder-Mead', tol=1e-6) /usr/local/lib/python3.4/dist-packages/scipy/optimize/_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options) 473 callback=callback, **options) 474 elif meth == 'nelder-mead': --> 475 return _minimize_neldermead(fun, x0, args, callback, **options) 476 elif meth == 'powell': 477 return _minimize_powell(fun, x0, args, callback, **options) /usr/local/lib/python3.4/dist-packages/scipy/optimize/optimize.py in _minimize_neldermead(func, x0, args, callback, maxiter, maxfev, disp, return_all, initial_simplex, xatol, fatol, **unknown_options) 546 xbar = numpy.add.reduce(sim[:-1], 0) / N 547 xr = (1 + rho) * xbar - rho * sim[-1] --> 548 fxr = func(xr) 549 doshrink = 0 550 /usr/local/lib/python3.4/dist-packages/scipy/optimize/optimize.py in function_wrapper(*wrapper_args) 290 def function_wrapper( *wrapper_args): 291 ncalls[0] += 1 --> 292 return function(*(wrapper_args + args)) 293 294 return ncalls, function_wrapper in loglikelihood(param) 28 for i in range(0, len(data)): 29 start_time = time.time() ---> 30 LL = LL - np.log(findYpdf(data.ix[i,'TC'], data.ix[i,'tsigt'])) 31 #print("--- data",i,"is processed in %s seconds ---" % (time.time() - start_time)) 32 print("processed in %s seconds ---" % (time.time() - start_time)) in findYpdf(x, tsigt) 20 def findYpdf(x, tsigt): 21 index = bisect(grid, tsigt) ---> 22 pdf_l = Y[index-1].pdf(x) 23 pdf_u = Y[index].pdf(x) 24 pdf_final = pdf_l + (pdf_u - pdf_l)*(tsigt - grid[index-1])/(grid[index] - grid[index-1]) #interpolate /usr/local/lib/python3.4/dist-packages/pacal/distr.py in pdf(self, x) 124 125 def pdf(self,x): --> 126 return self.get_piecewise_pdf()(x) 127 def cdf(self,x): 128 """Cumulative piecewise function.""" /usr/local/lib/python3.4/dist-packages/pacal/distr.py in get_piecewise_pdf(self) 78 """return PDF function as a PiecewiseDistribution object""" 79 if self.piecewise_pdf is None: ---> 80 self.init_piecewise_pdf() 81 return self.piecewise_pdf 82 def get_piecewise_cdf(self): /usr/local/lib/python3.4/dist-packages/pacal/distr.py in init_piecewise_pdf(self) 1126 return r1 + r2 1127 def init_piecewise_pdf(self): -> 1128 self.piecewise_pdf = conv(self.d1.get_piecewise_pdf(), self.d2.get_piecewise_pdf()) 1129 class SubDistr(SubRV, OpDistr): 1130 """Difference of distributions.""" /usr/local/lib/python3.4/dist-packages/pacal/distr.py in get_piecewise_pdf(self) 78 """return PDF function as a PiecewiseDistribution object""" 79 if self.piecewise_pdf is None: ---> 80 self.init_piecewise_pdf() 81 return self.piecewise_pdf 82 def get_piecewise_cdf(self): /usr/local/lib/python3.4/dist-packages/pacal/distr.py in init_piecewise_pdf(self) 1126 return r1 + r2 1127 def init_piecewise_pdf(self): -> 1128 self.piecewise_pdf = conv(self.d1.get_piecewise_pdf(), self.d2.get_piecewise_pdf()) 1129 class SubDistr(SubRV, OpDistr): 1130 """Difference of distributions.""" /usr/local/lib/python3.4/dist-packages/pacal/distr.py in get_piecewise_pdf(self) 78 """return PDF function as a PiecewiseDistribution object""" 79 if self.piecewise_pdf is None: ---> 80 self.init_piecewise_pdf() 81 return self.piecewise_pdf 82 def get_piecewise_cdf(self): /usr/local/lib/python3.4/dist-packages/pacal/distr.py in init_piecewise_pdf(self) 584 return self.f(self.d.rand(n, cache)) 585 def init_piecewise_pdf(self): --> 586 self.piecewise_pdf = self.d.get_piecewise_pdf().copyComposition(self.f, self.f_inv, self.f_inv_deriv, pole_at_zero = self.pole_at_zero) 587 588 class ShiftedScaledDistr(ShiftedScaledRV, OpDistr): /usr/local/lib/python3.4/dist-packages/pacal/distr.py in get_piecewise_pdf(self) 78 """return PDF function as a PiecewiseDistribution object""" 79 if self.piecewise_pdf is None: ---> 80 self.init_piecewise_pdf() 81 return self.piecewise_pdf 82 def get_piecewise_cdf(self): /usr/local/lib/python3.4/dist-packages/pacal/standard_distr.py in init_piecewise_pdf(self) 154 wrapped_pdf = wrap_pdf(self.pdf) 155 self.piecewise_pdf.addSegment(MInfSegment(b1, wrapped_pdf)) --> 156 self.piecewise_pdf.addSegment(Segment(b1, b2, wrapped_pdf)) 157 self.piecewise_pdf.addSegment(PInfSegment(b2, wrapped_pdf)) 158 def pdf(self, x): /usr/local/lib/python3.4/dist-packages/pacal/segments.py in addSegment(self, segment) 1016 if segment.isDirac(): 1017 segment = DiracSegment(seg.a, segment.f) -> 1018 assert seg>segment or segment>seg, "{}=={}=={}=={}=={}".format(self.segments, seg, segment, seg.a, segment.a) 1019 bisect.insort(self.segments, segment) 1020 self.breaks = unique(append(self.breaks, [segment.a, segment.b])) AssertionError: [MInfSegment, [-inf, 0.05648840406515476]]==MInfSegment, [-inf, 0.05648840406515476]==Segment, [0.05648840406515476, -0.08108631573253092]==-inf==0.05648840406515476 What does this error message mean? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#3>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABFcaY9oWd-B0PF7N4Q9J7utXLDdeZ-pks5t9XpXgaJpZM4Uqnxq> .

[michelleowen]

@jszymon thanks for the quick response.

I used the version for python 3.

Please see my code below:

import time
from bisect import bisect
from scipy.optimize import minimize
grid = np.linspace(0, 7.61659217e-01,num=100)

def loglikelihood(param):
    print("param:",param)
    start_time = time.time()
    #params = [sigma_bval, mu_tc, sigma_tc]
    X0 = NormalDistr(0, param[0])
    X1 = NormalDistr(param[1], param[2])
    X2 = exp(X1)
    X3 = X0 + X2

    Y = [None]*len(grid)
    for i in range(0, len(grid)):
        if (i==0):
            Y[i] = X3
        else:
            Y[i] = X3 + NormalDistr(0, grid[i])

    def findYpdf(x, tsigt):   
        index = bisect(grid, tsigt)
        pdf_l = Y[index-1].pdf(x)
        pdf_u = Y[index].pdf(x)
        #interpolate 
        pdf_final = pdf_l + (pdf_u - pdf_l)*(tsigt - grid[index-1])/(grid[index] - grid[index-1]) 
        return pdf_final 

    LL = 0
    for i in range(0, len(data)):
        #start_time = time.time()
        LL = LL - np.log(findYpdf(data.ix[i,'TC'], data.ix[i,'tsigt']))
        #print("--- data",i,"is processed in %s seconds ---" % (time.time() - start_time))
    print("--- processed in %s seconds ---" % (time.time() - start_time))
    print(LL)
    return LL

res = minimize(loglikelihood, [2.0,0,0.2], method='Nelder-Mead', tol=1e-6)

[michelleowen]

@jszymon I actually modified my code as follows, but the same error still raise.

import time
from bisect import bisect
from scipy.optimize import minimize
grid = np.linspace(0, 7.61659217e-01,num=100)

def loglikelihood(param):
    print("param:",param)
    start_time = time.time()
    #params = [sigma_bval, mu_tc, sigma_tc]
    X0 = NormalDistr(0, param[0])
    X1 = NormalDistr(param[1], param[2])
    X2 = exp(X1)
    X3 = X0 + X2

    Y = [None]*len(grid)
    for i in range(0, len(grid)):
        if (i==0):
            Y[i] = X3
        else:
            Y[i] = X3 + NormalDistr(0, grid[i])

    def findYpdf(x, tsigt):   
        try:
            index = grid.tolist().index(tsigt)
            pdf_final = Y[index].pdf(x)
        except:
            index = bisect(grid, tsigt)
            pdf_l = Y[index-1].pdf(x)
            pdf_u = Y[index].pdf(x)
            #interpolate 
            pdf_final = pdf_l + (pdf_u - pdf_l)*(tsigt - grid[index-1])/(grid[index] - grid[index-1]) 
        return pdf_final 

    LL = 0
    for i in range(0, len(data)):
        #start_time = time.time()
        LL = LL - np.log(findYpdf(data.ix[i,'TC'], data.ix[i,'tsigt']))
        #print("--- data",i,"is processed in %s seconds ---" % (time.time() - start_time))
    print("--- processed in %s seconds ---" % (time.time() - start_time))
    print(LL)
    return LL

res = minimize(loglikelihood, [2.0,0,0.2], method='Nelder-Mead', tol=1e-6)

[jszymon]

Hi, I used the version for python 3.

Atucally it has been merged to main branch, which is actually the newest version to use. Also, thanks for the code. Actually it does not run for me, since I don't have the 'data'. Tried to use some dummy data but it didn't work. There are 2 options: you can send me the data also, or if you don't want to do it, maybe you can send me the exact distributions with parameters and arguments in the function call which causes problems. If you choose the second option, a small functional program would make my life easier ;) All the best Szymon Please see my code below:

…

` import time from bisect import bisect from scipy.optimize import minimize grid = np.linspace(0, 7.61659217e-01,num=100) def loglikelihood(param): print("param:",param) start_time = time.time() #params = [sigma_bval, mu_tc, sigma_tc] X0 = NormalDistr(0, param[0]) X1 = NormalDistr(param[1], param[2]) X2 = exp(X1) X3 = X0 + X2 Y = [None]*len(grid) for i in range(0, len(grid)): if (i==0): Y[i] = X3 else: Y[i] = X3 + NormalDistr(0, grid[i]) def findYpdf(x, tsigt): index = bisect(grid, tsigt) pdf_l = Y[index-1].pdf(x) pdf_u = Y[index].pdf(x) #interpolate pdf_final = pdf_l + (pdf_u - pdf_l)*(tsigt - grid[index-1])/(grid[index] - grid[index-1]) return pdf_final LL = 0 for i in range(0, len(data)): #start_time = time.time() LL = LL - np.log(findYpdf(data.ix[i,'TC'], data.ix[i,'tsigt'])) #print("--- data",i,"is processed in %s seconds ---" % (time.time() - start_time)) print("--- processed in %s seconds ---" % (time.time() - start_time)) print(LL) return LL res = minimize(loglikelihood, [2.0,0,0.2], method='Nelder-Mead', tol=1e-6) ` — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#3 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABFcaRTEGYE8uDRfz5sFyo8o1vbqATjVks5t95eygaJpZM4Uqnxq> .

[michelleowen]

@jszymon
please see below:

import time
from bisect import bisect

params= [ 0.78440987, -0.01229896, -0.06878736]
grid = np.linspace(0, 7.61659217e-01,num=100)

tsigt = 0.10695393814651923
try:
    index = grid.tolist().index(tsigt)
except:
    index = bisect(grid, tsigt)


X0 = NormalDistr(0, params[0])
X1 = NormalDistr(params[1], params[2])
X2 = exp(X1)
X3 = X0 + X2

Y = [None]*len(grid)
for i in range(0, len(grid)):
    if (i==0):
        Y[i] = X3
    else:
        Y[i] = X3 + NormalDistr(0, grid[i])

pdf_l = Y[index-1].pdf(0.02049999999999841)
pdf_l

[jszymon]

Hi, it turned out to be quite simple. You have negative std. deviation for the normal distribution. You should use constrained optimization or some kind of parameter transformation. Pacal now prints a proper error message in this case. All the best, Szymon

…

On Tue, Jun 19, 2018 at 2:47 PM, michelleowen ***@***.***> wrote: @jszymon <https://github.com/jszymon> please see below: import time from bisect import bisect params= [ 0.78440987, -0.01229896, -0.06878736] grid = np.linspace(0, 7.61659217e-01,num=100) tsigt = 0.10695393814651923 try: index = grid.tolist().index(tsigt) except: index = bisect(grid, tsigt) X0 = NormalDistr(0, params[0]) X1 = NormalDistr(params[1], params[2]) X2 = exp(X1) X3 = X0 + X2 Y = [None]*len(grid) for i in range(0, len(grid)): if (i==0): Y[i] = X3 else: Y[i] = X3 + NormalDistr(0, grid[i]) pdf_l = Y[index-1].pdf(0.02049999999999841) pdf_l — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#3 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABFcaRrm0tEE2HWnjkXO2yE33snpiWw0ks5t-PLSgaJpZM4Uqnxq> .

[michelleowen]

@jszymon good catch!