Curve_fit 用于返回 numpy 数组的函数（sumproduct函数返回value）

25-04-28 1

针对Curve_fit用于返回numpy数组的函数和sumproduct函数返回value这两个问题，本篇文章进行了详细的解答，同时本文还将给你拓展2DGaussianCurve_fit与有缺陷的数据

针对Curve_fit 用于返回 numpy 数组的函数和sumproduct函数返回value这两个问题，本篇文章进行了详细的解答，同时本文还将给你拓展2D Gaussian Curve_fit 与有缺陷的数据、Curve_Fit 不准确、numpy.random.random & numpy.ndarray.astype & numpy.arange、numpy.ravel()/numpy.flatten()/numpy.squeeze()等相关知识，希望可以帮助到你。

本文目录一览：

Curve_fit 用于返回 numpy 数组的函数（sumproduct函数返回value）
2D Gaussian Curve_fit 与有缺陷的数据
Curve_Fit 不准确
numpy.random.random & numpy.ndarray.astype & numpy.arange
numpy.ravel()/numpy.flatten()/numpy.squeeze()

Curve_fit 用于返回 numpy 数组的函数（sumproduct函数返回value）

如何解决Curve_fit 用于返回 numpy 数组的函数

我知道 curve_fit 的库 scipy 及其拟合曲线的能力。我在这里和文档中阅读了许多示例，但我无法解决我的问题。例如，我有 10 个文件（化学结构，但无所谓）和 10 个实验能量值。我在一个类中有一个函数，它为每个结构计算某些参数的理论能量，并返回一个包含理论能量值的 numpy 数组。

我想找到理论值最接近实验值的最佳参数。我将在这里提供我的代码的最小示例

这是读取实验能量文件、提取正确子串并将值作为numpy数组返回的类函数。 self.path 只是目录和 self.nPoints = 10。这不是那么重要，但为了完整起见我提供

def experimentalValues(self):
        os.chdir(self.path)
        energy = np.zeros(self.nPoints)
        for i in range(1,self.nPoints):
            f = open("p_" + str(i + 1) + ".xyz","r")
            energy[i] = float(f.readlines()[1].split()[1])
            f.close()
        os.chdir(''..'')
        return energy

我用这个以两个 numpy 数组作为参数的类函数计算了理论值，比如说

sigma = np.full(nSubstrate,2.)
epsilon = np.full(nSubstrate,0.15)

哪里nSubstrate = 9

这里是类函数。它读取文件并执行两个嵌套循环来为每个文件计算理论值并将其返回到一个 numpy 数组。

def theoreticalEnergy(self,epsilon,sigma):
        os.chdir(self.path)
        cE = np.zeros(self.nPoints)
        for n in range(0,self.nPoints):
            filenameXYZ = "p_" + str(n + 1) + "_extended.xyz"

            allCoordinates = np.loadtxt(filenameXYZ,skiprows = 0,usecols = (1,2,3))
            substrate = allCoordinates[0:self.nSubstrate]
            surface = allCoordinates[self.nSubstrate:]
            for i in range(0,substrate.shape[0]):
                positionAtomI = np.array(substrate[i][:])
                for j in range(0,surface.shape[0]):
                    positionAtomJ = np.array(surface[j][:])
                    distanceIJ = self.distance(positionAtomI,positionAtomJ)
                    cE[n] += self.LennardJones(distanceIJ,epsilon[i],sigma[i])
                
        os.chdir(''..'')
        return cE

同样，为了完整起见，Lennard Jones 类函数定义为

def LennardJones(self,distance,sigma):
        repulsive = (sigma/distance) ** 12.
        attractive = (sigma/distance) ** 6.
        potential = 4. * epsilon* (repulsive - attractive)
        return potential

在这种情况下，所有参数都是标量作为返回值。总结问题陈述，我有 3 个要素：

一个带有实验数据的 numpy 数组
两个 numpy 数组，猜测参数 sigma 和 epsilon
一个函数，它接受最后一个参数并返回一个包含要拟合的值的 numpy 向量。

如何像文档 https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html 中描述的方法一样解决这个问题？

解决方法

曲线拟合

