5

我正在尝试为x,y点列表生成密度图，每个点都与给定的密度值相关联。看这张图片，看看我在追求什么。

我尝试在此答案中应用Joe Kington编写的代码，但它返回错误。numpy.linalg.linalg.LinAlgError: singular matrix

这是MWE我的一段代码（与 Joe 的代码基本相同，只是更改了数据数组）：

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate

x = np.array([0.005, 0.018, 0.008, 0.015, 0.016, 0.0135, 0.0155, 0.0155, 0.0105, 0.005, 0.0125, 0.0185, 0.0095, 0.003, 0.019, 0.0175, 0.0165, 0.011, 0.007, 0.0195, 0.017, 0.011, 0.0125, 0.0165, 0.0045, 0.0145, 0.02, 0.0185, 0.001, 0.015, 0.0105, 0.016, 0.0185, 0.0035, 0.0025, 0.0015, 0.0055, 0.0185, 0.005, 0.0135, 0.0175, 0.0095, 0.0095, 0.0115, 0.0025, 0.0105, 0.0015, 0.0045, 0.011, 0.009, 0.0045, 0.013, 0.016, 0.009, 0.018, 0.0145, 0.013, 0.0105, 0.019, 0.0145, 0.0145, 0.016, 0.014, 0.01, 0.018, 0.0075, 0.0195, 0.019, 0.012, 0.0035, 0.015, 0.0095, 0.0005, 0.0045, 0.011, 0.005, 0.019, 0.0015, 0.01, 0.0055, 0.0005, 0.018, 0.0155, 0.0065, 0.016, 0.002, 0.015, 0.0045, 0.0075, 0.0035, 0.0145, 0.018, 0.0145, 0.009, 0.0125, 0.005, 0.002, 0.0125, 0.013, 0.002, 0.007, 0.0125, 0.006, 0.015, 0.009, 0.0115, 0.0095, 0.016, 0.0045, 0.0035, 0.004, 0.0195, 0.0085, 0.0115, 0.011, 0.0175, 0.0115, 0.0085, 0.0185, 0.009, 0.007, 0.006, 0.0195, 0.018, 0.0115, 0.004, 0.017, 0.001, 0.0065, 0.018, 0.0175, 0.0115, 0.013, 0.0095, 0.0175, 0.0175, 0.012, 0.006, 0.015, 0.0075, 0.0015, 0.0085, 0.011, 0.0025, 0.001, 0.0105, 0.02, 0.0005, 0.0035, 0.0185, 0.0085, 0.011, 0.0005, 0.0055, 0.018, 0.019, 0.004, 0.01, 0.002, 0.0005, 0.0085, 0.0095, 0.0175, 0.0035, 0.0125, 0.0085, 0.0175, 0.011, 0.011, 0.0075, 0.0185, 0.0115, 0.0085, 0.005, 0.002, 0.003, 0.0095, 0.007, 0.011, 0.001, 0.0135, 0.003, 0.0125, 0.007, 0.0055, 0.0075, 0.019, 0.0055, 0.001, 0.0055, 0.003, 0.0085, 0.017, 0.01, 0.0065, 0.008, 0.013, 0.005, 0.0115, 0.005, 0.0055, 0.019, 0.001, 0.0095, 0.011, 0.008, 0.0165, 0.0195, 0.0025, 0.02, 0.0045, 0.0175, 0.018, 0.012, 0.0055, 0.008, 0.0025, 0.0155, 0.0055, 0.003, 0.0055, 0.0065, 0.011, 0.013, 0.0075, 0.0045, 0.005, 0.004, 0.0155, 0.0075, 0.0075, 0.0065, 0.0105, 0.0185, 0.0045, 0.0175, 0.0055, 0.0045, 0.0145, 0.015, 0.02, 0.007, 0.019, 0.019, 0.0075, 0.004, 0.009, 0.016, 0.0125, 0.0135, 0.012, 0.0145, 0.017, 0.006, 0.0085, 0.015, 0.0105, 0.0005, 0.003, 0.011, 0.0035, 0.0065, 0.011, 0.017, 0.003, 0.0145, 0.011, 0.0025, 0.0175, 0.011, 0.014, 0.01, 0.004, 0.0015, 0.0075, 0.0095, 0.009, 0.0195, 0.0025, 0.0135, 0.015, 0.006, 0.016, 0.016, 0.0105, 0.0065, 0.011, 0.019, 0.0145, 0.0065, 0.0185, 0.019, 0.0075, 0.0095, 0.0145, 0.013, 0.0175, 0.0085, 0.002, 0.005, 0.013, 0.0045, 0.0175, 0.01, 0.019, 0.004, 0.018, 0.0135, 0.0105, 0.0095, 0.0015, 0.0145, 0.019, 0.02, 0.012, 0.0145, 0.018, 0.0195, 0.001, 0.013, 0.0145, 0.0015, 0.0025, 0.0085, 0.0075, 0.0035, 0.019, 0.0125, 0.017, 0.0145, 0.003, 0.011, 0.0135, 0.003, 0.0085, 0.0195, 0.011, 0.015])
y = np.array([7.1, 9.35, 6.7, 9.9, 7.85, 7.4, 8.3, 8.8, 9.55, 8.55, 9.45, 8.7, 8.45, 8.0, 7.6, 8.45, 7.25, 7.95, 7.6, 7.75, 8.4, 7.55, 7.5, 7.1, 8.7, 7.0, 8.55, 9.45, 6.95, 9.6, 9.4, 7.65, 7.7, 9.15, 8.15, 9.95, 8.1, 8.6, 9.2, 8.8, 7.8, 7.85, 9.2, 6.9, 8.1, 9.8, 8.05, 7.45, 7.05, 7.3, 9.25, 8.8, 7.7, 9.25, 9.0, 6.7, 6.9, 9.15, 7.9, 7.35, 7.8, 7.35, 8.0, 7.0, 9.75, 8.95, 7.75, 9.25, 8.95, 7.6, 7.2, 7.6, 9.35, 7.3, 7.0, 6.65, 7.0, 9.9, 6.85, 9.8, 7.6, 9.35, 7.7, 8.2, 9.55, 7.0, 8.05, 9.95, 7.2, 9.7, 9.65, 8.85, 7.2, 9.9, 7.35, 6.9, 7.65, 9.8, 8.6, 8.75, 8.8, 7.1, 7.05, 7.8, 8.9, 8.15, 10.05, 9.95, 6.85, 7.1, 8.6, 9.45, 7.4, 7.45, 6.6, 9.9, 7.7, 9.95, 9.9, 8.1, 10.0, 7.3, 7.85, 9.95, 8.6, 8.55, 9.7, 8.6, 8.8, 9.6, 8.45, 6.65, 7.05, 8.5, 8.85, 9.55, 9.75, 8.95, 6.9, 6.8, 7.15, 9.95, 9.05, 8.1, 9.4, 8.05, 7.55, 8.85, 9.9, 9.65, 6.65, 8.6, 7.15, 8.95, 8.45, 8.8, 9.3, 9.4, 8.3, 9.85, 7.45, 6.85, 7.25, 7.55, 7.35, 9.5, 9.85, 9.75, 7.75, 7.55, 6.65, 6.6, 7.35, 8.25, 7.5, 8.2, 9.0, 7.95, 8.3, 8.15, 7.7, 7.45, 7.75, 7.8, 8.7, 9.9, 8.2, 7.1, 8.9, 8.85, 6.8, 7.2, 7.1, 7.65, 8.7, 6.9, 9.4, 7.25, 9.8, 8.4, 7.6, 8.5, 7.95, 6.7, 8.45, 9.2, 8.8, 7.85, 7.95, 8.7, 7.55, 9.6, 8.85, 8.9, 8.1, 8.25, 8.1, 8.3, 8.9, 7.1, 9.8, 8.25, 8.75, 6.85, 8.9, 6.95, 9.0, 8.35, 9.0, 7.15, 8.9, 8.2, 9.15, 6.65, 9.35, 8.85, 6.85, 7.8, 8.4, 7.75, 7.55, 7.85, 7.6, 8.2, 7.15, 8.55, 7.8, 8.8, 9.75, 9.0, 9.65, 7.15, 7.3, 7.1, 9.7, 7.75, 8.85, 9.75, 7.75, 7.1, 9.8, 9.95, 7.0, 9.0, 6.65, 7.55, 6.7, 7.65, 9.7, 7.15, 8.6, 8.55, 7.0, 9.4, 7.25, 9.0, 9.45, 8.2, 7.9, 8.95, 8.05, 8.9, 7.7, 7.35, 7.55, 9.75, 8.8, 7.35, 8.2, 8.7, 8.7, 8.2, 7.6, 8.4, 7.15, 8.8, 7.25, 7.4, 7.65, 9.2, 7.3, 7.05, 8.45, 7.0, 8.55, 8.2, 8.45, 7.4, 9.15, 8.45, 7.15, 8.75, 7.05, 7.5, 9.45, 8.85, 7.15, 7.85, 8.9, 8.8, 9.2, 8.1, 9.95, 7.55, 7.4, 9.65, 6.85, 8.85, 8.9, 7.0, 7.2, 8.6, 7.4, 8.55, 7.45, 8.15, 9.45, 6.85])
z = np.array([3021.0279029149683, 4975.3037400799076, 2166.9841077494534, 6289.9297927621046, 1769.826392967929, 1718.8244103972752, 1762.4826301548458, 2892.0488281847693, 6271.3213266755065, 2065.6752057097788, 3376.244940630062, 2082.5656535205321, 1823.5812514071088, 1973.1591378311682, 1797.3251073485019, 2612.0911561842113, 2249.1162757223706, 2398.027412668795, 2502.482089998005, 1819.1869918508887, 2377.09819196745, 1781.4988210953811, 1706.0247161815421, 1909.0435719934635, 1662.0553564384486, 2132.6588030625549, 2136.2280746624447, 5130.5751044254675, 1793.6247353368949, 6337.4932727294181, 4462.7694292877422, 2110.5308215132864, 1867.4707084026049, 2088.1839351230669, 2461.8645333625827, 5419.1889489642499, 2129.4105626134383, 2005.6244970468119, 3827.1395925866591, 2513.2456828786935, 2164.255723310996, 2206.3593513204733, 2790.194999425913, 1877.3040520904212, 1879.8192626127952, 5842.4922672912217, 2439.7109266628595, 1825.4377748685583, 1992.6098863537443, 1903.4690337855423, 4266.0357702742913, 1910.2633981988065, 1881.8280410083426, 3366.6103402944168, 2647.1141831688333, 1875.7829962178498, 1918.1350622011096, 2053.737955354547, 1828.9518511755123, 1884.073961574501, 1865.9550644370991, 1763.9694644756794, 1926.8518278729125, 2097.1403545248913, 6198.7268871463293, 1708.009992247851, 1818.0568126526525, 4452.1016078799939, 2047.2960066111493, 1963.8657274866737, 2047.4948739008223, 1768.3096547617561, 2447.5028998716234, 2482.7108962691, 2128.0018288469919, 1806.7301882321021, 1772.3246603000719, 5803.738518186321, 1798.572994122713, 6342.6198323644467, 2719.0407047465369, 5872.4448971310285, 1718.0287344031356, 1975.8553803398997, 4089.0418254578526, 2827.6461855175298, 2047.5731199798217, 6182.5306850708821, 1799.2216814306257, 5821.8657988768718, 5261.9920805401234, 1728.5326365076648, 2498.4478787599023, 5711.0176941257841, 2132.2157852666737, 2305.8342275164105, 1742.9260379707273, 4495.0605229412349, 2307.025813218379, 2267.0830186364001, 1850.0099188084214, 1788.4569797768027, 1799.9473765786229, 1701.3152693662535, 2714.7002916944525, 1832.2386812881207, 6033.4178712040912, 5795.5855254086073, 2354.2479143787036, 3136.1334489510668, 2053.8469973983347, 6322.4332297897745, 1981.3654435442081, 1819.4461046157055, 2257.6059180738293, 4709.7418781823208, 1961.9279449958558, 6243.2423074873122, 6175.1119528389308, 2177.0613286219436, 6187.6795632100693, 1830.8592800553179, 1811.1712359480662, 6106.7472822509098, 1952.8809811764972, 1870.965173064684, 6228.8248244690431, 2699.164454284873, 2295.7322910170833, 3807.5109993951892, 2188.8297091094282, 1813.4642741233388, 2413.5448089794727, 1989.4306848559088, 2201.1048218029914, 6295.4190187566846, 6243.2423074873122, 2733.0328883407497, 2158.2110618263923, 2050.959769289871, 2543.0427216636222, 6233.4918957305254, 1919.4414073039704, 1703.7336685448629, 3858.7713493494748, 1943.9159186025965, 2079.3275013804068, 2319.7081124813781, 5994.9649670064709, 6251.8593579480848, 2139.6291622278786, 2162.581672501492, 2841.7605426723185, 2048.5187718112393, 2193.2036553224912, 2251.9333157773744, 2629.2190178228229, 3690.6316021925595, 2637.030920619778, 5842.4922672912217, 1884.2560494399575, 1843.1577358091208, 2163.1556153631163, 2651.9755344260984, 2209.7526452687962, 3108.0919418402241, 4387.2789976920412, 4698.6052131550832, 1947.2382624692616, 1721.7122378504125, 1883.567174876209, 1791.5387081525541, 1808.3840823374492, 2027.1259150370804, 3008.0670290601133, 2783.8860603503094, 3548.6751272484839, 2119.2582438340964, 1842.8306356209096, 1811.8758696043149, 1993.8806420387857, 2230.114543957538, 2524.6101073467694, 2369.3233933185866, 2317.7688836929988, 6134.6748393553416, 1885.4706325047302, 3154.9375544439058, 2314.5555039102414, 1766.5362521033173, 2819.8071270945788, 2689.7822575815635, 2204.2481365143462, 1971.6766514902765, 2272.7389701563638, 1756.1110070448294, 2800.2434013912875, 2093.2498937672021, 6187.6795632100693, 1841.0113159699797, 1860.7680087210833, 1982.8860785887384, 2801.9132242057758, 2744.0105956647167, 2025.3324054638258, 2722.3349116785962, 2684.0320952587235, 1877.8708415793587, 1997.7146695460535, 2104.1892082880063, 2725.8442629042142, 3625.9437800615337, 2145.9331536509681, 1964.5714745836017, 1941.6269985047763, 1840.9640642894606, 1858.1655072711599, 2122.3339549617058, 2076.5760478221541, 2185.2013499877644, 6148.6869234687683, 1996.7536642289936, 2231.2474284717237, 2098.8536029973343, 1742.4203420876274, 2306.2812497850091, 2721.5057732393784, 1875.4024744482288, 3064.375291034848, 1829.6452820628654, 1794.8335411785936, 1915.9536359891322, 2502.9808788537543, 1768.7422023935612, 6233.4309130206575, 2399.5441046522633, 1771.295681640004, 2095.2163400865065, 2170.1556654381593, 2121.6059394283293, 2207.0814350651663, 1809.6000363096575, 1908.4235529291209, 1856.0817854726031, 1862.2542145669672, 1778.8206888476191, 1754.4167167356909, 2747.6536379759777, 4392.6756654129649, 3565.8650546652166, 5476.6856895960746, 2786.4985484893873, 1915.735638108521, 2041.4366888494503, 6187.6795632100693, 2341.2486995710688, 2103.793796652817, 5865.8139190498505, 1699.9828833421673, 1993.3785877622747, 5964.8913804791482, 6156.2099074387725, 2413.3395548571748, 2099.4689601305859, 2915.1159717066835, 1782.5651767836252, 1850.0197312477569, 2104.4166049482315, 6206.2249413381178, 2990.8556200784769, 2235.4837197118086, 1952.9282245926095, 2414.4380556536712, 3224.5742704799181, 1712.3537287234201, 3377.4874074163745, 3562.3822787444183, 1674.4781633929065, 1770.3781243649291, 2512.7270580207683, 2519.4737433197356, 2600.7028461763307, 1687.6149740077039, 1827.2126312245898, 1765.946329287706, 5842.4922672912217, 2133.4553630401206, 1987.9591394650151, 1818.3518579455942, 1999.0936725599697, 1988.5131267757788, 2059.2467809359555, 1741.9951505552397, 2018.1693071570285, 2441.0117858276699, 2031.5820629850468, 2021.9955740863736, 2121.6430818723975, 1789.2775504127987, 5638.9892789948635, 2326.4686662613717, 2959.3842095472223, 1984.1449776696347, 1774.9347978340554, 2327.371938051715, 1881.5330953605135, 2139.1151745789821, 1669.0467180215926, 2525.0241948857788, 1906.6525847893236, 2557.5130308053958, 1900.4062816966011, 2095.8573933740408, 1773.6179090742717, 6266.4557619192337, 1935.0743271817748, 2135.1532490883042, 1738.9633786795393, 2115.9327681390764, 2580.6955754655746, 2720.535179681423, 1860.1178987262649, 5333.5674851558724, 1778.2603949578265, 2081.653719365047, 6048.5456782611936, 1904.0520140723243, 2135.8697908711961, 2297.7416719290022, 2290.7393924215753, 3308.5077581082523, 1689.2883297381291, 1754.1733576638849, 2134.7146778225488, 1856.944696747086, 1820.2830994284482, 6183.9447504226655, 2285.18744670974])

