c++ ヒストグラムにガウスフィッティングを行いたい

前提・実現したいこと

タイトル通りヒストグラムにガウスフィッティングを行いたいのですがうまくできません。

データは以下の通りです。ピークが8個できるので、まずは左側から一個ずつフィッティングしたいと考えています。
どなたかご教授１いただけると幸いです。

0.311 1
0.313 1
0.318 2
0.321 1
0.329 1
0.33 1
0.333 17
0.335 23
0.336 25
0.337 87
0.338 107
0.339 162
0.34 315
0.341 484
0.342 623
0.343 832
0.344 1140
0.345 1311
0.346 1465
0.347 1671
0.348 1787
0.349 1814
0.35 1376
0.351 1437
0.352 1416
0.353 1066
0.354 905
0.355 863
0.356 783
0.357 827
0.358 750
0.359 507
0.36 655
0.361 531
0.362 612
0.363 497
0.364 502
0.365 487
0.366 470
0.367 484
0.368 562
0.369 564
0.37 772
0.371 872
0.372 963
0.373 1277
0.374 872
0.375 2012
0.376 1724
0.377 1866
0.378 1665
0.379 1663
0.38 1809
0.381 1208
0.382 1309
0.383 911
0.384 946
0.385 847
0.386 656
0.387 698
0.388 645
0.389 520
0.39 558
0.391 466
0.392 455
0.393 486
0.394 364
0.395 408
0.396 412
0.397 345
0.398 419
0.399 131
0.4 494
0.401 382
0.402 290
0.403 332
0.404 273
0.405 279
0.406 323
0.407 334
0.408 358
0.409 443
0.41 497
0.411 575
0.412 796
0.413 916
0.414 1159
0.415 1352
0.416 1292
0.417 1877
0.418 1813
0.419 1727
0.42 1708
0.421 1857
0.422 1487
0.423 1379
0.424 1155
0.425 1086
0.426 948
0.427 905
0.428 679
0.429 635
0.43 560
0.431 568
0.432 551
0.433 487
0.434 446
0.435 447
0.436 415
0.437 363
0.438 348
0.439 323
0.44 391
0.441 364
0.442 339
0.443 364
0.444 313
0.445 353
0.446 291
0.447 332
0.448 353
0.449 301
0.45 509
0.451 514
0.452 681
0.453 822
0.454 1007
0.455 1392
0.456 1765
0.457 1827
0.458 1948
0.459 2131
0.46 2262
0.461 1960
0.462 1739
0.463 1328
0.464 1414
0.465 990
0.466 821
0.467 745
0.468 634
0.469 596
0.47 519
0.471 485
0.472 394
0.473 382
0.474 415
0.475 347
0.476 317
0.477 341
0.478 345
0.479 350
0.48 371
0.481 357
0.482 371
0.483 323
0.484 325
0.485 319
0.486 317
0.487 312
0.488 323
0.489 276
0.49 325
0.491 470
0.492 503
0.493 484
0.494 432
0.495 1194
0.496 1233
0.497 1826
0.498 89
0.5 6082
0.501 195
0.502 3433
0.503 1869
0.504 2757
0.505 896
0.506 1362
0.507 1025
0.508 1097
0.509 570
0.51 479
0.511 442
0.512 495
0.513 457
0.514 400
0.515 413
0.516 372
0.517 377
0.518 392
0.519 311
0.521 389
0.522 361
0.523 373
0.524 379
0.525 438
0.526 383
0.527 396
0.528 490
0.529 462
0.53 510
0.531 535
0.532 740
0.533 769
0.534 796
0.535 1068
0.536 1337
0.537 1558
0.538 1714
0.539 1854
0.54 2266
0.541 2014
0.542 2090
0.543 1854
0.544 1656
0.545 1311
0.546 1177
0.547 985
0.548 715
0.549 520
0.55 592
0.551 475
0.552 455
0.553 410
0.554 391
0.555 331
0.556 416
0.557 365
0.558 410
0.559 321
0.56 440
0.561 336
0.562 324
0.563 387
0.564 416
0.565 374
0.566 390
0.567 501
0.568 489
0.569 576
0.57 568
0.571 661
0.572 859
0.573 836
0.574 1131
0.575 1153
0.576 1409
0.577 1492
0.578 1929
0.579 1511
0.58 1994
0.581 2022
0.582 1915
0.583 1502
0.584 1457
0.585 1267
0.586 1048
0.587 895
0.588 746
0.589 603
0.59 648
0.591 490
0.592 457
0.593 508
0.594 433
0.595 506
0.596 491
0.597 502
0.598 607
0.599 31
0.6 830
0.601 659
0.602 602
0.603 719
0.604 690
0.605 918
0.606 1108
0.607 1129
0.608 1495
0.609 1600
0.61 1606
0.611 2107
0.612 2007
0.613 1984
0.614 2071
0.615 1659
0.616 1478
0.617 1319
0.618 974
0.619 912
0.62 722
0.621 446
0.622 334
0.623 163
0.624 94
0.625 158
0.626 53
0.627 41
0.628 16
0.629 10
0.63 8
0.631 8
0.632 3
0.633 4
0.634 2
0.635 3
0.876 478.021

発生している問題・エラーメッセージ

正規分布のようにフィッティングしたいですが、どうしても直線になってしまいます。

該当のソースコード

