CODEWITHSUNDEEP
pythonfunctionnumpy
Homer Jay Simpson
Homer Jay Simpson

Reputation: 1232

invalid value encountered in log in Python

I am trying to impement the random walk metropolis hastings algorithm which my code is :

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)

But when i run the created function i receive a warning that says

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)))

i checked the parenthesis but it seem to me ok.Probably some 0 (zero) creates the warning ?

what is wrong with my code?

I would appreciate any help from anybody.

Upvotes: 0

Views: 1229

Answers (1)

Nameless
Nameless

Reputation: 61

When you call the function with this variables. The result of 1+omega*(data-mucan)**2 has a lot of negative number and when calculate np.log(1+omega*(data-mucan)**2) code calculate np.log(negative number) and this is invalid value encountered in log. negative number for log is out of the range

Upvotes: 1

Related Questions