curve_fit 通过找到使 f(w,x[i]) 最小化的 y[i] 使函数 w 拟合到点 sum((f(w,x[i] - y[i])**2 for i in range(n))。正如您将在函数定义后的第一行中读到的那样

[它使用] 非线性最小二乘法将函数 f 拟合到数据。

它指的是 least_squares 的地方

给定残差 f(x)（n 个实变量的 mD 实函数）和损失函数 rho(s)（标量函数），least_squares 找到成本函数的局部最小值 F(x)：

曲线拟合是一种凸成本多目标优化。由于每个单独的成本都是凸的，您可以将所有成本相加，这仍然是一个凸函数。请注意，决策变量（要优化的参数）在每个点都是相同的。

你的问题

根据我对每个能级的理解，你有一组不同的参数，如果你把它写成曲线拟合问题，目标函数可以表示为 sum((f(w[i],x[i]) - y[i])**2 ...),where y[i]is determined by the energy level. Since each of the terms in the sum is independent on the other terms,this is equivalent to finding each group of parametersw [i]separately minimizing(f(w[i],x[i]) - y[i])**2`。

凸性

凸性是一种非常方便的优化属性，因为它确保参数空间中只有一个最小值。我没有做详细的分析，但对你的能量函数的凸性有合理的怀疑。

Lennard Jones 函数具有排斥力和吸引力之差，两者在距离上都具有负偶数指数，仅此一项不太可能是凸的。
以不同位置为中心的多个局部函数的和没有定义的凸性。
众所周知，分子能量、晶体能量或蛋白质折叠是非凸的。

几天前（骑自行车）我在想这个问题，分子将如何在全局最小能量下配置，我想知道它是否会因为量子隧道效应而如此迅速地找到这种配置。

非凸优化

非凸（全局）优化与（非线性）最小二乘法不同，从某种意义上说，当找到局部最小值时，过程不会立即返回，而是开始在不同区域进行新的尝试搜索空间。如果函数是平滑的，你仍然可以利用基于梯度的局部优化方法，但复杂度仍然是 NP。

一个经典的全局优化方法是 Simulated annenaling，如果你有化学背景，我想你会对它有一些见解。 Once upon a time，在 scipy.optimize 中提供了模拟退火。

您会在 scipy.optimize 中找到一些 global optimization 方法。我鼓励您尝试 Basin hopping，因为它已成功应用于类似的问题，您可以在 references 中阅读。

我希望这能让您找到解决方案的正确方法。但是，请注意，您可能需要花时间学习如何使用该功能，并且需要做出一些决定。您需要在准确性、简单性和效率之间找到平衡点。

如果您想要更好的解决方案，请花时间推导出成本函数的梯度（您可以返回两个值 f 和 df，其中 df 是 f 相对于决策变量的梯度）。

2D Gaussian Curve_fit 与有缺陷的数据

如何解决2D Gaussian Curve_fit 与有缺陷的数据

我一直在尝试在我的数据上拟合 2D 高斯分布。我有一个 (x,y) 地图和每个坐标的 z 值。但是，我想删除一些 z 值，但是由于 curve_fit 要求将数组作为 (x,y) 参数，因此很难删除 (x,y) 数组中的对应 (x,y) 值。

我已经看到有很多解决方案（掩码数组、lmfit），但我选择了一个更简单的解决方案，因为我对 python 不太擅长。

这是我的代码：

def gaussienne(Var,xo,yo,backgr,gx,gy,theta):
    
    x,y = Var
    
    a = np.cos(theta)**2/(2*gx**2) + np.sin(theta)**2/(2*gy**2)
    b = -np.sin(2*theta)/(4*gx**2) + np.sin(2*theta)/(4*gy**2)
    c = np.sin(theta)**2/(2*gx**2) + np.cos(theta)**2/(2*gy**2)
    
    res=backgr+np.exp(-(a*(x-xo)**2 + 2*b*(x-xo)*(y-yo)+ c*(y-yo)**2))
    
    return res.ravel() 
x=np.linspace(start_x,start_x+(taille_x)*step_x,taille_x)
y=np.linspace(start_y,start_y+(taille_y)*step_y,taille_y)
mask=[]
intensite=intensite.reshape(taille_x,taille_y)
for i in range(taille_x):
     for j in range (taille_y):
         if intensite[i,j]!=0:
               mask.append([x[i],y[j],intensite[i,j]])
mask=np.array(mask)
p0 = a,b,80,0.2,0.001
popt,pcov = curve_fit(gaussienne,mask[:,:2],2],p0,maxfev=5000)

基本上，对于我想要取的每个 z 值（即不等于 0），我添加以屏蔽 z 值和对应的 (x,y) 坐标。掩码是一个元组列表。我把它变成一个数组并将它交给curve_fit。

然后我得到错误：ValueError：解包的值太多（预期为 2）

我不明白我是怎么得到这个错误的，因为 mask[:,:2] 是 (2,n) 而 mask[:,2] 是 (1,n)

感谢您的帮助

解决方法

我不明白我是怎么得到这个错误的，因为 mask[:,:2] 是 (2,n)

不完全 - mask[:,:2].shape 是 (n,2)，所以它必须是 x,y = Var.T - 或者当然你可以将 mask[:,:2].T 传递给 curve_fit()。

Curve_Fit 不准确

如何解决Curve_Fit 不准确

我试图尽可能好地拟合随时间变化很大的数据。所以首先我平滑了工作正常的数据。我从中得到的平滑数据应该进一步从拟合中表示出来，以得到更多的峰值。正如您在代码中看到的，我想使用 log-tanh 函数来拟合数据。我很清楚这个问题已经在一些线程中出现了，但我已经尝试过，而且数据也不是很小或很大，我知道这也会导致问题。

如您所见，我尝试的多项式拟合也很有效，但它并没有消除所有波浪值。它们对下面的导数造成了非常糟糕的问题。

import tkinter as tk
from tkinter import filedialog
import numpy as np
import scipy.signal
from scipy.optimize import curve_fit
from numpy import diff
import matplotlib.pyplot as plt
from lmfit.models import StepModel,LinearModel



def loghypfunc(x,A,B,C,D,E):
    return A*np.log(1+x)+B*np.tanh(C*x)+D*x+E

def expfunc(t,c0,c1,c2,c3):
    return c0+c1*t-c2*np.exp(-c3*t)

def expdecay(x,a,b,c):
    return a * np.exp(-b * x) + c



path="C:/Users/Sammy/Documents/Masterarbeit WT/CSM und Kriechdaten/Kriechen/Creep_10mN_00008_LC_20210406_2121_DYN.txt"

dataFile = np.loadtxt(path,delimiter=''\\t'',skiprows=2,usecols=(0,1,2,3,29,30),dtype=float)
num_rows,num_cols = dataFile.shape

# time column
time = dataFile[:,[0]].transpose()
time = time.flatten()
refTime = time[0]  # get first time in column (reference)
# genullte Testzeit
timeNull = time - refTime
print("time",time)
flatTimeNull = timeNull.flatten()  # jetzt ein 1D array (one row)
##################################################################################
# indent displacement column
indentdis = dataFile[:,[4]].transpose()
indentdis = indentdis.flatten()
indentdis = indentdis - indentdis[0]
# the indendt data has to be smoothed so there is not such a big fluctuation
indentSmooth = scipy.signal.savgol_filter(indentdis,2001,3)
# null the indent Smooth data
indentSmooth_Null = indentSmooth - indentSmooth[0]
hind_Smooth_flat = indentSmooth_Null.flatten()  # jetzt ein 1D array
print(''indent smooth'',indentSmooth)
######################################################################

p0 = [100,0.1,100,0.1]
c,cov = curve_fit(expfunc,time,indentSmooth,p0)
y_indent = expfunc(indentSmooth,*c)


p0 = [70,0.5,50,100]
popt,pcov = curve_fit(loghypfunc,p0,maxfev = 10000)
y_indentTan = loghypfunc(indentSmooth,*popt)

modelh_t = np.poly1d(np.polyfit(time,8))

plt.plot(time,''r'',label="Data smoothed")
plt.scatter(time,modelh_t(time),s=0.1,label="polyfit")
plt.plot(time,y_indentTan,label="Curve fit Tangens function")
plt.plot(time,y_indent,label="Curve fit exp function")
plt.legend(loc="lower right")
plt.xlabel("time")
plt.ylabel("indent")

plt.show()

这是我从中获取数据的两个数组

time [  6.299596   6.349592   6.399589 ... 608.0109   608.060897 608.110894] 
indent smooth [120.81411822 121.07093706 121.32748184 ... 476.78825661 476.89357473 476.99915287]

这里是情节 Plots

我现在的问题是如何解决它。是因为要拟合的错误优化参数吗？但是我猜python应该自动做到足够好？我的第二个猜测是数据被定时沿这个轴压缩，因为数组大约有 12000 个值。这可能是一个原因吗？

我将非常感谢任何有关合身的建议。

问候 Hndrx

numpy.random.random & numpy.ndarray.astype & numpy.arange

今天看到这样一句代码：

xb = np.random.random((nb, d)).astype(''float32'') #创建一个二维随机数矩阵（nb行d列）

xb[:, 0] += np.arange(nb) / 1000. #将矩阵第一列的每个数加上一个值

要理解这两句代码需要理解三个函数

1、生成随机数

numpy.random.random(size=None)

size为None时，返回float。

size不为None时，返回numpy.ndarray。例如numpy.random.random((1,2))，返回1行2列的numpy数组

2、对numpy数组中每一个元素进行类型转换

numpy.ndarray.astype(dtype)

返回numpy.ndarray。例如 numpy.array([1, 2, 2.5]).astype(int)，返回numpy数组 [1, 2, 2]

3、获取等差数列

numpy.arange([start,]stop,[step,]dtype=None)

功能类似python中自带的range()和numpy中的numpy.linspace

返回numpy数组。例如numpy.arange(3)，返回numpy数组[0, 1, 2]

numpy.ravel()/numpy.flatten()/numpy.squeeze()

numpy.ravel(a, order=''C'')

　　Return a flattened array

numpy.chararray.flatten(order=''C'')

　　Return a copy of the array collapsed into one dimension

numpy.squeeze(a, axis=None)

　　Remove single-dimensional entries from the shape of an array.

相同点: 将多维数组降为一维数组

不同点:

　　ravel() 返回的是视图(view)，意味着改变元素的值会影响原始数组元素的值;

　　flatten() 返回的是拷贝，意味着改变元素的值不会影响原始数组;

　　squeeze()返回的是视图(view)，仅仅是将shape中dimension为1的维度去掉;

ravel()示例:

 1 import matplotlib.pyplot as plt
 2 import numpy as np
 3 
 4 def log_type(name,arr):
 5     print("数组{}的大小:{}".format(name,arr.size))
 6     print("数组{}的维度:{}".format(name,arr.shape))
 7     print("数组{}的维度:{}".format(name,arr.ndim))
 8     print("数组{}元素的数据类型:{}".format(name,arr.dtype))
 9     #print("数组:{}".format(arr.data))
10     
11 a = np.floor(10*np.random.random((3,4)))
12 print(a)
13 log_type(''a'',a)
14 
15 a1 = a.ravel()
16 print("a1:{}".format(a1))
17 log_type(''a1'',a1)
18 a1[2] = 100
19 
20 print(a)
21 log_type(''a'',a)

flatten()示例

 1 import matplotlib.pyplot as plt
 2 import numpy as np
 3 
 4 def log_type(name,arr):
 5     print("数组{}的大小:{}".format(name,arr.size))
 6     print("数组{}的维度:{}".format(name,arr.shape))
 7     print("数组{}的维度:{}".format(name,arr.ndim))
 8     print("数组{}元素的数据类型:{}".format(name,arr.dtype))
 9     #print("数组:{}".format(arr.data))
10     
11 a = np.floor(10*np.random.random((3,4)))
12 print(a)
13 log_type(''a'',a)
14 
15 a1 = a.flatten()
16 print("修改前a1:{}".format(a1))
17 log_type(''a1'',a1)
18 a1[2] = 100
19 print("修改后a1:{}".format(a1))
20 
21 print("a:{}".format(a))
22 log_type(''a'',a)

squeeze()示例：

1. 没有single-dimensional entries的情况

 1 import matplotlib.pyplot as plt
 2 import numpy as np
 3 
 4 def log_type(name,arr):
 5     print("数组{}的大小:{}".format(name,arr.size))
 6     print("数组{}的维度:{}".format(name,arr.shape))
 7     print("数组{}的维度:{}".format(name,arr.ndim))
 8     print("数组{}元素的数据类型:{}".format(name,arr.dtype))
 9     #print("数组:{}".format(arr.data))
10     
11 a = np.floor(10*np.random.random((3,4)))
12 print(a)
13 log_type(''a'',a)
14 
15 a1 = a.squeeze()
16 print("修改前a1:{}".format(a1))
17 log_type(''a1'',a1)
18 a1[2] = 100
19 print("修改后a1:{}".format(a1))
20 
21 print("a:{}".format(a))
22 log_type(''a'',a)

从结果中可以看到，当没有single-dimensional entries时，squeeze()返回额数组对象是一个view,而不是copy。

2. 有single-dimentional entries 的情况

 1 import matplotlib.pyplot as plt
 2 import numpy as np
 3 
 4 def log_type(name,arr):
 5     print("数组{}的大小:{}".format(name,arr.size))
 6     print("数组{}的维度:{}".format(name,arr.shape))
 7     print("数组{}的维度:{}".format(name,arr.ndim))
 8     print("数组{}元素的数据类型:{}".format(name,arr.dtype))
 9     #print("数组:{}".format(arr.data))
10 
11 a = np.floor(10*np.random.random((1,3,4)))
12 print(a)
13 log_type(''a'',a)
14 
15 a1 = a.squeeze()
16 print("修改前a1:{}".format(a1))
17 log_type(''a1'',a1)
18 a1[2] = 100
19 print("修改后a1:{}".format(a1))
20 
21 print("a:{}".format(a))
22 log_type(''a'',a)

今天关于Curve_fit 用于返回 numpy 数组的函数和sumproduct函数返回value的讲解已经结束，谢谢您的阅读，如果想了解更多关于2D Gaussian Curve_fit 与有缺陷的数据、Curve_Fit 不准确、numpy.random.random & numpy.ndarray.astype & numpy.arange、numpy.ravel()/numpy.flatten()/numpy.squeeze()的相关知识，请在本站搜索。

本文标签：