14

CCD のスペクトログラフで検出されたライン プロファイルを適合させようとしています。検討しやすいように、解決した場合に実際に解決したいものと非常によく似たデモンストレーションを含めました。

私はこれを見てきました： https://stats.stackexchange.com/questions/46626/fitting-model-for-two-normal-distributions-in-pymc と他のさまざまな質問と回答ですが、根本的に異なることをしていますやりたいことより。

import pymc as mc
import numpy as np
import pylab as pl
def GaussFunc(x, amplitude, centroid, sigma):
    return amplitude * np.exp(-0.5 * ((x - centroid) / sigma)**2)

wavelength = np.arange(5000, 5050, 0.02)

# Profile 1
centroid_one = 5025.0
sigma_one = 2.2
height_one = 0.8
profile1 = GaussFunc(wavelength, height_one, centroid_one, sigma_one, )

# Profile 2
centroid_two = 5027.0
sigma_two = 1.2
height_two = 0.5
profile2 = GaussFunc(wavelength, height_two, centroid_two, sigma_two, )

# Measured values
noise = np.random.normal(0.0, 0.02, len(wavelength))
combined = profile1 + profile2 + noise

# If you want to plot what this looks like
pl.plot(wavelength, combined, label="Measured")
pl.plot(wavelength, profile1, color='red', linestyle='dashed', label="1")
pl.plot(wavelength, profile2, color='green', linestyle='dashed', label="2")
pl.title("Feature One and Two")
pl.legend()

真値と測定値のプロット

私の質問: PyMC (またはいくつかのバリアント) は、上記で使用された 2 つのコンポーネントの平均、振幅、およびシグマを教えてくれますか?

実際の問題に実際に適合する関数はガウス関数ではないことに注意してください-したがって、「組み込み」の pymc.Normal() 型ではなく、汎用関数 (私の例では GaussFunc など) を使用して例を提供してください関数。

また、モデルの選択も別の問題であることは理解しています。したがって、現在のノイズでは、1 つのコンポーネント (プロファイル) だけが統計的に正当化される可能性があります。しかし、1、2、3 などのコンポーネントの最適なソリューションが何であるかを確認したいと思います。

また、PyMC を使用するという考えにも慣れていません。scikit-learn、astroML、またはその他のパッケージが完璧だと思われる場合は、お知らせください。

編集：

私は多くの方法で失敗しましたが、正しい軌道に乗っていたと思うことの 1 つは次のとおりです。

sigma_mc_one = mc.Uniform('sig', 0.01, 6.5)
height_mc_one = mc.Uniform('height', 0.1, 2.5)
centroid_mc_one = mc.Uniform('cen', 5015., 5040.)

しかし、機能する mc.model を構築できませんでした。

4

1 に答える 1

16

最も簡潔なPyMCコードではありませんが、読者を助けるためにその決定を下しました。これが実行され、（本当に）正確な結果が得られるはずです。

何をモデル化しているのか本当にわからないので、自由な範囲で均一な事前分布を使用することにしました。しかし、おそらく重心の位置についてのアイデアがあり、そこでより良い事前分布を使用できます。

### Suggested one runs the above code first.
### Unknowns we are interested in


est_centroid_one = mc.Uniform("est_centroid_one", 5000, 5050 )
est_centroid_two = mc.Uniform("est_centroid_two", 5000, 5050 )

est_sigma_one = mc.Uniform( "est_sigma_one", 0, 5 )
est_sigma_two = mc.Uniform( "est_sigma_two", 0, 5 )

est_height_one = mc.Uniform( "est_height_one", 0, 5 ) 
est_height_two = mc.Uniform( "est_height_two", 0, 5 ) 

#std deviation of the noise, converted to precision by tau = 1/sigma**2
precision= 1./mc.Uniform("std", 0, 1)**2

#Set up the model's relationships.

@mc.deterministic( trace = False) 
def est_profile_1(x = wavelength, centroid = est_centroid_one, sigma = est_sigma_one, height= est_height_one):
    return GaussFunc( x, height, centroid, sigma )


@mc.deterministic( trace = False) 
def est_profile_2(x = wavelength, centroid = est_centroid_two, sigma = est_sigma_two, height= est_height_two):
    return GaussFunc( x, height, centroid, sigma )


@mc.deterministic( trace = False )
def mean( profile_1 = est_profile_1, profile_2 = est_profile_2 ):
    return profile_1 + profile_2


observations = mc.Normal("obs", mean, precision, value = combined, observed = True)


model = mc.Model([est_centroid_one, 
              est_centroid_two, 
                est_height_one,
                est_height_two,
                est_sigma_one,
                est_sigma_two,
                precision])

#always a good idea to MAP it prior to MCMC, so as to start with good initial values
map_ = mc.MAP( model )
map_.fit()

mcmc = mc.MCMC( model )
mcmc.sample( 50000,40000 ) #try running for longer if not happy with convergence.

重要

アルゴリズムはラベリングに依存しないため、結果はprofile1からのすべての特性で表示される可能性がprofile2あり、その逆も同様です。

于 2013-03-03T19:58:11.023 に答える