scipy.statsパッケージを使用してリストデータをガウシアンフィッティングし、ピーク値から1/e^2になる所の幅を知りたい。

ガウシアンの形を持つデータをリストに格納し、このデータをガウシアンフィットし、図に表現したいです。
更に、フィッティングしたガウシアンのピーク値から1/e^2に落ちる値の幅を知りたいです。
以下にデータをリストに格納しピーク値でリスト全体を割り、ピーク値を１にするまでのコードをのせます。

dat=np.genfromtxt('test.csv',delimiter=',',filling_values=(0))
tate = np.sum(dat, axis=0)
yoko = np.sum(dat, axis=1)
matate=max(tate)
mayoko=max(yoko)
tate_1=tate/matate
yoko_1=yoko/mayoko
from scipy.stats import norm
param=norm.fit(tate_1)

下に示すのがガウシアン形のデータです。

[114, 113, 112, 112, 114, 116, 114, 113, 112, 115, 112, 111, 109, 106, 105, 105, 105, 104, 106, 103, 106, 108, 105, 104, 102, 103, 106, 108, 108, 109, 111, 111, 115, 116, 119, 119, 122, 122, 125, 125, 128, 133, 135, 138, 142, 147, 149, 150, 155, 158, 161, 167, 168, 172, 171, 177, 179, 181, 184, 187, 192, 195, 200, 210, 214, 216, 222, 226, 233, 236, 244, 250, 257, 266, 273, 278, 287, 299, 302, 314, 321, 328, 333, 343, 357, 366, 378, 382, 390, 401, 410, 431, 434, 446, 443, 459, 461, 476, 488, 505, 511, 517, 530, 536, 549, 562, 578, 592, 599, 601, 612, 629, 637, 659, 671, 681, 694, 709, 727, 743, 764, 790, 802, 814, 831, 835, 844, 866, 892, 922, 947, 979, 994, 1009, 1025, 1051, 1066, 1089, 1113, 1153, 1183, 1197, 1240, 1269, 1298, 1316, 1353, 1378, 1406, 1417, 1444, 1465, 1495, 1514, 1549, 1597, 1631, 1656, 1688, 1715, 1740, 1755, 1778, 1803, 1838, 1853, 1872, 1902, 1924, 1927, 1968, 1986, 2027, 2051, 2085, 2116, 2151, 2185, 2209, 2246, 2262, 2297, 2337, 2389, 2434, 2468, 2503, 2560, 2582, 2626, 2659, 2688, 2716, 2778, 2819, 2864, 2926, 2954, 3000, 3023, 3075, 3120, 3153, 3192, 3240, 3299, 3358, 3415, 3446, 3494, 3534, 3581, 3640, 3677, 3704, 3726, 3770, 3815, 3871, 3932, 3992, 4026, 4055, 4074, 4089, 4100, 4131, 4166, 4215, 4274, 4333, 4408, 4472, 4533, 4575, 4621, 4661, 4681, 4725, 4770, 4826, 4871, 4927, 4980, 5070, 5170, 5231, 5254, 5274, 5294, 5334, 5377, 5442, 5493, 5568, 5613, 5649, 5689, 5685, 5740, 5813, 5885, 5937, 5993, 6051, 6110, 6130, 6145, 6209, 6292, 6381, 6407, 6465, 6502, 6526, 6586, 6681, 6759, 6801, 6871, 6924, 6944, 6999, 7074, 7141, 7216, 7280, 7321, 7356, 7412, 7444, 7489, 7508, 7568, 7622, 7694, 7781, 7819, 7862, 7861, 7937, 8006, 8051, 8073, 8122, 8182, 8257, 8285, 8317, 8335, 8352, 8382, 8423, 8465, 8488, 8542, 8573, 8595, 8612, 8669, 8718, 8736, 8791, 8852, 8914, 8942, 8993, 9048, 9084, 9116, 9139, 9170, 9225, 9245, 9272, 9274, 9322, 9394, 9429, 9461, 9516, 9543, 9537, 9589, 9617, 9635, 9649, 9680, 9700, 9718, 9726, 9774, 9819, 9802, 9823, 9801, 9812, 9804, 9851, 9875, 9874, 9891, 9951, 9943, 9915, 9916, 9920, 9945, 9947, 9958, 9937, 9953, 9983, 9958, 9931, 9947, 9988, 9986, 9996, 9978, 10000, 9996, 9963, 9954, 9952, 9960, 9932, 9892, 9889, 9849, 9835, 9837, 9834, 9837, 9795, 9772, 9754, 9748, 9728, 9723, 9675, 9601, 9580, 9558, 9534, 9531, 9500, 9495, 9493, 9468, 9460, 9413, 9386, 9378, 9315, 9286, 9236, 9224, 9217, 9209, 9198, 9162, 9129, 9086, 9050, 9029, 8993, 9004, 8985, 8955, 8865, 8827, 8775, 8763, 8753, 8721, 8703, 8667, 8650, 8601, 8558, 8515, 8466, 8411, 8365, 8314, 8272, 8234, 8213, 8169, 8134, 8099, 8056, 8006, 7932, 7867, 7814, 7772, 7725, 7679, 7637, 7566, 7512, 7499, 7471, 7416, 7343, 7294, 7231, 7171, 7121, 7070, 7005, 6949, 6898, 6829, 6774, 6729, 6665, 6587, 6523, 6492, 6431, 6352, 6309, 6267, 6184, 6145, 6119, 6073, 6033, 5986, 5929, 5899, 5844, 5811, 5749, 5704, 5645, 5591, 5534, 5454, 5403, 5345, 5305, 5258, 5228, 5179, 5144, 5093, 5051, 5010, 4959, 4909, 4861, 4815, 4793, 4739, 4699, 4643, 4603, 4546, 4502, 4459, 4401, 4380, 4330, 4286, 4243, 4204, 4171, 4141, 4102, 4053, 4015, 3970, 3930, 3891, 3830, 3795, 3749, 3710, 3656, 3617, 3581, 3532, 3494, 3475, 3433, 3403, 3348, 3315, 3279, 3237, 3218, 3160, 3110, 3068, 3036, 3003, 2998, 2960, 2913, 2849, 2803, 2750, 2732, 2704, 2674, 2660, 2613, 2581, 2536, 2491, 2452, 2419, 2381, 2364, 2329, 2294, 2253, 2215, 2190, 2167, 2123, 2082, 2033, 2002, 1963, 1938, 1908, 1870, 1830, 1801, 1768, 1737, 1710, 1672, 1644, 1610, 1589, 1568, 1547, 1510, 1481, 1443, 1404, 1388, 1360, 1329, 1319, 1294, 1271, 1251, 1233, 1209, 1180, 1156, 1138, 1123, 1108, 1088, 1075, 1058, 1047, 1021, 1004, 987, 972, 961, 945, 937, 920, 903, 879, 872, 859, 847, 832, 818, 798, 788, 768, 748, 732, 725, 702, 692, 670, 659, 644, 625, 612, 593, 578, 573, 568, 550, 536, 518, 504, 492, 486, 477, 466, 455, 447, 438, 431, 419, 413, 413, 402, 395, 393, 386, 379, 376, 374, 371, 367, 360, 361, 356, 355, 355, 359, 353, 351, 340, 343, 338, 341, 339, 340, 333, 335, 329, 329, 330, 332, 335, 334, 333, 331, 336, 332, 334, 334, 332, 326, 331, 334, 335, 337, 336, 338, 339, 341, 339, 341, 342, 347, 343, 349, 348, 346, 345, 347, 345, 340, 341, 341, 342, 342, 342, 335, 335, 337, 333, 333, 333, 335, 333, 328, 326]

