如何舍入模型中的PyMC3先验参数并与之进行比较?

2024-06-26 14:00:51 发布

您现在位置:Python中文网/ 问答频道 /正文

如果事情看起来比需要的复杂,我提前道歉。试图使它比这更简单,但后来一些问题更难制定。这样就尽可能接近我的处境。我有一些数据,我想用一个类似于下面描述的函数来拟合。但是在实现pymc3fit时,就像在第二个代码部分中一样,我遇到了与这两个问题相关的错误。你知道吗

# some constants and fixed things
rin = 1.
rout = 1000.
tmin = 10.
x = np.arange(0.1, 10, 0.2)

# parameters to be fitted
a = 4.0
b = 0.5
c = 10

# define r_struct, depends on b, but is one value
r_struct = 0.1**(1./b)

# nrad below, the number of radial points,  
# depends on r_struct and the constant rout
# QUESTION 1 here
nrad = round(r_struct, -1) + round((rout - r_struct)/10., -1) * 2

# define the radial bins, logarithmically spaced
rbins = np.logspace(np.log10(rin), np.log10(rout), num=nrad+1, endpoint=True)
rad = 10**(( np.log10(rbins[1:]) + np.log10(rbins[:-1]))/2.)
t = rad**(-b)

# QUESTION 2 here
# t cannot go below tmin, so I do this
t[t < tmin] = tmin

# then I want to know where the limit of r_struct is here
iout = rad > r_struct
iin = rad <= r_struct

# depending on if rad is past r_struct or not
# I do two different things.
c1  = c_struct * (rad[iin]/50.)**(-1.*a)
taper = np.exp( -(rad[iout]/rad[iout][0])**(2.-a) )
c2 = c_taper * taper/taper[0]
tau = np.append(c1, c2) 

y_true = J0( x[:, np.newaxis] * rad) * (t * (1.0 - np.exp(-1.*tau) ) )
y_true = y_true.sum(axis=1)

y_obs = y_true + np.random.normal(size=len(x))*50
y_error = np.ones_like(y_true)*0.1*y_true.max()

plt.errorbar(x,y_obs, yerr= y_error, color='r', marker='.', ls='None', label='Observed')
plt.plot(x,y_true, 'b', marker='None', ls='-', lw=1, label='True')
plt.xlabel('x'); plt.ylabel('y'); plt.legend()

enter image description here

好的,现在我想把这个放到PyMC3中。我从以下几点开始。你知道吗

with pm.Model() as basic_model:
    # parameters
    a = pm.Uniform("a", 0, 10)
    b = pm.Uniform("b", 0, 10)
    c = pm.Uniform("c", 0,100)

    r_struct = 0.1**(1./b)
    nrad = round(r_struct, -1) + round((rout - r_struct)/10., -1) * 2

    rbins = np.logspace(np.log10(rin), np.log10(rout), num=nrad+1, endpoint=True)
    rad = 10**(( np.log10(rbins[1:]) + np.log10(rbins[:-1]))/2.)
    t = rad**(-b)

    t[t < tmin] = tmin
    iout = rad > r_struct
    iin = rad <= r_struct

    c1  = c_struct * (rad[iin]/50.)**(-1.*a)
    taper = np.exp( -(rad[iout]/rad[iout][0])**(2.-a) )
    c2 = c_taper * taper/taper[0]

    tau = np.append(c1, c2) 
    y_true = J0( x[:, np.newaxis] * rad) * (t * (1.0 - np.exp(-1.*tau) ) )
    y_true = y_true.sum(axis=1)

    func=pm.Deterministic('gauss',y_true)

    y = pm.Normal('y', mu=func, tau=1.0/y_error**2, observed=y_obs)

    trace = pm.sample(5000)

但是,这在第一次“round()”调用时就已经失败了。参见下面的问题。你知道吗

问题1.“nrad”需要是一个整数,它取决于pymc3前面的“b”(通过r\u struct),运行pymc3时如何将其舍入到和整数? 作为参考,我得到以下错误:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-11-fcf897a21b89> in <module>()
     15     r_struct = 0.1**(1./b)
     16
---> 17     nrad = round(r_struct, -1) + round((rout - r_struct)/10., -1) * 2
     18     rbins = np.logspace(np.log10(rin), np.log10(rout), num=nrad+1, endpoint=True)

TypeError: type TensorVariable doesn't define __round__ method

问题2:如何与之前使用PyMC3创建的数组进行比较?代码在比较时失败。你知道吗

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-12-7f77e70d147a> in <module>()
     26     t = rad**(-b)
     27
---> 28     t[t < tmin] = tmin
     29     iout = rad > r_struct
     30     iin = rad <= r_struct

任何关于如何澄清问题的建议,包括题目都欢迎!你知道吗


Tags: truenpstructtauradroundpmiin