Q&A
> フィッティングの初期パラメータを適当になるよう何度か試して変更しました、しかし綺麗なフィットになりませんでした。
https://hackmd.io/@tenoto/Skl_wArqD
の
「そこそこ良さそうな初期値を見つけたら、Chi-square で最適化してみます。」
と書かれてるところの下のグラフのように、データと(フィッティング後の)ガウス関数を重ね書きしたグラフを、質問に追加してください
リンク内容### 前提
スペクトルのガウスフィットを行う為に、Minuitライブラリを用いて最小二乗法でグラフをガウス関数でフィッティングを行い適当なパラメータを推定するプログラムを作成してます。
実現したいこと
現在、ライブラリ上ではエラーなく出力することは十分できたのですがスペクトルよりも鋭いピークでフィッティングというのには不適当な状態になってます。
Minuitのライブラリを使用してグラフのガウシアンフィットを適切に行うことが目的です。
コード
import numpy as np
import matplotlib as mtb
import matplotlib.pyplot as plt
import pandas as pd
from iminuit import Minuit
from probfit import BinnedChi2,Extended,gaussian,doublegaussian,Chi2Regression
df = pd.read_csv('Na1copy.txt',delim_whitespace=True,names=['channel','number'])
error = np.sqrt(df['number'])
fig,ax = plt.subplots(1,1,figsize=(13,7))
ax.plot(df['channel'] ,df['number'],'o')
ax.set_ylabel('count')
ax.set_xlabel('channel[ch]')
fig.savefig('out')
ext_gauss = Extended(gaussian)
#bc2 = BinnedChi2(ext_gauss,df['channel'],df['number'],error)
#これより最小二乗法でフィッティングする。
#データのx,y,yeを用いて線形のフィッティングの準備をする
x2r = Chi2Regression(ext_gauss,np.array(df['channel']),np.array(df['number']),np.array(error))
#初期値をm,cのようにしてx2rを解析する
fit = Minuit(x2r,mean=400,sigma=50,N=100,errordef=1.0)
#フィッティングを行う
fit.migrad()
#フィットして非対称的にエラーを評価する
fit.minos()
#結果を出力する
print(fit.params)
#最適化した値らを出力
print(fit.values)
#最適化したエラーを出力
print(fit.errors)
mean_fit = fit.values[0]
sigma_fit = fit.values[1]
fig = plt.figure(figsize=(10,5))
plt.subplot(111)
x2r.draw(fit)
fig.savefig('gauss_plot')
###出力
return Minos error 0 , 0
return Minos error 0 , 0
return Minos error 0 , 0
<ValueView of Minuit at 7f941ff3fb20>
mean: 400.0
sigma: 0.0021726906777423363
N: 82.26998476159888
<ErrorView of Minuit at 7f941ff3fb20>
mean: 0.00500187446840473
sigma: 0.00016493120008058114
N: 12.796241733594902
※スペクトルのsigmaは目分量で大体、100はありますが以上のような出力になっています。またグラフに適当な高さになっていません。
###使用したデータ
ファイル名、Na1copy.txt
300 10419
301 10401
302 10302
303 10321
304 10302
305 10427
306 10418
307 10477
308 10355
309 10373
310 10530
311 10264
312 10479
313 10384
314 10628
315 10686
316 10537
317 10710
318 10602
319 10734
320 10846
321 10789
322 10820
323 10993
324 11101
325 10940
326 11132
327 11358
328 11380
329 11395
330 11616
331 11599
332 11640
333 11759
334 11958
335 12066
336 12357
337 12462
338 12435
339 12721
340 12732
341 12878
342 13331
343 13471
344 13585
345 13644
346 13819
347 14322
348 14279
349 14421
350 14818
351 14787
352 15172
353 15255
354 15820
355 16006
356 15763
357 16286
358 16513
359 16658
360 17103
361 17332
362 17595
363 17812
364 17967
365 18402
366 18634
367 18914
368 18925
369 19432
370 19792
371 20251
372 20288
373 20185
374 20813
375 20818
376 21383
377 21762
378 21909
379 22354
380 22251
381 22257
382 23075
383 22870
384 22959
385 22988
386 23916
387 23831
388 23915
389 24090
390 24184
391 24372
392 24300
393 24635
394 24693
395 24945
396 24953
397 24473
398 24937
399 24990
400 25106
401 25239
402 24742
403 25287
404 24446
405 24702
406 24957
407 24519
408 24371
409 24368
410 24487
411 24068
412 23992
413 24073
414 23884
415 23418
416 23377
417 22849
418 22562
419 22727
420 22267
421 21637
422 21386
423 21520
424 20860
425 20657
426 20412
427 20179
428 19361
429 19441
430 19190
431 18843
432 18696
433 18016
434 17743
435 17562
436 16919
437 16672
438 16226
439 15981
440 15529
441 15289
442 15007
443 14400
444 14052
445 13431
446 13433
447 13141
448 12627
449 12250
試したこと
フィッティングの初期パラメータを適当になるよう何度か試して変更しました、しかし綺麗なフィットになりませんでした。
どうかよろしくお願いします。
グラフの画像を送付したいと思ってたのですが、方法が探したところ見つかりませんでした…。設定ってどうすればいいでしょうか…。
画像を送付するのできないなら、Twitterで画像を投稿してそのリンクを送ろうかななんて思ってますが…
形としては高さのあるガウシアンの両端の足の高さが非対称な元のスペクトルのピークのところに、
デルタ関数のように鋭いピークをもった関数でフィッティングされました。
> グラフの画像を送付したいと思ってたのですが、方法が探したところ見つかりませんでした…。
これ見てください
https://teratail.com/questions/31272
ボタンの配置やアイコンは、上記質問がされた当時と変わってるかもしれませんが、そのボタンは質問を作成・編集する画面のどこかにあります
丁寧にありがとうございます！
職場に行ってパソコンでその画像が見れるのが今週の土曜日なので、またその時に改めて送るので、よろしくお願いできますか？
回答、というところに画像を送付しました！よろしくお願いします。
画像は、回答ではなくて、質問を編集して追加してください
google colabで、質問のコードの先頭に
!pip install iminuit
!pip install probfit
を追加し、データを
df = pd.read_csv('Na1copy.txt',delim_whitespace=True,names=['channel','number'])
↓ 変更
from scipy.stats import norm
x_vals = np.linspace(start=300, stop=450, num=150)
density = norm.pdf(x=x_vals, loc=400.0, scale=25.0) * 1000000
df = pd.DataFrame(np.hstack([x_vals.reshape((-1, 1)), density.reshape((-1, 1))]), columns=['channel', 'number'])
とでっち上げて、他は質問のコードそのままで実行したら、結果は
<ValueView of Minuit at 3839b90>
mean: 399.9999352586955
sigma: 24.999972945071544
N: 1000000.3183422928
となり、でっち上げたデータに合ってました
表示されたグラフの赤線も、データに合ってました
上記の通り、似たようなデータをでっち上げたらうまくいきました
質問者さんが実際に使ってるデータではうまくいかない理由は、質問を編集して、他人が現象を再現できるようなデータを提示していただけたら、何か分かるかもしれません
(現象が再現できれば、本物のデータとは別のデータでもかまいません)
申し訳あるません。慣れてなくて簡単に画像を送付してしまいました。
ですがjbpb0様を勿論BAにするので、そこは安心してください！
さて、問題のコードですが、頂いた変更点を適用すると良いフィッティングが確かに確認できました！
いま、実際に使用したデータを質問のテキストに追加しましたが、どうでしょうか？
引き続きよろしくお願いします！
回答2件
0
ベストアンサー
###使用したデータ
ファイル名、Na1copy.txt
と同じテキストファイルを作成してgoogle colabにアップロードして、google colabで
python
1!pip install iminuit 2!pip install probfit
を実行してから、質問のコードをそのまま実行したら、結果は下記の通りで、
スペクトルよりも鋭いピーク
にはなりませんでした
┌───┬───────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐ │ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │ ├───┼───────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤ │ 0 │ mean │ 395.70 │ 0.07 │ -0.07 │ 0.07 │ │ │ │ │ 1 │ sigma │ 61.03 │ 0.10 │ -0.10 │ 0.10 │ │ │ │ │ 2 │ N │ 3.420e6 │ 0.004e6 │ -0.004e6 │ 0.004e6 │ │ │ │ └───┴───────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘ <ValueView of Minuit at 3ad0690> mean: 395.70307408983706 sigma: 61.03244776048152 N: 3420304.956424334 <ErrorView of Minuit at 3ad0690> mean: 0.0727042015252482 sigma: 0.09935932511199844 N: 4220.3636166626175
投稿2022/12/26 09:05
総合スコア7658
ありがとうございます！
一番フィッティングさせたいのは、全体的なスペクトルというより、ピークであるx=400の部分とその近傍のピークが半分になる程度までのxの領域に集中してフィッティングさせたいですが、このフィッティングの範囲をより限定的にする方法は分かりますでしょうか？
因みに私は、VSCodeで入力してAnacondaを使用してます。引き続きよろしくお願いします。
> フィッティングの範囲をより限定的にする方法は分かりますでしょうか？
は、この質問
> スペクトルよりも鋭いピークでフィッティングというのには不適当な状態になってます。
とは別内容なので、別の質問にしてください
わかりました！ではこの質問はここで閉じますが、また対応して頂けますか？
BAにしようとまた思いますので！
0
出力はこんな具合になっています。
投稿2022/12/24 01:26
総合スコア2
あなたの回答
tips
プレビュー
