Pastebin is 300% more awesome when you are logged in. Sign Up, it's FREE!
Guest

Untitled

By: a guest on Sep 4th, 2012  |  syntax: None  |  size: 1.84 KB  |  hits: 18  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
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
  6. from scipy.integrate import quad,quadrature
  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])
  58.   print fmin_l_bfgs_b(pi0,x0,approx_grad=1),pi0(p[0])
  59.   x0=p[1]
  60.   p1 = float(fmin_l_bfgs_b(pi1,x0,bounds=[(0,pmax)],approx_grad=1)[0])
  61.   print fmin_l_bfgs_b(pi1,x0,approx_grad=1),pi1(p[1])
  62.   return [p0,p1]
  63.  
  64. def eqP(p,q):
  65.   print fixed_point(br,x0=p,args=(q,))
  66.   return fixed_point(br,x0=p,args=(q,))
  67.  
  68. print eqP(p,q)