c++
1#include <iostream>
2#include <fstream>
3#include <string>
4#include <stdlib.h>
5#include <math.h>
6#include <vector>
7#include <iomanip>
8
9using namespace std;
10
11int main(){
12    
13    FILE *fp = popen("gnuplot", "w" );
14    
15    fputs("set terminal png\n",fp);
16    fputs("set out 'gaussian_x76.png'\n",fp);
17    fputs("set xrange [0.33:0.365] \n" , fp);
18    fputs("set xl \'xcoo\'\n",fp);
19    fputs("set yl \'Event Counts\'\n",fp);
20    fputs("f(x)=a*exp(-(x-b)**2/(2*c**2)) \n",fp);
21    fputs("a=1;b=1;c=1 \n",fp);
22    fputs("fit f(x) 'x76.txt' using 1:2 via a,b,c\n" ,fp);
23    fputs("plot 'x76.txt' with boxes fill solid lc rgb 'yellow',f(x) with line lc rgb 'black' \n",fp);
24    
25            
26            
27    fflush(fp);
28    cin.get();
29    pclose(fp);
30    
31    cout << "done" << endl;
32}
33

試したこと

違うデータならフィッティングができるので、コードというよりはデータが悪いのでしょうか

補足情報（FW/ツールのバージョンなど）

mac

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

回答2件

ベストアンサー

ガウス分布は単峰性の分布 (山のピークが1つ) なので、1つ目の画像のように多峰性 (山のピークが複数) の分布のモデルとするのには向いていません。この場合、

データの8個の峰が性質の異なるデータが混じり合った結果としてできた分布であるならば、層別でグループごとに分布を分けて8個の単峰型の分布に帰着してから、個別にガウス分布でフィッティングする

もしくは

8個のガウス分布とわかっているなら、混合ガウス分布を使って推定する

といった方法が考えられます。

以下 Python の scikit-learn を使った実装例を貼ります。
C++ のライブラリはわからないですが、探せばなにかしらあると思います。

Python での実装例

sklearn.mixture.GaussianMixture — scikit-learn 0.23.2 documentation

python
1import numpy as np
2from matplotlib import pyplot as plt
3from scipy.stats import norm
4from sklearn.mixture import GaussianMixture
5
6# データを読み込む。
7data, count = np.loadtxt("sample.csv", unpack=True)
8count = count.astype(int)
9
10# 集計前のデータに戻す。
11data = np.repeat(data, count)
12
13# 外れ値を削除する (なんらかの方法で外れ値を検出して削除できたと仮定)
14data = data[data <= data.max() * 0.9]
15
16# ガウス混合モデルで推定する。
17gmm = GaussianMixture(n_components=8, covariance_type="spherical")
18gmm.fit(data.reshape(-1, 1))
19
20menas = gmm.means_.ravel()  # 平均
21sigmas = np.sqrt(gmm.covariances_)  # 標準偏差
22weights = gmm.weights_
23
24fig, ax = plt.subplots(figsize=(10, 4))
25
26# ヒストグラムを描画する。
27ax.hist(data, bins=50, density=True, alpha=0.3, ec="k")
28
29# 推定した混合ガウス分布を描画する。
30x = np.linspace(data.min() - 0.05, data.max() + 0.05, 100)
31for w, mean, sigma in zip(weights, menas, sigmas):
32    y = norm.pdf(x, loc=mean, scale=sigma)
33    ax.plot(x, w * y, zorder=10)
34plt.show()

データに外れ値があるようなので、なんらかの方法で事前に削除してください (以下の赤丸部分)

混合ガウス分布による推定結果

投稿2020/09/06 17:30

編集2020/09/06 17:32

tiitoi

総合スコア21956

fana

2020/09/07 02:05

> 外れ値そのデータだけ整数になっていないみたいなので，何か要らない物な気がしますね

kai_programing

2020/09/07 03:26

多分必要なデータは0.3~0.7の範囲なので外れ値はすぐに削除できると思います。やり方は①そもそも範囲を決めて８個のヒストグラムに分けてフィッティングする②一気に混合ガウスでフィッティングするかのどちらかなのですね。①のほうは手間と時間がかかりそうなので②のほうで行いたいと思います。pythonで作っていただいたガウスフィッティングはまさに僕が作りたかった図です。自分でもやってみようと思います。丁寧に教えてくださりありがとうございました。

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

ヒューリスティックな話ですが……

(1)まず……

適当な方法(bin値の平均とか最大値とか何かに基づいて)で閾値を決め，

bin値がその閾値以下であるbinの値を全て0にしたデータを考えます．
例示のデータの場合だと，例えば，閾値=1000付近とか．
2. 1.によって，分布が孤立した複数の山に分かれることが期待できるので，山ごとに正規分布の当てはめを行います．
3. 2.のうち，当てはまりが悪いやつがある場合，そこには山が複数ある（足きりがうまく行っていない）と思われるので，そいつを再度Step1.に回す（再帰的にやる）

といった方法で，分布の個数とそれらの大まかなパラメータを推定します．
（各分布のパラメータは相応にCoarseでしょうが，分布の個数がうまく推定できればとりあえずＯＫな感じで．）

(2)次に，(1)の結果を出発点(初期値)として，大元のデータに対しての当てはめ（各分布のパラメータのrefinement）を行います．

投稿2020/09/07 02:15

fana

総合スコア11708