0

I am trying to surface fit 3d data (z is a function of x and y). I have assymetrical error bars for each point. I would like the fit to take this uncertainty into account.

I am using scipy.linalg.lstsq(). It does not have any option for uncertainties in its arguments.

I am trying to adapt some code found on this page.

import numpy as np
import scipy.linalg
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

# Create data with x and y random over [-2, 2], and z a Gaussian function of x and y.
np.random.seed(12345)
x = 2 * (np.random.random(500) - 0.5)
y = 2 * (np.random.random(500) - 0.5)

def f(x, y):
    return np.exp(-(x + y ** 2))

z = f(x, y)

data = np.c_[x,y,z]

# regular grid covering the domain of the data
mn = np.min(data, axis=0)
mx = np.max(data, axis=0)
X,Y = np.meshgrid(np.linspace(mn[0], mx[0], 20), np.linspace(mn[1], mx[1], 20))
XX = X.flatten()
YY = Y.flatten()

# best-fit quadratic curve (2nd-order)
A = np.c_[np.ones(data.shape[0]), data[:,:2], np.prod(data[:,:2], axis=1), data[:,:2]**2]
C,_,_,_ = scipy.linalg.lstsq(A, data[:,2])

# evaluate it on a grid
Z = np.dot(np.c_[np.ones(XX.shape), XX, YY, XX*YY, XX**2, YY**2], C).reshape(X.shape)

# plot points and fitted surface using Matplotlib
fig = plt.figure(figsize=(10, 10))
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=0.2)
ax.scatter(data[:,0], data[:,1], data[:,2], c='r', s=50)
plt.xlabel('X')
plt.ylabel('Y')
ax.set_zlabel('Z')
ax.axis('equal')
ax.axis('tight')
Tommy95
  • 21
  • 4
  • 2
    Hi, did you read the docu? **scipy.linalg.lstsq...Compute least-squares solution to equation Ax = b.** So, I'd say it is not really suited for your exponential. – mikuszefski Jun 24 '19 at 12:35
  • Also have a look here: https://stackoverflow.com/q/19116519/803359 – mikuszefski Jun 24 '19 at 12:39
  • ... thanks but change it to z=x² + y² if you prefer – Tommy95 Jun 24 '19 at 12:39
  • Least squares fitting, like any other tool, works if your data satisfies, or at least almost satisfies certain assumptions. Your data does not, so least squares is not the crowbar for your nail. – Mad Physicist Jun 24 '19 at 12:39
  • Thanks Mikuszefski. Thanks Mad Physicist but sorry, how do you know my data does not satisfy certain assumptions? Which assumptions are you referring to? – Tommy95 Jun 24 '19 at 12:43
  • 1
    Well, you don't have a Gaussian to fit, you have random data with a Gaussian distribution. – mikuszefski Jun 24 '19 at 12:45
  • in 1D this can be handled by scipy.stats.norm.fit` but not in 2D – mikuszefski Jun 24 '19 at 12:46
  • This data is simulated for the purpose of this question. And if you actually run the code, the fit is perfect. Which other fitting methodology do you believe would be more appropriate to use for asteroseismology (fmin and fmax in x and y, based on an fft, and temperature in z) ? – Tommy95 Jun 24 '19 at 12:48
  • 1
    Sorry, got it wrong, so you really build a Gaussian over random (x,y), OK. if you want to fit this with asymmetric errors you could try to build your own cost funvtion and use `scipy.optimize.leastsq` – mikuszefski Jun 24 '19 at 12:52
  • scipy.optimize.curvefit has the possibility of using error bars so I might use that instead of scipy.linalg.lstsq – Tommy95 Jun 24 '19 at 13:12
  • @Tommy95. That's a much better idea. – Mad Physicist Jun 24 '19 at 14:54

0 Answers0