【问题标题】:Python least squares with scipy.integrate.quad带有 scipy.integrate.quad 的 Python 最小二乘法
【发布时间】:2013-09-21 00:47:20
【问题描述】:

我无法为我的代码破译错误消息以找到适合两个参数(eps 和 sig)的复杂最小二乘的一些参数。

from pylab import *
import scipy
import numpy as np

from scipy import integrate, optimize

# Estimate parameters with least squares fit   
T = [90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250, 260, 270, 280, 290, 300]
B = [-0.2221, -0.18276, -0.15348, -0.13088, -0.11293, -0.09836, -0.086301, -0.076166, -0.067535, -0.060101, -0.053636, -0.047963, -0.04295, -0.038488, -0.034494, -0.030899, -0.027648, -0.02469, -0.022, -0.019534, -0.017268, -0.015181]


def funeval(Temp,eps,sig):
    return -2.*np.pi*scipy.integrate.quad( lambda x: np.exp(4.*eps/Temp*((sig/x)**6.-(sig/x)**12.)*(x**2)) ,0.0,Inf )[0]

def residuals(p,y,Temp):
    eps,sig = p
    err = y-(funeval(Temp,eps,sig) )
    return err

print funeval(90.,0.001, 0.0002) 

plsq = scipy.optimize.leastsq(residuals, [0.00001, 0.0002], args=(B, T))

funeval 给出了一个合理的浮点数,但是当我运行它返回的代码时:

error: Supplied function does not return a valid float.

错误似乎对初始条件不敏感。我是 python 新手,所以任何帮助或帮助指南将不胜感激。谢谢。

【问题讨论】:

    标签: python numpy scipy least-squares


    【解决方案1】:

    funeval(90.,0.001, 0.0002)Temp为奇异值;但是,当您调用scipy.optimize 时,您会将整个T 数组传递给funeval,从而导致scipy.integrate 崩溃。

    快速解决方法是:

    def funeval(Temp,eps,sig):
        out=[]
        for T in Temp:
            val = scipy.integrate.quad( lambda x: np.expm1( ((4.*eps)/T)* ((sig/x)**12.-(sig/x)**6.)* (x**2.) ), 0.0, np.inf )[0]
            out.append(val)
        return np.array(out)
    
    def residuals(p,y,Temp):
        eps,sig = p
        err = y-(funeval(Temp,eps,sig) )
        return err
    
    print funeval([90],0.001, 0.0002)
    
    plsq = scipy.optimize.leastsq(residuals, [0.00001, 0.0002], args=(B, T))
    (array([  3.52991175e-06,   9.04143361e-02]), 1)
    

    不幸的是,这并不能很好地收敛。你能解释一下你想做什么吗?

    【讨论】:

    • 好的,感谢您到目前为止的帮助。这个想法是在给定第二维里系数的数据的情况下计算 Argon 的 Lennard Jones 势参数。所以描述这种关系的方程是i.imgur.com/VfCl5EE.gif 所以这个想法是对epsilon 和sigma 进行微调。 (应该是120和341e-9左右)(还在学习SO界面)
    • @doubletheron 请查看我的编辑,您忘记了-1。你也应该使用np.expm1,请查看页面以获得解释。您可能需要考虑加权不同的温度。
    • 我发现教科书在撒谎。第 1 册(骗子):i.imgur.com/fAO4Fdl.jpg 第 2 册(可能说的是真话):i.imgur.com/GVJVSaT.jpg 所以我在funeval 中修复了 eqn,但仍然没有收敛。
    • B 的单位是什么?它通常是cm^3/mol,但不存在这种转换。这是作业题吗?
    • B 在m^3/kmol 中我在 REFPROP 中计算了 1 bar 的氩气。我希望这是一个硬件问题,然后我就会有信心跳过它;不,这只是我试图开始工作以试图弄清楚物理并在 python 上改进的东西。
    猜你喜欢
    • 2017-09-22
    • 2015-08-01
    • 1970-01-01
    • 1970-01-01
    • 2014-04-28
    • 2012-02-10
    • 2017-03-18
    • 2013-05-20
    • 1970-01-01
    相关资源
    最近更新 更多