#Estimating 3_D Gaussian distribution
import numpy as np
import lmfit
from lmfit import Model
x = np.linspace(-1., 1., 100)
y = x
z = x
c00 = 1. #parameters must have names
c01 = 0.
c02 = 0.
c11 = 1.
c12 = 0.
c22 = 1.
C = np.matrix(((c00, c01, c02), (c01, c11, c12), (c01, c12, c22)))
print(C) #just to verify, that it C the unit matric
#three dimensional Gaussian
def gaussian_nd(x, y, z, c00, c01, c02, c11, c12, c22):
C = np.array(((c00, c01, c02), (c01, c11, c12), (c01, c12, c22)))
X = np.array(((x), (y), (z)))
return np.array(1./np.sqrt(np.pi**3 / np.linalg.det(C)) * np.exp(X.T @ C @ X))
gauss = gaussian_nd(x, y, z, c00, c01, c02, c11, c12, c22) # true values
Z = gauss + np.random.normal(size=(100,100), scale=0.1) #measured values
print(Z.shape, type(Z))
model = lmfit.Model(gaussian_nd, independent_vars=['x', 'y', 'z'])
params = model.make_params()
params['c00'].set(value=1. + 0.01 * np.random.normal(size=1, scale=0.1))
params['c01'].set(value=0. + 0.01 * np.random.normal(size=1, scale=0.1))
params['c02'].set(value=0. + 0.01 * np.random.normal(size=1, scale=0.1))
params['c11'].set(value=1. + 0.01 * np.random.normal(size=1, scale=0.1))
params['c12'].set(value=0. + 0.01 * np.random.normal(size=1, scale=0.1))
params['c22'].set(value=1. + 0.01 * np.random.normal(size=1, scale=0.1))
result = model.fit(Z, x=x, y=y, z=z, params=params)
Now I get this error message:
LinAlgError: Last 2 dimensions of the array must be square