# Set up a regular grid of interpolation points
xi, yi = np.linspace(x.min(), x.max(), 100), np.linspace(y.min(), y.max(), 100)
xi, yi = np.meshgrid(xi, yi)
# Interpolate
rbf = scipy.interpolate.Rbf(x, y, z, function='linear')
zi = rbf(xi, yi)

plt.imshow(zi, vmin=z.min(), vmax=z.max(), origin='lower',
           extent=[x.min(), x.max(), y.min(), y.max()])
plt.scatter(x, y, c=z)
plt.colorbar()
plt.show()

（对不起，混乱的x,y,z数组，但这就是我的原始代码从代码的另一部分获取它们的方式）

鉴于我的数据格式，有没有办法绕过该错误？这个答案表明 usingnumpy.linalg.lstsq可以作为一种解决方法，但我无法使其工作。

4

2 回答 2

2

问题是您有多个具有相同 (x,y) 值和不同 z 值的点（参见例如点 1 和 81）。一个简单但有点小技巧的解决方案就是在 (x,y) 点上添加噪声。

x = x + np.random.normal(scale=1e-8, size=x.shape)
y = y + np.random.normal(scale=1e-8, size=y.shape)

如果你这样做，你将解决关于不可逆矩阵的问题。

我注意到的另一个问题是最终表面看起来真的被拉伸了，因为 RBF 函数假设 x 和 y 坐标的缩放比例相似，但在你的情况下它们不是。将它们重新缩放到零和一之间可以解决此问题。您仍然可以使用它们的原始比例进行绘图。

