2

在 Python 中，我试图规范化两个数组，然后取它们重叠的区域的平均值来创建一个新的复合数组。

为此，我认为我必须：

  1. 找到重叠区域，
  2. 插值重叠的 y 值，
  3. 遍历以找到最佳拟合的归一化常数，然后
  4. 将碎片粘贴在一起以形成我的新曲线

使用一些半随机值，如下所示：

在此处输入图像描述

此代码适用于 y 值相距不太远的小型数据集，但当 Y1 和 Y2 之间存在数量级时 Python 会崩溃（显然是由于迭代）。这是代码：

X1o = [x for x in X1 if x > X2[0]]
X2o = [x for x in X2 if x < X1[-1]]
Y1o = [y for y in Y1[(len(Y1)-len(X1o)):]]
Y2o = [y for y in Y2[:len(X2o)]]
Y2o = list(interp(X1o,X2o,Y2o))

c = abs(min(Y1o)-max(Y2o))
Y2test = [y2+c for y2 in Y2o]
Y2s = []
d = 0.01*min(Y2test)
while min(Y2test) < max(Y1o):
  Y2test = [y+d for y in Y2test]
  Y2s.append(Y2test)
  plot(X1o,Y2test,c='k',alpha=0.5)

idx = min(map(lambda i: (u.squaredError(Y1o, i), i, Y2s.index(i)), Y2s))[-1]         
Yavg = [(y1+y2)/2 for y1,y2 in zip(Y1o,Y2s[idx])]
diff = Y2s[idx][0]-Y2o[0]

X = [x for x in X1 if x < X2[0]] + X1o + [x for x in X2 if x > X1[-1]]
Y = [y for x,y in zip(X1,Y1) if x < X2[0]] + Yavg + [y+diff for x,y in zip(X2,Y2) if x > X1[-1]]

我真的需要用具有数千个数据点的恒星光谱来做到这一点，并且在 y 值之间的分布高达 20 个数量级。

任何建议将不胜感激！

4

1 回答 1

1

您的代码将从 numpy 和使用较少的 python 列表中受益匪浅，这些列表效率低下，尤其是您的 line Y2s.append(Y2test)。当你的while周期太长时，你只会追加一个很长的列表，这既慢又低效。

话虽如此，代码的瓶颈是最小化。您目前正在使用 python 列表进行暴力破解。您将从使用scipy.optimize函数之一中受益匪浅。

以下是我会做的一些广泛的建议：

  1. 找到两个光谱的 x 坐标极值，将两者插入到公共 x 值的网格中。
  2. 使用 scipy.optimize.fmin 为您进行最小化并计算最佳标准化。
  3. 将归一化光谱的部分内插到公共网格

这是 fmin 的一些示例代码（未经测试）：

import numpy as np
import scipy.optimize as opt

# y1 = interpolated values for one of the spectra
# y2 = interpolated values for the other spectra, normalise this one 

def errfunc(p, a1, a2):
    return np.sum(a1 - a2 * p)

p0 = 1.  # initial guess
norm_factor = opt.fmin(errfunc, p0, args=(y1, y2))

这应该给你最合适的norm_factor。

于 2012-12-12T22:00:34.773 回答