よろしくお願いします。

行動規範の内容に同意します

回答1件

ベストアンサー

こんな感じでしょうか？

とりあえず

X軸をどのように取るのかわからなかったので(0～745) を振ってます
Fitting は scipy.optimize.curve_fit() を使用
パラメータとしては (A, μ, σ) の３つを導く
1 / e^2 の箇所は単に２シグマ区間で導く

としております

Python
1import numpy as np
2from scipy.stats import norm
3from scipy.optimize import curve_fit
4import matplotlib.pyplot as plt
5
6# データ読み込み
7dat=np.genfromtxt('test.csv',delimiter=',',filling_values=(0))
8
9N = len(dat)
10x = range(N)
11y = dat[:]
12
13def gauss(x, A, mu, sigma):
14    return A * norm(mu, sigma).pdf(x)
15
16# パラメータの初期値(A, mu, sigma)
17print(N)
18p0 = [1., 0., 1.]
19parameter,_ = curve_fit(gauss, x, y, p0=p0)
20print(parameter)
21
22# ２シグマ区間
23rmin = parameter[1] - 2 * parameter[2]
24rmax = parameter[1] + 2 * parameter[2]
25
26
27xx = np.linspace(0,N, 100)
28yy = gauss(xx, *parameter)
29
30ax = plt.subplot()
31ax.plot(x, y, 'c.')
32ax.plot(xx, yy, 'r-')
33ax.axvline(rmin, linestyle='--', color='blue', label='X = {}'.format(rmin))
34ax.axvline(rmax, linestyle='-.', color='blue', label='X = {}'.format(rmax))
35ax.legend(loc='upper right')
36plt.savefig('out.png')
37plt.show()