# Numerical Problems

Hi all,

I’m trying to write a simple code for an MCMC and I get into numerical errors when comparing two numbers in the metropolis step. The code for the metropolis step is
if rand() < exp( lokNew - logOld)
oldPar = NewPar
end
The problem is that my posterior distribution equals -1.0e10 when a given parameter is outside an specific range and in this case the draw should not be accepted. However, julia accepts the draw because when comparing
if rand()= < exp(-1.0e10 - logOld)
oldPar = NewPar
end
Where is my mistake?
I can share the code if that helps.
Andres

Do you mean `if rand() <= exp(-1.0e10 - logOld)`?

Yes

Is it not very strange? u is never equal to zero while the exp equals 0.0

What’s `logOld`? any chance it’s around 10^10?

logOld is a Float64 number close to -617.0

You don’t have enough context to reproduce this error `exp(-1.0e10 - logOld)=0.0`, so something else is the problem.

Hi Oscar, I can share the code. It is simple code. I’m sure it is my error but I cannot see it.
A.

using Plots
using StatsBase, StatsFuns, DataFrames
using Distributions
using NLsolve
using CSV
using LinearAlgebra
using Statistics

## Assumes balance panel

function DataTransPanel(X, N, T)
Xt = Array{Float64}(size(X)[1], size(X)[2])
K = size(X)[2]

``````for j=1:K
lb = 1;ub = 0;
for i=1:N
ub = lb + (T-1) ;
Xt[lb:ub, j] = X[lb:ub,j]-mean(X[lb:ub,j])
lb = ub + 1
end
end
return Xt
``````

end

function MeanPanel(X, N, T)
Xt = Array{Float64}(size(X)[1], size(X)[2])
K = size(X)[2]
for j=1:K # Number of variables in the VAR
lb = 1;ub = 0;
for i=1:N
ub = lb + (T-1) ;
Xt[lb:ub, j] = mean(X[lb:ub,j])
lb = ub + 1
end
end
# Returns a vector of means across time for each unit
return Xt
end

## Logistic variable

function LSTAR(q,theta, T)
T = convert(Int, T)
#tr = Array{Float64}(T, 1)
gam = exp(theta[1])
c = theta[2]

``````tr = (ones(T) +exp.((-gam.*(q-c*ones(T))))).^-1.0

#Returns a vector of means across time for each unit
return tr
``````

end

function PosteriorLSTAR(data, par)
T = size(data)[1] # Sample size
y = data[:,1];
x = data[:,2:3];
q = data[:, 4];

``````# Parameters

sigmaP = par[6];
beta = par[1:3];
theta = par[4:5];
tr =  LSTAR(q,theta, T);

xt = x[:,2].*tr;
xs = [x xt];

# prior for c
ra = quantile(q, [0.2, 0.8])

# #Prior for gamma
pgamma = log((1.0+exp(par[4])^2)^-1)

if (theta[2]<ra[1])
return(-1000.0-(theta[2]-ra[1])^2)
end
if (theta[2]>ra[2])
return(-1000.0-(theta[2]-ra[2])^2)
end
if (exp(theta[1])>100.0)
return(-1000.0-(theta[1]-100)^2)
end

sigmaP = exp(sigmaP);
res = y-xs*beta;

lik = -0.5*log(2*pi*sigmaP^2) - 0.5*dot(res,res)/(convert(Float64, T)*sigmaP^2) + pgamma;
return lik
``````

end

function mcmcLSTAR(mc, data, ini)
mc = convert(Int, mc)
T = size(data)[1]
StoreDraws = zeros(mc, 6) ; #Array{Float64}(mc, 6);

``````oldlik = PosteriorLSTAR(data, ini);

tryPar = ini;
StoreDraws[1,:] = ini;
d2 = Normal(0,0.1); # Candidate draw
d3 = Uniform(0,1);

for i in 2:mc
for j in 1:6
newbeta = tryPar;
newbeta[j] =  tryPar[j] + rand(d2, 1)[1]; # propose
ltry = PosteriorLSTAR(data, newbeta);
u = rand(d3,1)[1]
if u < exp(ltry - oldlik)
print([u exp(ltry - oldlik) ltry  oldlik], "\n")
tryPar = newbeta
end
end
oldlik = PosteriorLSTAR(data, tryPar);
StoreDraws[i,:] = tryPar;
end

Draws = StoreDraws[1:size(StoreDraws)[1],:];
``````

return Draws;
end

d = Normal(0,0.2)
dq = Normal(0,1)
theta = [2.0,0.0];
T = 50
q = rand(dq, T);
ra = quantile(q, [0.2, 0.8])
print(ra)
tr = LSTAR(q,theta, T);
x = ones(T,2);
x[:,2] = rand(dq, T);
y = 0.1*ones(T,1) + 0.3.*x[:, 2] + 1.5.*x[:,2].*tr + rand(d, T);

data = [y x q];
ini = [0.1;0.3;0.5;2.0;0.0;0.2];
#lik = PosteriorLSTAR(data, par);
mc = 10000
draws = mcmcLSTAR(mc, data, ini);
draws[:,4] = exp.(draws[:,4]);
draws[:,6] = exp.(draws[:,6]);
mean!([1. 1. 1. 1. 1. 1. ], draws)

Not sure what is going on but if you are forcing the loglikelihood to take on some specific value when a parameter moves out of bounds, why not just look for that specifically?

``````if rand() < exp(lokNew - lokOld) & lokNew  != -1e10
oldPar = newPar
end
``````

Note that all practical MCMC code works in log densities, where numerical overflow/underflow is much less acute of an issue. I see that you are doing this to some extent, but you may want to use somethng like

``````using Random
if randexp() < -(logNew - logOld)
...
``````

Is

``````randexp() < (lokNew - lokOld)
``````

really equivalent to

``````rand() < exp(lokNew - lokOld)
``````

?
`log(rand())` is always a negative number whereas `randexp()` is always a positive number?

1 Like

Thanks, I missed the `-`, fixed now.

1 Like