Check out the Pastebin Gadgets Shop. We have thousands of fun, geeky & affordable gadgets on sale :-)Want more features on Pastebin? Sign Up, it's FREE!
tweet

# Untitled

By: a guest on Sep 4th, 2012  |  syntax: None  |  size: 1.84 KB  |  views: 18  |  expires: Never
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
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]])
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])
54.     return -value
55.   pmax = v(q[0],gb,0)
56.   x0=p[0]