Numerical Problems

General Usage
1

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

2

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

3

Yes

4

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

5

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

6

logOld is a Float64 number close to -617.0

7

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

8

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

9

using Plots
using StatsBase, StatsFuns, DataFrames
using Distributions
using NLsolve
using CSV
using LinearAlgebra
using Statistics
#Pkg.free(“Dataframes”, “nl/readtable”)

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)

10

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
11

@Andres_Gonzalez, please quote your code, this makes it much easier to read.

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

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
13

Thanks, I missed the -, fixed now.

1 Like