Python中的日志中遇到的无效值

Question

我正在尝试实现我的代码是随机游走大都会黑斯廷斯算法：

import numpy as np
def rwmetrop(data,mu0=0,kappa0=1,alpha0=1,lambda0=1,nburn=1000,ndraw=10000,vmu=1,vomega=1):
    
    
    n = len(data)
    alpha1 = (n/2) + alpha0 - 1


    stdvmu = np.sqrt(vmu)
    stdvomega = np.sqrt(vomega)


    mu = np.random.normal(mu0,
                          np.sqrt(1/kappa0),
                          size =1)

    omega = np.random.gamma(shape = alpha0,
                                scale = 1/alpha0,
                                size  = 1)
    
    draws  =  np.empty([ndraw,2])

    acceptmu = 0
    acceptomega = 0


    it = -nburn

    while it < ndraw:
        
        
        mucan = np.random.normal(mu,stdvmu,1)
        logp = (kappa0/2)*((mu-mu0)**2-((mucan-mu0)**2)) + np.sum( np.log( 1+omega*(  (data-mu)**2 )  - np.log(1+omega*(data-mucan)**2))) 


        u  = np.random.uniform(0,1,1)

        if np.log(u) < logp:
            acceptmu = acceptmu + 1
            mu  = mucan


            omegacan = np.random.normal(omega,stdvomega,1)

            if omegacan > 0:
                logp = alpha1*(np.log(omegacan)-np.log(omega)) + lambda0*(omega-omegacan)+ sum(np.log(1+omega*( ((data-mu)**2)))- np.log(1+omegacan*( (data-mu)**2)) )
                u = np.random.uniform(0,1,1)

        if np.log(u) < logp:
            acceptomega = acceptomega + 1
            omega = omegacan


  
        if it>0:
            draws[it,0] = mu
            draws[it,1] = omega


        it = it+1
  
        
    return(draws)

from scipy.stats import cauchy
data = cauchy.rvs(loc=1, scale=1, size=100)

但是当我运行创建的函数时，我会收到一条警告

draws4=rwmetrop(data=data,mu0=0,kappa0=1,alpha0=1,lambda0=1,nburn=1000,ndraw=10000,vmu=1,vomega=1)
mu=draws4[:,0];mu

 RuntimeWarning: invalid value encountered in log
  logp = (kappa0/2)*((mu-mu0)**2-((mucan-mu0)**2)) + np.sum( np.log( 1+omega*(  (data-mu)**2 )  - np.log(1+omega*(data-mucan)**2)))

我检查了括号，但在我看来没问题。可能一些 0（零）会产生警告？

我的代码有什么问题？

感谢任何人的帮助。

Answer 1

当你用这个变量调用函数时。 1+omega*(data-mucan)**2 的结果有很多负数，当计算 np.log(1+omega*(data-mucan)**2) 代码计算 np.log（负数），这是 invalid value encountered in log。 log 的负数超出范围