• API
• FAQ
• Tools
• Trends
• Archive
daily pastebin goal
76%
SHARE
TWEET

# Untitled

a guest Sep 4th, 2012 20 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. from __future__ import division
2.
3. from scipy.optimize import *
4. from scipy.misc import derivative
5. from numpy import inf,linspace
7. #from numpy.testing import *
8. from csv import DictWriter
9.
10. #normalize transportation cost to t=1
11. t=.5
12. gu=.9
13. gb=1.1
14.
15. q=[1.,1.]
16. p=[.5,.5]
17.
18. def v(q,gamma,p):
19.   return max(q*gamma-q**2/2-p,0)
20.
21. def d0(gamma,q,p):
22.   if v(q[0],gamma,p[0])+v(q[1],gamma,p[1])<t:
23.     return v(q[0],gamma,p[0])/t
24.   elif v(q[0],gamma,p[0])-v(q[1],gamma,p[1])>t:
25.     return 1
26.   elif abs(v(q[0],gamma,p[0])-v(q[1],gamma,p[1]))<=t:
27.     return 1/2+(v(q[0],gamma,p[0])-v(q[1],gamma,p[1]))/(2*t)
28.   elif v(q[0],gamma,p[0])-v(q[1],gamma,p[1])<-t:
29.     return 0
30.
31. def d1(gamma,q,p):
32.   if v(q[0],gamma,p[0])+v(q[1],gamma,p[1])<t:
33.     return v(q[1],gamma,p[1])/t
34.   elif v(q[0],gamma,p[0])-v(q[1],gamma,p[1])>t:
35.     return 0
36.   elif abs(v(q[0],gamma,p[0])-v(q[1],gamma,p[1]))<=t:
37.     return 1/2+(v(q[1],gamma,p[1])-v(q[0],gamma,p[0]))/(2*t)
38.   elif v(q[0],gamma,p[0])-v(q[1],gamma,p[1])<-t:
39.     return 1
40.
41.
42. def br(p,q):
43.   def pdf(gamma):
44.     return 1/(gb-gu)
45.   def pi0(p0):
46.     def integrand(gamma):
47.       return max(p0-q[0]**2/2,0)*d0(gamma,q,[p0,p[1]])
48.     value, error = quad(integrand,gu,gb)
49.     return -value
50.   def pi1(p1):
51.     def integrand(gamma):
52.       return max(p1-q[1]**2/2,0)*d1(gamma,q,[p[0],p1])
53.     value, error = quad(integrand,gu,gb)
54.     return -value
55.   pmax = v(q[0],gb,0)
56.   x0=p[0]
57.   p0 = float(fmin_l_bfgs_b(pi0,x0,bounds=[(0,pmax)],approx_grad=1)[0])