title: Estimate the standard deviation of fitted parameters in scipy.odr?
body: <p>(<em>Somewhat related to this question <a href="https://stackoverflow.com/questions/23951876/linear-fit-including-all-errors-with-numpy-scipy">Linear fit including all errors with NumPy/SciPy</a>, and borrowing code from this one <a href="https://stackoverflow.com/questions/22670057/linear-fitting-in-python-with-uncertainty-in-both-x-and-y-coordinates">Linear fitting in python with uncertainty in both x and y coordinates</a></em>)</p>
  
  <p>I fit a linear model (<code>y=a*x+b</code>) using fixed errors in <code>x,y</code> using <a href="https://docs.scipy.org/doc/scipy/reference/odr.html" rel="nofollow noreferrer">scipy.odr</a> (code is below), and I get:</p>
  
  <pre><code>Parameters (a, b): [ 5.21806759 -4.08019995]
  Standard errors: [ 0.83897588  2.33472161]
  Squared diagonal covariance: [ 1.06304228  2.9582588 ]
  </code></pre>
  
  <p>What is the correct standard deviation values for the fitted <code>a, b</code> parameters? I'm assuming these must be obtained from the <code>Squared diagonal covariance</code> values, but then how are these values related to the <code>Standard errors</code>?</p>
  
  <p><strong>ADD</strong></p>
  
  <p>As mentioned in the answer to <a href="https://stackoverflow.com/questions/41028846/how-to-compute-standard-error-from-odr-results">How to compute standard error from ODR results?</a> by <code>ali_m</code>, this is apparently related to <a href="https://github.com/scipy/scipy/issues/6842" rel="nofollow noreferrer">a bug in scipy.odr</a>. If one uses</p>
  
  <pre><code>np.sqrt(np.diag(out.cov_beta * out.res_var))
  </code></pre>
  
  <p>(i.e: multiply the covariance by the residual variance) instead of just</p>
  
  <pre><code>np.sqrt(np.diag(out.cov_beta))
  </code></pre>
  
  <p>the result now coincides with <code>out.sd_beta</code>.</p>
  
  <p>So now my question is: which is the proper standard deviation for the fitted parameters <code>(a, b)</code>? Is it <code>out.sd_beta</code> (equivalently: <code>np.sqrt(np.diag(out.cov_beta * out.res_var))</code>) or <code>np.sqrt(np.diag(out.cov_beta))</code>?</p>
  
  <p><br><br><br></p>
  
  <hr>
  
  <pre><code>import numpy as np
  from scipy.odr import Model, RealData, ODR
  import random
  
  random.seed(9001)
  np.random.seed(117)
  
  def getData(c):
      """Initiate random data."""
      x = np.array([0, 1, 2, 3, 4, 5])
      y = np.array([i**2 + random.random() for i in x])
      xerr = c * np.array([random.random() for i in x])
      yerr = c * np.array([random.random() for i in x])
      return x, y, xerr, yerr
  
  def linear_func(p, x):
      """Linear model."""
      a, b = p
      return a * x + b
  
  def fitModel(x, y, xerr, yerr):
      # Create a model for fitting.
      linear_model = Model(linear_func)
      # Create a RealData object using our initiated data from above.
      data = RealData(x, y, sx=xerr, sy=yerr)
      # Set up ODR with the model and data.
      odr = ODR(data, linear_model, beta0=[0., 1.])
      # Run the regression.
      out = odr.run()
  
      # Estimated parameter values
      beta = out.beta
      print("Parameters (a, b): {}".format(beta))
      # Standard errors of the estimated parameters
      std = out.sd_beta
      print("Standard errors: {}".format(std))
      # Covariance matrix of the estimated parameters
      cov = out.cov_beta
      stddev = np.sqrt(np.diag(cov))
      print("Squared diagonal covariance: {}".format(stddev))
  
  
  # Generate data and fit the model.
  x, y, xerr, yerr = getData(1.)
  fitModel(x, y, xerr, yerr)
  </code></pre>