title: How to perform sigma clipping on data set?
body: <p>*Context: I am looking at the change of velocity of an object periodically, where the period is 1.846834days and I am expecting a sinusoidal fit to my set of data.
  Suppose I have a set of data which look like this:</p>
  
  <pre><code>#days        vel         error
  5725.782701 0.195802    0.036312
  5729.755560 -0.006370   0.041495
  5730.765352 -0.071253   0.030760
  5745.710214 0.092082    0.036094
  5745.932853  0.238030   0.040097
  5749.705307 0.196649    0.037140
  5741.682112 0.186664    0.028075
  5742.681765 -0.262104   0.038049
  6186.729146 -0.243796   0.031687
  6187.742803 -0.009394   0.054541
  6190.717317 -0.001821   0.033684
  6192.716356 0.117557    0.037807
  6196.704736 0.093935    0.032336
  6203.683879 0.076051    0.033085
  6204.679898 -0.301463   0.033483
  6208.659585 -0.409340   0.036002
  6209.669701 0.180807    0.041666
  </code></pre>
  
  <p>There are only one or two data points observed in each cycle, so what I wanted to do is to phase fold my data, plot them and fit my data using chi-square minimization. This is what I have so far:</p>
  
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  from scipy import optimize as op
  import csv
  
  with open('testdata.dat') as fin:
    with open('testdata_out.dat', 'w') as fout:   #write output in a new .dat file
      o=csv.writer(fout)
      for line in fin:
          o.writerow(line.split())
  
  
  #create a 2d array and give names to columns
  data = np.genfromtxt('testdata_out.dat',delimiter=',',dtype=('f4,f4,f4'))
  data.dtype.names = ('bjd','rv','err')
  
  
  # parameters
  x = data['bjd']
  
  y = data['rv']
  
  e = data['err']
  
  P = 1.846834    #orbital period
  
  T = 3763.85112  #time from ephemeris
  
  q = (data['bjd']-T)%P
  
  #print(q)
  
  
  def niceplot():
      plt.xlabel('BJD')
      plt.ylabel('RV (km/s)')
      plt.tight_layout()
  
  def model(q,A,V):
      return A*np.cos(np.multiply(np.divide((q),P),np.multiply(2,np.pi))-262) + V
  
  def residuals((A,V),q,y,e):     #for least square
      return (y-model(q,A,V)) / e
  
  def chisq((A,V),q,y,e):
      return ((residuals((A,V),q,y,e))**2).sum()
  
  result_min = op.minimize(chisq, (0.3,0), args=(q,y,e))
  print(result_min)
  
  
  A,V = result_min.x
  xf = np.arange(q.min(), q.max(), 0.1)
  yf = model(xf,A,V)
  
  print(xf)
  print(yf)
  
  plt.errorbar(q, y, e, fmt='ok')
  plt.plot(xf,yf)
  niceplot()
  plt.show()
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/yVyHm.jpg" alt="Chi-squared minimized fit"></p>
  
  <p>In my plot, the shape of the sine curve seems to fit my data, but it does not go through all the data points. 
  My question is: How can I perform sigma clipping such that I can get a better fit for my set of data? I know the curve_fit module in Scipy may do the job. But I would like to know if it is possible to perform sigma clipping while using minimize module??</p>
  
  <p>Many thanks in advance!</p>