def error_ellipse(mue1,mue2,sig1,sig2,rho,sample_size):

  """creates correlated binormal distributed variates from
  uncorrelated ones, plots a corresponding sample and
  overplots the 1,2 and 3 sigma error-ellipses.
  input: mu_x, mu_y, sig_x, sig_y, corr. coeff, sample_size
  note: sig_x and sig_y refer to the correlated variates
  """
  import numpy as np
  import matplotlib.pyplot as plt
  import matplotlib.patches as patches

# values for the CORRELATED variables  
#  mue1=2.
#  mue2=2.
#  sig1=1.
#  sig2=np.sqrt(2.)
  
# we want to have ss x and ss y coordinates
  ss=2*sample_size 

#  np.random.seed(5)

  dist=np.random.normal(size=ss)
  
  dist1=dist[0:ss/2]
  dist2=dist[ss/2:ss]
  
  if sig1 != sig2:
    tan2theta=2.*sig1*sig2*rho/(sig1**2-sig2**2)  
  else:
    tan2theta=1.e20
  
  theta=np.arctan(tan2theta)/2.
  print 'rotation angle',theta/np.pi*180.
  
  sig1_2=sig1**2
  sig2_2=sig2**2
  
  sinth=np.sin(theta)
  costh=np.cos(theta)
  
# values for the UNCORRELATED variables  
  sig1p_2=(sig1_2*sig2_2*(1.-rho**2)/
    (sig1_2*sinth**2+sig2_2*costh**2-2.*rho*sig1*sig2*sinth*costh))
  
  sig2p_2=(sig1_2*sig2_2*(1.-rho**2)/
    (sig1_2*costh**2+sig2_2*sinth**2+2.*rho*sig1*sig2*sinth*costh))

  sig1p=np.sqrt(sig1p_2)  
  sig2p=np.sqrt(sig2p_2)  

  print 'sig_x(uncorrelated):',sig1p
  print 'sig_y(uncorrelated):',sig2p
  
  dist1=dist1*sig1p
  dist2=dist2*sig2p

#transform uncorrelated variates to correlated ones, but
#orthogonal rotation. The rotation matrix corresponds to the
#inverse(=transposed) of the one given in the lecture.
  
  x=costh*dist1-sinth*dist2
  y=sinth*dist1+costh*dist2
  
  x=x+mue1
  y=y+mue2
  
  plt.clf()
  
  plt.ion()
  
  ax = plt.subplot(111, aspect='equal')
  ax.set_xlim(mue1-4*sig1,mue1+4*sig1)
  ax.set_ylim(mue2-4*sig2,mue2+4*sig2)

#plot error-ellipses with axis (lengths = 2*sig(uncorrelated)) and rotation angle from above
#overplot corresponding error-rectangle (lengths = 2*sig(correlated))
  for i in xrange(1,4):
# i =1,3 corresponds to 1 to 3 sig error ellipses
    ells = patches.Ellipse(xy=(mue1,mue2), width=2*i*sig1p, height=2*i*sig2p, angle=theta/np.pi*180, 
      fill=0,linewidth=2,edgecolor='red')
    ax.add_patch(ells)
    rect = patches.Rectangle((mue1-i*sig1,mue2-i*sig2), 2*i*sig1, 2*i*sig2,
      fill=0,linestyle='dotted',edgecolor='black')
    ax.add_patch(rect)

#overplot correlated random-variates
  plt.scatter(x,y,s=2,c='blue',edgecolors='blue') #s=size
  plt.show()