如何为任意函数定义卡方值函数?
我正在使用pyminuit这个Python库来进行数据拟合,它是minuit最小化代码的一个绑定(可以在http://code.google.com/p/pyminuit/找到)。这个最小化工具接受一个函数,并通过反射来提取需要最小化的参数。一般来说,我想要最小化的是给定数据集的卡方值,以便用一个特定的函数来描述这个数据集。
我想问的是:有没有办法定义一个卡方函数,这个函数可以接受一个参数数量不固定的任意函数,并返回一个只包含需要最小化的参数的函数,这个返回的函数可以计算出卡方值?
举个例子:
from scipy import *
import minuit
# Generate some data to fit
data_x = arange(50)
noise = 0.3
data_y = data_x**3 + normal(0.0, noise)
# Fit function, e.g. a cubic
fit_func = lambda x, a1, a2, a3, a4: a1 + a2*x + a3*x**2 + a4*x**3
# Minimisation function e.g. chi squared
# Note this has only the parameters to be minimised in the definition (eg not data_x)
min_func = lambda a1, a2, a3, a4: sum( (fit_func(data_x, a1, a2, a3, a4) - data_y)**2 / noise**2 )
在这里,我想写类似于min_func = make_chi2(fit_func)
的代码。但是我不知道该怎么做,因为data_x
和data_y
只在函数外部定义。为了完整性,最小化的其余部分看起来是这样的:
# Initialise minimiser object with initial values
m = minuit.Minuit(min_func, {'a1': 1.0, 'a2': 1.0, 'a3': 1.0, 'a4': 1.0})
# Run minimiser
m.migrad()
# Print minimised values - example output
print m.values
>>> {'a1': 0.000, 'a2': 0.000, 'a3': 0.000, 'a4': 1.000}
提前感谢你的帮助!
1 个回答
1
因为 PyMinuit 使用了自省(也就是程序可以检查自己的结构和状态),所以你也需要用自省的方式来处理。make_chi_squared()
可以这样实现:
import inspect
chi_squared_template = """
def chi_squared(%(params)s):
return (((f(data_x, %(params)s) - data_y) / errors) ** 2).sum()
"""
def make_chi_squared(f, data_x, data_y, errors):
params = ", ".join(inspect.getargspec(f).args[1:])
exec chi_squared_template % {"params": params}
return chi_squared
使用示例:
import numpy
def f(x, a1, a2, a3, a4):
return a1 + a2*x + a3*x**2 + a4*x**3
data_x = numpy.arange(50)
errors = numpy.random.randn(50) * 0.3
data_y = data_x**3 + errors
chi_squared = make_chi_squared(f, data_x, data_y, errors)
print inspect.getargspec(chi_squared).args
打印输出
['a1', 'a2', 'a3', 'a4']