Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- from __future__ import division
- from scipy.optimize import *
- from scipy.misc import derivative
- from numpy import inf,linspace
- from scipy.integrate import quad,quadrature
- #from numpy.testing import *
- from csv import DictWriter
- #normalize transportation cost to t=1
- t=.5
- gu=.9
- gb=1.1
- q=[1.,1.]
- p=[.5,.5]
- def v(q,gamma,p):
- return max(q*gamma-q**2/2-p,0)
- def d0(gamma,q,p):
- if v(q[0],gamma,p[0])+v(q[1],gamma,p[1])<t:
- return v(q[0],gamma,p[0])/t
- elif v(q[0],gamma,p[0])-v(q[1],gamma,p[1])>t:
- return 1
- elif abs(v(q[0],gamma,p[0])-v(q[1],gamma,p[1]))<=t:
- return 1/2+(v(q[0],gamma,p[0])-v(q[1],gamma,p[1]))/(2*t)
- elif v(q[0],gamma,p[0])-v(q[1],gamma,p[1])<-t:
- return 0
- def d1(gamma,q,p):
- if v(q[0],gamma,p[0])+v(q[1],gamma,p[1])<t:
- return v(q[1],gamma,p[1])/t
- elif v(q[0],gamma,p[0])-v(q[1],gamma,p[1])>t:
- return 0
- elif abs(v(q[0],gamma,p[0])-v(q[1],gamma,p[1]))<=t:
- return 1/2+(v(q[1],gamma,p[1])-v(q[0],gamma,p[0]))/(2*t)
- elif v(q[0],gamma,p[0])-v(q[1],gamma,p[1])<-t:
- return 1
- def br(p,q):
- def pdf(gamma):
- return 1/(gb-gu)
- def pi0(p0):
- def integrand(gamma):
- return max(p0-q[0]**2/2,0)*d0(gamma,q,[p0,p[1]])
- value, error = quad(integrand,gu,gb)
- return -value
- def pi1(p1):
- def integrand(gamma):
- return max(p1-q[1]**2/2,0)*d1(gamma,q,[p[0],p1])
- value, error = quad(integrand,gu,gb)
- return -value
- pmax = v(q[0],gb,0)
- x0=p[0]
- p0 = float(fmin_l_bfgs_b(pi0,x0,bounds=[(0,pmax)],approx_grad=1)[0])
- print fmin_l_bfgs_b(pi0,x0,approx_grad=1),pi0(p[0])
- x0=p[1]
- p1 = float(fmin_l_bfgs_b(pi1,x0,bounds=[(0,pmax)],approx_grad=1)[0])
- print fmin_l_bfgs_b(pi1,x0,approx_grad=1),pi1(p[1])
- return [p0,p1]
- def eqP(p,q):
- print fixed_point(br,x0=p,args=(q,))
- return fixed_point(br,x0=p,args=(q,))
- print eqP(p,q)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement