向量化具有数组参数的函数

Question

我目前正在尝试对某些数据应用卡方分析。我想根据模型的两个系数绘制不同值的颜色图

def f(x, coeff):
    return coeff[0] + numpy.exp(coeff[1] * x)

def chi_squared(coeff, x, y, y_err):
    return numpy.sum(((y - f(x, coeff) / y_err)**2)

us = numpy.linspace(u0, u1, n)
vs = numpy.linspace(v0, v1, n)
rs = numpy.meshgrid(us, vs)

chi = numpy.vectorize(chi_squared)
chi(rs, x, y, y_error)

我尝试对函数进行矢量化，以便能够传递不同系数的网格来生成颜色图。

x,y的值y_err都是长度为n的一维数组。而u、v是各种变化系数。

但是这不起作用，导致 IndexError: invalid index to scalar variable.

这是因为 coeff 是作为标量而不是向量传递的，但是我不知道如何更正它。

更新

我的目标是获取坐标数组

rs = [[[u0, v0], [u1, v0],..,[un, v0]],...,[[u0, vm],..,[un,vm]]

其中每个坐标都是传递给卡方法的系数参数。这应该 return 一个二维数组，其中填充了适当坐标的卡方值

chi = [[c00, c10, ..., cn0], ..., [c0m, c1m, ..., cnm]]

然后我可以使用此数据使用 imshow 绘制颜色图

Answer 1

这是我第一次尝试运行您的代码：

In [44]: def f(x, coeff):
    ...:     return coeff[0] + numpy.exp(coeff[1] * x)
    ...: 
    ...: def chi_squared(coeff, x, y, y_err):
    ...:     return numpy.sum((y - f(x, coeff) / y_err)**2)

（我不得不删除最后一行的 (）

首先猜测可能的数组值：

In [45]: x = np.arange(3)
In [46]: y = x
In [47]: y_err = x
In [48]: us = np.linspace(0,1,3)
In [49]: rs = np.meshgrid(us,us)
In [50]: rs
Out[50]: 
[array([[ 0. ,  0.5,  1. ],
        [ 0. ,  0.5,  1. ],
        [ 0. ,  0.5,  1. ]]), 
 array([[ 0. ,  0. ,  0. ],
        [ 0.5,  0.5,  0.5],
        [ 1. ,  1. ,  1. ]])]

In [51]: chi_squared(rs, x, y, y_err)
/usr/local/bin/ipython3:5: RuntimeWarning: divide by zero encountered in true_divide
  import sys
Out[51]: inf

哎呀，y_err 不应该是 0。再试一次：

In [52]: y_err = np.array([1,1,1])
In [53]: chi_squared(rs, x, y, y_err)
Out[53]: 53.262865105526018

如果我将 rs 列表转换为数组，它也有效：

In [55]: np.array(rs).shape
Out[55]: (2, 3, 3)
In [56]: chi_squared(np.array(rs), x, y, y_err)
Out[56]: 53.262865105526018

现在，vectorize 的目的是什么？

f函数returns一个(n,n)数组：

In [57]: f(x, rs)
Out[57]: 
array([[ 1.        ,  1.5       ,  2.        ],
       [ 1.        ,  2.14872127,  3.71828183],
       [ 1.        ,  3.21828183,  8.3890561 ]])

让我们修改 chi_squared 给 sum 一个轴

In [61]: def chi_squared(coeff, x, y, y_err, axis=None):
    ...:     return numpy.sum((y - f(x, coeff) / y_err)**2, axis=axis)

In [62]: chi_squared(np.array(rs), x, y, y_err)
Out[62]: 53.262865105526018

In [63]: chi_squared(np.array(rs), x, y, y_err, axis=0)
Out[63]: array([  3.        ,   6.49033483,  43.77253028])

In [64]: chi_squared(np.array(rs), x, y, y_err, axis=1)
Out[64]: array([  1.25      ,   5.272053  ,  46.74081211])

我很想将 coeff 更改为 coeff0, coeff1，以便从一开始就对该参数的传递方式给予更多控制，但这可能没有什么不同。

更新

既然您已经更具体地说明了 coeff 值与 x、y 等的关系，我认为这可以通过简单的广播来解决。无需使用 np.vectorize.

首先，定义一个不同大小的网格；这样我们和代码就不会认为 coeff 网格的每个维度都与 x、y 值有任何关系。

In [134]: rs = np.meshgrid(np.linspace(0,1,4), np.linspace(0,1,5), indexing='ij')
In [135]: coeff=np.array(rs)
In [136]: coeff.shape
Out[136]: (2, 4, 5)

现在看看给定 coeff 和 x 时 f 的样子。

In [137]: f(x, coeff[...,None]).shape
Out[137]: (4, 5, 3)

coeff 实际上是 (4,5,1)，而 x 是 (1,1,3)，结果是 (4,5,3)（通过广播规则）

同样的事情发生在 chi_squared 内部，最后一步 sum 在最后一个轴上（尺寸 3）：

In [138]: chi_squared(coeff[...,None], x, y, y_err, axis=-1)
Out[138]: 
array([[  2.        ,   1.20406718,   1.93676807,   8.40646968,
         32.99441808],
       [  2.33333333,   2.15923164,   3.84810347,  11.80559574,
         38.73264336],
       [  3.33333333,   3.78106277,   6.42610554,  15.87138846,
         45.13753532],
       [  5.        ,   6.06956056,   9.67077427,  20.60384785,
         52.20909393]])
In [139]: _.shape
Out[139]: (4, 5)

每 coeff 对值一个值，(4,5) 网格。

向量化具有数组参数的函数

Vectorizing A Function With Array Parameter

python

numpy

vectorization

chi-squared

更新