;*****************************************************
;  Please sent only good vars into this subroutine
;*****************************************************

pro linear_fit_corr_bias,x,y,true_vector,$
  linear_slope,linear_intercept,lcor,bias,bias_sdv,$
  norm_dev

;  Sort them into increasing x to improve linear fit
xx=x[sort(x)]
yy=y[sort(x)]

;  Linear fit
linear_params=linfit(xx,yy)
linear_slope=linear_params[1]
linear_intercept=linear_params[0]

;  Linear correlation
lcor=correlate(x,y)

;  Bias of an estimator (from wikipedia)
;  In statistics, the difference between an estimator's 
;  expected value and the true value of the parameter 
;  being estimated is called the bias. An estimator 
;  having nonzero bias is said to be biased.
;  Suppose we are trying to estimate the parameter T 
;  using an estimator E (that is, some function of the 
;  observed data). Then the bias of E is defined to be
;  E-T.  In words, this would be "the expected value 
;  of the difference between the estimator E minus the 
;  true value T." 
if true_vector eq 'x' then begin
  bias_var=y-x
endif else if true_vector eq 'y' then begin
  bias_var=x-y
endif
bias_stats=moment(bias_var)
bias=bias_stats[0]
bias_sdv=sqrt(bias_stats[1])

;----------------------
;Normal Deviation (norm_dev)
;----------------------
residual_array = (linear_intercept + linear_slope * xx) - yy

residual_sq_sums = total(residual_array^2)

norm_dev = (residual_sq_sums/n_elements(xx))^(0.5)


return
end