问将curve_fit或多拟合限制为单调函数
EN

Stack Overflow用户

提问于 2018-07-09 03:35:53

回答 1查看 895关注 0票数 1

我试图创建一个函数，该函数接受一组观察到的和预期的数据点，确定最佳的校准函数，并将此校准应用于整个数据集(其中数据点是子集)。但是，我想确保由scipy.optimize.curve_fit或numpy.polyfit拟合的多项式函数是单调的(一阶导数不改变符号)。

下面是我的当前测试代码：

#! /usr/bin/env python
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import numpy as np
import math
import sys
import random
random.seed(28)

#############
# Test Data #
#############
measuredMaximaMZ = [100.0, 201.0, 222.0, 401.0, 500.0]
presentCalibrants = [100.5, 200.5, 300.5, 400.5, 500.5]

Fake_Raw_X = np.linspace(100,500,1000)
Fake_Raw_Y = random.sample(xrange(100000), 1000)

#############
# Functions #
#############
def powerLaw(x,a,b,c):
    penalty = 0
    if b > 2.:
        penalty = abs(b-1.)*10000
    if b < 0.:
        penalty = abs(2.-b)*10000
    return a*x**b + c + penalty

class powerLawCall:
    def __init__(self,x,a,b,c):
        self.a = a
        self.b = b
        self.c = c
        self.f = lambda x: a*x**b+c
    def __call__(self,x):
        return self.f(x)
    def describe(self):
        return str(self.a)+"*X^"+str(self.b)+"+"+str(self.c)

def performCalibration(measured, expected):
    RMS = sys.maxint
    func = None

    # Power Law
    z = curve_fit(powerLaw, measured, expected)
    RMS_buffer = []
    for index, i in enumerate(measured):
        RMS_buffer.append((powerLaw(i,*z[0])-expected[index])**2)
    RMS_buffer = np.mean(RMS_buffer)
    RMS_buffer = math.sqrt(RMS_buffer)
    if RMS_buffer < RMS:
        RMS = RMS_buffer
        func = powerLawCall(0,*z[0])

    # Polynomials between 1 and len(expected)
    for i in range(1,len(expected)):
        RMS_buffer = []
        z = np.polyfit(measured, expected, i)
        f = np.poly1d(z)
        for index, j in enumerate(measured):
            RMS_buffer.append((f(j) - expected[index])**2)
        RMS_buffer = np.mean(RMS_buffer)
        RMS_buffer = math.sqrt(RMS_buffer)
        if RMS_buffer < RMS:
            RMS = RMS_buffer
            func = f
    return func

############# 
# Test Main #
#############
f = performCalibration(measuredMaximaMZ, presentCalibrants)

if isinstance(f,powerLawCall):
    label = f.describe()
else:
    label = ""
    for index,i in enumerate(f):
        if index < len(f):
            label += "{0:.2e}".format(i)+"x^"+str(len(f)-index)+" + "
        else:
            label += "{0:.2e}".format(i)

# Write results
with open("UPC Results.txt",'w') as fw:
    fw.write("Function: "+str(label)+"\n")
    fw.write("Expected\tOriginal\tCalibrated\n")
    for index, i in enumerate(measuredMaximaMZ):
        fw.write(str(presentCalibrants[index])+"\t"+str(i)+"\t"+str(f(i))+"\n")

# Plotting
X_data = []
Y_data = []
for index, i in enumerate(measuredMaximaMZ):
    X_data.append(presentCalibrants[index])
    Y_data.append(f(i))
newX = np.linspace(X_data[0],X_data[-1],1000)
newY = f(newX)

fig =  plt.figure(figsize=(8,6))
ax = fig.add_subplot(211)
plt.scatter(X_data,measuredMaximaMZ,c='b',label='Raw',marker='s',alpha=0.5)
plt.scatter(X_data,Y_data,c='r',label='Calibrated',marker='s',alpha=0.5)
plt.plot(newX,newY,label="Fit, Function: "+str(label))
plt.legend(loc='best')
plt.title("UPC Test")
plt.xlabel("Expected X")
plt.ylabel("Observed X")

ax = fig.add_subplot(212)
plt.plot(Fake_Raw_X,Fake_Raw_Y,c='b')
plt.plot(f(Fake_Raw_X),Fake_Raw_Y,c='r')
plt.xlabel("X")
plt.ylabel("Y")
plt.show()