0

我正在尝试（在 python 中）将一系列任意数量的高斯函数（由仍在改进的简单算法确定）拟合到数据集。对于我当前的样本数据集，我有 174 个高斯函数。我有一个适合的程序，但它基本上是复杂的猜测和检查，并且消耗了所有 4GB 的可用内存。

有没有办法使用 scipy 或 numpy 中的东西来实现这一点？

这是我要使用的，其中波长 [] 是 x 坐标列表，而通量c [] 是 y 坐标列表：

#选择一个高斯
在范围（0,2）中重复：
    对于范围内的 f（0，len（质心））：
        #遍历所有其他高斯
        对于我在范围内（0，len（质心））：
            如果我！= f：
                #对于每个波长，
                对于 w 波长：
                    #将每个的值附加到一个列表中，称为其他
                    others.append(height[i]*math.exp(-(w-centroid[i])**2/(2*width[i]**2)))

    #优化当前高斯的质心
        上一页 = 质心 [f]
        最佳 = 质心 [f]
        #选择一个数量级
        对于范围内的 p (int(round(math.log10(centroid[i]))-3-repeat),int(round(math.log10(centroid[i])))-6-repeat,-1)：
            #选择那个数量级的值
            对于范围内的 m (-5,9)：
                #改变当前项的值
                质心[f] = prev + m * 10 **(p)
                #在所有波长上递增，列出新值
                变异 = 0
                残差 = 0
                测试 = []
                #在每个波长上增加并评估此变化是否使 R^2 更大
                对于范围内的k（0，len（波长））：
                    test.append(height[i]*math.exp(-(wavelength[k]-centroid[f])**2/(2*width[i]**2)))
                    残差 += (test[k]+others[k]-cflux[k])**2
                    变异性 += (test[k]+others[k]-avgcflux)**2
                rsquare = 1-（残差/变异）
                #检查这个新拟合的 R^2 值
                如果 rsquare > bestr：
                    bestr = rsquare
                    最佳 = 质心 [f]
        质心[f] = 最佳

    #优化当前高斯的高度
        上一页 = 高度[f]
        最佳 = 高度 [f]
        #选择一个数量级
        对于范围内的 p (int(round(math.log10(height[i]))-repeat),int(round(math.log10(height[i])))-3-repeat,-1)：
            #选择那个数量级的值
            对于范围内的 m (-5,9)：
                #改变当前项的值
                高度[f] = prev + m * 10 **(p)
                #在所有波长上递增，列出新值
                变异 = 0
                残差 = 0
                测试 = []
                #在每个波长上增加并评估此变化是否使 R^2 更大
                对于范围内的k（0，len（波长））：
                    test.append(height[f]*math.exp(-(wavelength[k]-centroid[i])**2/(2*width[i]**2)))
                    残差 += (test[k]+others[k]-cflux[k])**2
                    变异性 += (test[k]+others[k]-avgcflux)**2
                rsquare = 1-（残差/变异）
                #检查这个新拟合的 R^2 值
                如果 rsquare > bestr：
                    bestr = rsquare
                    最佳 = 高度 [f]
        高度[f] = 最佳

    #优化当前高斯的宽度
        上一页 = 宽度[f]
        最佳 = 宽度 [f]
        #选择一个数量级
        对于范围内的 p (int(round(math.log10(width[i]))-repeat),int(round(math.log10(width[i])))-3-repeat,-1)：
            #选择那个数量级的值
            对于范围内的 m (-5,9)：
                如果 prev + m * 10**(p) == 0：
                    m+=1
                #改变当前项的值
                宽度[f] = prev + m * 10 **(p)
                #在所有波长上递增，列出新值
                变异 = 0
                残差 = 0
                测试 = []
                #在每个波长上增加并评估此变化是否使 R^2 更大
                对于范围内的k（0，len（波长））：
                    test.append(height[i]*math.exp(-(wavelength[k]-centroid[i])**2/(2*width[f]**2)))
                    残差 += (test[k]+others[k]-cflux[k])**2
                    变异性 += (test[k]+others[k]-avgcflux)**2
                rsquare = 1-（残差/变异）
                #检查这个新拟合的 R^2 值
                如果 rsquare > bestr：
                    bestr = rsquare
                    最佳 = 宽度 [f]
        宽度[f] = 最佳
        计数 += 1
        #print '{} of {} peaks 优化，迭代 {} of {}'.format(f+1,len(centroid),repeat+1,2)
        完成=圆形（100 *（计数/（浮点数（len（质心））* 2）），2）
        print '{}% 完成'.format（完成）
    print 'New R^2 = {}'.format(bestr)
4

1 回答 1

2

是的，使用 scipy 可能会做得更好（更容易）。但首先，将您的代码重构为更小的函数；它只是让阅读和理解正在发生的事情变得更加容易。

至于内存消耗：你可能在某个地方过度扩展了一个列表（others是一个候选者：我从来没有看到它被清除（或初始化！），而它被填充在一个四重循环中）。那，或者你的数据就是那么大（在这种情况下，你真的应该使用 numpy 数组，只是为了加快速度）。我不知道，因为您在不知道大小的情况下引入了各种变量（有多大wavelengths？others得到多大？数据数组的所有初始化是什么以及在哪里？）

此外，拟合 174 个高斯函数有点疯狂。要么寻找另一种方法来确定你想从数据中得到什么，要么把事情分开。从wavelengths变量来看，您似乎正在尝试在高分辨率光谱中拟合线；也许隔离大部分线并分别拟合这些孤立的组会更好。如果它们都重叠，我怀疑任何正常的拟合技术都会对你有所帮助。

最后，也许像pandas这样的包可以提供帮助（例如，计算子包）。

也许最后一点，因为我看到代码中有很多可以改进的地方。在某些时候，codereview也可能有用。虽然现在我猜你的内存使用是最有问题的部分。

于 2013-04-19T09:28:38.437 回答