使用 celerite2 (pymc3)

问题描述

我最近从 celerite 切换到 celerite2 来模拟恒星光曲线。我已经密切关注 celerite 2 tutorial，使其适应我的需要。我的模型由平坦均值、抖动项、颗粒化项和 SHO 项组成，用于模拟非径向振荡。

数据输入：

观察结果从 lightkurve.KeplerLightCurve 对象通过以下方式转换：

import numpy as np
import lightkurve as lk

import pymc3 as pm
import theano.tensor as tt
import theano

import celerite2
from celerite2.theano import terms

datalist = lk.search_lightcurvefile('KIC005002732',cadence='long')
data = datalist[2:3].download_all()
lc = data.stitch().remove_outliers().remove_nans()
lc = lc[0:2001] # cut lightcurve

# Convert to ppm
lc.flux = (lc.flux - 1) * 1e6
lc.flux_err *= 1e6

#Find numax prior:
pg = lc.to_periodogram(method='lombscargle',normalization='psd',minimum_frequency=50,maximum_frequency=300,oversample_factor = 3)
snr = pg.flatten()
seis = snr.to_seismology()
numax_guess = seis.estimate_numax(window_width = 1)

time =     np.ascontiguousarray(lc.time,dtype=np.float64)
flux =     np.ascontiguousarray(lc.flux,dtype=np.float64)
flux_err = np.ascontiguousarray(lc.flux_err,dtype=np.float64)

先验：

先验（prior_*_m 和 prior_*_sd）完全是 numpy float64 标量。

# MEAN FUNCTION
prior_mean_m  = np.float64(0.0)
prior_mean_sd = np.float64(np.std(flux))

# JITTER TERM
prior_jitter_m  = np.float64(10)
prior_jitter_sd = np.float64(5)

# GRANULATION TERM(S)
prior_w_gran_m  = np.float64([2*np.pi*numax_guess.value/2])
prior_w_gran_sd = np.float64(100.0)

prior_a_gran_m  = np.float64([np.var(flux)])
prior_a_gran_sd = np.float64(np.var(flux)/100*50)

prior_Q_gran_m =  np.float64(1/np.sqrt(2))
prior_Q_gran_sd = np.float64(4)

prior_S_gran_m =  prior_a_gran_m/(prior_w_gran_m * prior_Q_gran_m)
prior_S_gran_sd=  np.abs(prior_Q_gran_sd/(prior_w_gran_m*prior_Q_gran_m)
                        - (prior_a_gran_m *prior_w_gran_sd)/(prior_w_gran_m**2*prior_Q_gran_m)
                        - (prior_a_gran_m*prior_Q_gran_sd)/(prior_w_gran_m*prior_Q_gran_m**2))

# NONRADIAL OSCILLATION TERM(S)
prior_rho_osc_m  = np.float64(1/numax_guess.value)
prior_rho_osc_sd = prior_rho_osc_m/100*10

prior_tau_osc_m  = np.float64(2*4/(2*np.pi*numax_guess.value)) # deterministic result for Q = 4,numax = numax_guess
prior_tau_osc_sd = prior_tau_osc_m/100*50

prior_sig_osc_m =  np.float64(np.sqrt(np.var(flux)))
prior_sig_osc_sd = prior_sig_osc_m/100*50

模型构建：

我为大多数分布都包含了一个形状参数，因为我希望我的代码能够将向量作为输入处理，以防我需要多个粒度/振荡项。

with pm.Model() as gp_model:
# MEAN FUNCTION ----------------------------------------------------------------------------------------------------------------
    mean_function = pm.normal("mean",mu=prior_mean_m,sd=prior_mean_sd)

# JITTER TERM ------------------------------------------------------------------------------------------------------------------
    jitter = pm.Lognormal("jitter",mu=prior_jitter_m,sigma=prior_jitter_sd)
    
# GRANULATION TERM(S) ----------------------------------------------------------------------------------------------------------
    num_gran_terms = np.size(prior_w_gran_m) # used later on
    #Set pymc3 priors for parameters
    Q_gran = prior_Q_gran_m
    a_gran = pm.Lognormal("a_gran",mu = np.log(prior_a_gran_m),sigma = prior_a_gran_sd,shape = np.size(prior_a_gran_m)) #parameter a of SHO term. Used to link w0 and S0. 
    w_gran = pm.Lognormal("w_gran",mu = np.log(prior_w_gran_m),sigma = prior_w_gran_sd,shape = np.size(prior_w_gran_m))
    S_gran = pm.Deterministic("S_gran",a_gran/(w_gran * Q_gran)) #S0 = a/(w0*Q)

    # Build terms
    if np.size(prior_w_gran_m) is 1:
        gran_terms = terms.SHOTerm(S0=S_gran,w0=w_gran,Q=Q_gran)
    else:
        gran_terms = terms.SHOTerm(S0=S_gran[0],w0=w_gran[0],Q=Q_gran)
        for i in range(1,np.size(prior_a_gran_m)):
            gran_terms += terms.SHOTerm(S0=S_gran[i],w0=w_gran[i],Q=Q_gran)
    