rescale = lambda x: (x - x.min()) / (x.max() - x.min())
xs = rescale(x)
ys = rescale(y)

我建议在添加噪声之前重新缩放，以便两个轴上的噪声缩放相同。

最后，您应该添加aspect='auto'到您的imshow调用中，以便生成的绘图是一个合理的形状。通过这些更改，我认为您的代码应该做您想做的事情。

于 2014-01-31T22:52:06.687 回答
2

正如@hunse 正确指出的那样，您遇到了 aLinAlgError因为您尝试解决的拟合问题是超定的-本质上，因为您对某些 x,y 对有多个 z 值，所以不可能找到一个可以通过的精确解决方案通过所有 x,y,z 点。这就是调用linalg.solve失败的原因。

尽管您找不到确切的解决方案，但您仍然可以通过替换为以最小二乘的方式解决此问题，正如Muhammad Alkarouri 的回答所暗示的那样。linalg.solvelinalg.lstsq

scipy.interpolate.Rbf这是解决超定情况的修补版本示例：

class LSQ_Rbf(scipy.interpolate.Rbf):

    def __init__(self, *args, **kwargs):
        self.xi = asarray([asarray(a, dtype=float_).flatten()
                           for a in args[:-1]])
        self.N = self.xi.shape[-1]
        self.di = asarray(args[-1]).flatten()

        if not all([x.size == self.di.size for x in self.xi]):
            raise ValueError("All arrays must be equal length.")

        self.norm = kwargs.pop('norm', self._euclidean_norm)
        r = self._call_norm(self.xi, self.xi)
        self.epsilon = kwargs.pop('epsilon', None)
        if self.epsilon is None:
            self.epsilon = r.mean()
        self.smooth = kwargs.pop('smooth', 0.0)

        self.function = kwargs.pop('function', 'multiquadric')

        # attach anything left in kwargs to self
        #  for use by any user-callable function or
        #  to save on the object returned.
        for item, value in kwargs.items():
            setattr(self, item, value)

        self.A = self._init_function(r) - eye(self.N)*self.smooth
        # use linalg.lstsq rather than linalg.solve to deal with 
        # overdetermined cases
        self.nodes = np.linalg.lstsq(self.A, self.di)[0]

我刚刚__init__从原始类中复制粘贴了函数并修改了最后一行。

这是输出的样子（我还根据 hunse 的建议重新调整了输入 x,y 值）：

在此处输入图像描述

如果您添加的噪声规模选择得当，hunse 的答案可能会给出与我非常相似的结果，但是在我看来，明确地找到最小二乘解决方案更有意义，并且几乎肯定会在数值上更加稳定。

于 2014-02-02T00:18:39.180 回答