python-2.7 - 如何使用 pymc 查找初始状态的变量以将最终状态与观察结果匹配?
问题描述
我有一个可以用高斯函数表示的初始状态,比如说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
我查看了此链接和其他各种问题和答案:使用 pymc 将两个正态分布(直方图)与 MCMC 拟合? 并编写了相同的代码,但我不知道应该在哪里输入前三个值。
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
谢谢
解决方案
推荐阅读
- java - 将数组排序为子数组
- python - 删除重复行,但保留其数据的并集
- css - 使用 CSS 格式化标签
- javascript - 选项卡处于活动状态时隐藏和显示 div(选中)
- android - How to link android resource values together?
- python - 如何从命令行“上传和部署”到 Elastic Beanstalk?
- php - 如何为何时发送电子邮件创建全局规则
- android - 浮动标签 TextInputLayout insideTextInputEdittext 添加清除图标右侧角与 onclick 功能
- git - 如何在 git 中管理构建 ID?作为消息,分支,两者还是其他?
- android - 在 Android 上创建 JAXBContext 会引发 ClassCastException