# NONRADIAL OSCILLATION TERM(S) ------------------------------------------------------------------------------------------------    
    num_osc_terms = np.size(prior_rho_osc_m) # used later on
    #Set pymc3 priors for parameters
    tau_osc = pm.Lognormal("tau_osc",mu = np.log(prior_tau_osc_m),sigma = prior_tau_osc_sd,shape = np.size(prior_tau_osc_m))
    rho_osc = pm.Lognormal("rho_osc",mu = np.log(prior_rho_osc_m),sigma = prior_rho_osc_sd,shape = np.size(prior_rho_osc_m))
    sig_osc = pm.Lognormal("sig_osc",mu = np.log(prior_sig_osc_m),sigma = prior_sig_osc_sd,shape = np.size(prior_sig_osc_m))

    # Build terms
    if np.size(prior_rho_osc_m) is 1:
        osc_terms = terms.SHOTerm(rho=rho_osc,tau=tau_osc,sigma=sig_osc)
    else:
        osc_terms = terms.SHOTerm(rho=rho_osc[0],tau=tau_osc[0],sigma=sig_osc[0])
        for i in range(1,np.size(prior_rho_osc_m)):
            osc_terms += terms.SHOTerm(rho=rho_osc[i],tau=tau_osc[i],sigma=sig_osc[i])       

# BUILD GP ---------------------------------------------------------------------------------------------------------------------
    kernel = terms.TermSum(gran_terms,osc_terms)

    gp = celerite2.theano.GaussianProcess(kernel,mean = mean_function)
    gp.compute(time,diag=flux_err ** 2 + jitter,check_sorted = True) # Compute cholesky factorisation of covariance matrix
    gp.marginal("obs",observed=flux) # Add marginal likelihood to model (observations)
    print("Initial log likelihood: {0}".format(gp.log_likelihood(flux)))
    pm.Potential("loglike",gp.log_likelihood(flux - mean_function)) # Define Potential to be used for optimization

现在，当我运行此代码时，出现以下错误：

TypeError: The two branches should have identical types,but they are TensorType(float64,vector) and TensorType(float64,(True,True)) respectively. This error Could be raised if for example you provided a one element list on the `then` branch but a tensor on the `else` branch.
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-16-b80a96105cf6> in <module>
     16     # Build terms
     17     if np.size(prior_w_gran_m) is 1:
---> 18         gran_terms = terms.SHOTerm(S0=S_gran,Q=Q_gran)
     19     else:
     20         gran_terms = terms.SHOTerm(S0=S_gran[0],Q=Q_gran)
~\Anaconda3\lib\site-packages\celerite2\terms.py in wrapped(target,*args,**kwargs)
    604                             break
    605 
--> 606             return to_wrap(target,**kwargs)
    607 
    608         return wrapped
~\Anaconda3\lib\site-packages\celerite2\theano\terms.py in __init__(self,eps,**kwargs)
    444     def __init__(self,*,eps=1e-5,**kwargs):
    445         self.eps = tt.as_tensor_variable(eps).astype("float64")
--> 446         super().__init__(**kwargs)
    447 
    448     def overdamped(self):
~\Anaconda3\lib\site-packages\celerite2\theano\terms.py in __init__(self,dtype)
     27     def __init__(self,dtype="float64"):
     28         self.dtype = dtype
---> 29         self.coefficients = self.get_coefficients()
     30 
     31     def __add__(self,b):
~\Anaconda3\lib\site-packages\celerite2\theano\terms.py in get_coefficients(self)
    480     def get_coefficients(self):
    481         m = self.Q < 0.5
--> 482         return [
    483             ifelse(m,a,b)
    484             for a,b in zip(self.overdamped(),self.underdamped())
~\Anaconda3\lib\site-packages\celerite2\theano\terms.py in <listcomp>(.0)
    481         m = self.Q < 0.5
    482         return [
--> 483             ifelse(m,self.underdamped())
    485         ]
~\Anaconda3\lib\site-packages\theano\ifelse.py in ifelse(condition,then_branch,else_branch,name)
    367             if then_branch_elem.type != else_branch_elem.type:
    368                 # If the types still don't match,there is a problem.
--> 369                 raise TypeError(
    370                     'The two branches should have identical types,but '
    371                     'they are %s and %s respectively. This error Could be '
TypeError: The two branches should have identical types,True)) respectively. This error Could be raised if for example you provided a one element list on the `then` branch but a tensor on the `else` branch.

