Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

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

Closed
michelleowen opened this issue Jun 16, 2018 · 7 comments
Closed

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

michelleowen opened this issue Jun 16, 2018 · 7 comments

Comments

@michelleowen
Copy link

michelleowen commented Jun 16, 2018 •
edited
Loading

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
@jszymon
Copy link
Owner

jszymon commented Jun 17, 2018 via email

@michelleowen
Copy link
Author

michelleowen commented Jun 18, 2018 •
edited
Loading

@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
Copy link
Author

@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
Copy link
Owner

jszymon commented Jun 19, 2018 via email

@michelleowen
Copy link
Author

@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
Copy link
Owner

jszymon commented Jun 20, 2018 via email

@michelleowen
Copy link
Author

@jszymon good catch!

@jszymon jszymon closed this as completed Jul 2, 2018
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@jszymon @michelleowen