Reputation: 1232
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
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