这对我来说毫无意义（老实说，就像大多数 theano 错误消息一样）。为什么分支之一是布尔值？

我的尝试

我已经尝试在代码中省略形状参数，因为我之前遇到过这个问题。

a_gran = pm.Lognormal("a_gran",sigma = prior_a_gran_sd)
...

这会导致以下错误：

TypeError: For compute_test_value,one input test value does not have the requested type.

The error when converting the test value to that variable type:
Wrong number of dimensions: expected 0,got 1 with shape (1,).
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-70-d0ae027fe60a> in <module>
     10     num_gran_terms = np.size(prior_w_gran_m)
     11     Q_gran = prior_Q_gran_m
---> 12     a_gran = pm.Lognormal("a_gran",sigma = prior_a_gran_sd) #parameter a of SHO term. Used to link w0 and S0.
     13     w_gran = pm.Lognormal("w_gran",sigma = prior_w_gran_sd)
     14     S_gran = pm.Deterministic("S_gran",a_gran/(w_gran * Q_gran)) #S0 = a/(w0*Q)
~\Anaconda3\lib\site-packages\pymc3\distributions\distribution.py in __new__(cls,name,**kwargs)
     45             total_size = kwargs.pop('total_size',None)
     46             dist = cls.dist(*args,**kwargs)
---> 47             return model.Var(name,dist,data,total_size)
     48         else:
     49             raise TypeError("Name needs to be a string but got: {}".format(name))
~\Anaconda3\lib\site-packages\pymc3\model.py in Var(self,total_size)
    924             else:
    925                 with self:
--> 926                     var = TransformedRV(name=name,distribution=dist,927                                         transform=dist.transform,928                                         total_size=total_size,~\Anaconda3\lib\site-packages\pymc3\model.py in __init__(self,type,owner,index,distribution,model,transform,total_size)
   1654 
   1655             self.transformed = model.Var(
-> 1656                 transformed_name,transform.apply(distribution),total_size=total_size)
   1657 
   1658             normalRV = transform.backward(self.transformed)
~\Anaconda3\lib\site-packages\pymc3\distributions\transforms.py in apply(self,dist)
    110     def apply(self,dist):
    111         # avoid circular import
--> 112         return Transformeddistribution.dist(dist,self)
    113 
    114     def __str__(self):
~\Anaconda3\lib\site-packages\pymc3\distributions\distribution.py in dist(cls,**kwargs)
     55     def dist(cls,**kwargs):
     56         dist = object.__new__(cls)
---> 57         dist.__init__(*args,**kwargs)
     58         return dist
     59 
~\Anaconda3\lib\site-packages\pymc3\distributions\transforms.py in __init__(self,**kwargs)
    139         self.dist = dist
    140         self.transform_used = transform
--> 141         v = forward(FreeRV(name="v",distribution=dist))
    142         self.type = v.type
    143 
~\Anaconda3\lib\site-packages\pymc3\model.py in __init__(self,total_size,model)
   1368             self.tag.test_value = np.ones(
   1369                 distribution.shape,distribution.dtype) * distribution.default()
-> 1370             self.logp_elemwiset = distribution.logp(self)
   1371             # The logp might need scaling in minibatches.
   1372             # This is done in `Factor`.
~\Anaconda3\lib\site-packages\pymc3\distributions\continuous.py in logp(self,value)
   1832         mu = self.mu
   1833         tau = self.tau
-> 1834         return bound(-0.5 * tau * (tt.log(value) - mu)**2
   1835                      + 0.5 * tt.log(tau / (2. * np.pi))
   1836                      - tt.log(value),~\Anaconda3\lib\site-packages\theano\gof\op.py in __call__(self,*inputs,**kwargs)
    623             for i,ins in enumerate(node.inputs):
    624                 try:
--> 625                     storage_map[ins] = [self._get_test_value(ins)]
    626                     compute_map[ins] = [True]
    627                 except AttributeError:
~\Anaconda3\lib\site-packages\theano\gof\op.py in _get_test_value(cls,v)
    560             # ensure that the test value is correct
    561             try:
--> 562                 ret = v.type.filter(v.tag.test_value)
    563             except Exception as e:
    564                 # Better error message.
~\Anaconda3\lib\site-packages\theano\tensor\type.py in filter(self,strict,allow_downcast)
    174 
    175         if self.ndim != data.ndim:
--> 176             raise TypeError("Wrong number of dimensions: expected %s,"
    177                             " got %s with shape %s." % (self.ndim,data.ndim,178                                                         data.shape))
TypeError: For compute_test_value,).

我也尝试过手动设置测试值，但会产生同样的错误。

a_gran = pm.Lognormal("a_gran",testval = np.log(prior_a_gran_m))
...

这是我第一次在这里发帖，如果我忘记了任何重要信息，请见谅。 预先感谢您的帮助！

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）