使用数组参数向量化函数

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]]

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

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 函数返回一个 (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)

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

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) 网格。