Как использовать Pymc, чтобы найти переменные в начальном состоянии, чтобы сопоставить конечное состояние с наблюдением? - PullRequest
Новая
       26

Как использовать Pymc, чтобы найти переменные в начальном состоянии, чтобы сопоставить конечное состояние с наблюдением?

0 голосов
/ 05 июля 2019

У меня есть начальное состояние, которое может быть представлено функцией Гаусса, скажем, F(x), с тремя свободными параметрами "amplitude, centroid, sigma". Я должен суммировать F(x) с функцией преобразования, G(x)= exp(a*x)+b, с двумя свободными параметрами a и b.

Я хочу использовать pymc, чтобы найти эти пять свободных параметров, чтобы конечное состояние F(x)+G(x) представляло гауссову функцию с: 

amplitude=3.05
centroid=5.45
sigma=5.47

Я посмотрел на эту ссылку и другие вопросы и ответы: Установите два нормальных распределения (гистограммы) с MCMC, используя pymc? и написал тот же код, но я не знаю, куда мне вводить предыдущие три значения. 

import numpy as np
import matplotlib.pyplot as pl
from scipy.optimize import curve_fit
import pymc as mc



def GaussFunc(x, A, mu, sigma):
    return A * np.exp(-0.5 * ((x - mu) / sigma)**2)

def trFunc(x,a,b):
    return np.exp(a*x)+b


interval=np.arange(0, 10, 0.05)
A_ini=2.0
mu_ini=3.0
sigma_ini=1.0
initial=GaussFunc(interval, A_ini, mu_ini, sigma_ini, )



intervalT=np.arange(0, np.pi, np.pi/200.0)
a_t=0.2
b_t=-2.0
transf= trFunc(intervalT,a_t,b_t,)

final=np.zeros(200)
final=initial+transf

est_centroid_one = mc.Uniform("est_centroid_one", 0, 10 )
est_sigma_one = mc.Uniform( "est_sigma_one", 0, 5 )
est_height_one = mc.Uniform( "est_height_one", 0, 5 )

est_a_two = mc.Uniform("est_a_two", 0, 1 )
est_b_two = mc.Uniform("est_b_two", -3, 3 )

precision= 1./mc.Uniform("std", 0, 1)**2



@mc.deterministic( trace = False) 
def est_profile_1(x = interval, mu = est_centroid_one, sigma = est_sigma_one, A= est_height_one):
    return GaussFunc( x, A, mu, sigma )


@mc.deterministic( trace = False) 
def est_profile_2(x = intervalT, a = est_a_two, b = est_b_two):
    return trFunc( x, a, b )

@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 = final, observed = True)

model = mc.Model([est_centroid_one, 
                est_height_one,
                est_sigma_one,
                est_a_two,
                est_b_two,
                precision])

map_ = mc.MAP( model )
map_.fit()

mcmc = mc.MCMC( model )
mcmc.sample( 50000,40000 )


print est_centroid_one,est_height_one,est_sigma_one,est_a_two,est_b_two

Спасибо

...