# Numerical Problems

**URL:** <https://discourse.julialang.org/t/numerical-problems/37364>\
**Category:** General Usage\
**Tags:** question\
**Created:** [April 10, 2020, 6:57pm UTC](https://discourse.julialang.org/t/numerical-problems/37364 "2020-04-10T18:57:00Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Andres\_Gonzalez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andres_gonzalez/32/13993_2.png) [@Andres\_Gonzalez](https://discourse.julialang.org/u/Andres_Gonzalez)\
**Post date:** [April 10, 2020, 6:57pm UTC](https://discourse.julialang.org/t/numerical-problems/37364/1 "2020-04-10T18:57:00Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [April 10, 2020, 7:11pm UTC](https://discourse.julialang.org/t/numerical-problems/37364/2 "2020-04-10T19:11:37Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Andres\_Gonzalez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andres_gonzalez/32/13993_2.png) [@Andres\_Gonzalez](https://discourse.julialang.org/u/Andres_Gonzalez)\
**Post date:** [April 10, 2020, 7:11pm UTC](https://discourse.julialang.org/t/numerical-problems/37364/3 "2020-04-10T19:11:55Z")

</div>

Yes

---

<div class="post-metadata">

**Author:** ![Andres\_Gonzalez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andres_gonzalez/32/13993_2.png) [@Andres\_Gonzalez](https://discourse.julialang.org/u/Andres_Gonzalez)\
**Post date:** [April 10, 2020, 7:13pm UTC](https://discourse.julialang.org/t/numerical-problems/37364/4 "2020-04-10T19:13:14Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [April 10, 2020, 7:13pm UTC](https://discourse.julialang.org/t/numerical-problems/37364/5 "2020-04-10T19:13:38Z")

</div>

> [@Andres\_Gonzalez](#):
>
> exp(-1.0e10 - logOld)

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

---

<div class="post-metadata">

**Author:** ![Andres\_Gonzalez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andres_gonzalez/32/13993_2.png) [@Andres\_Gonzalez](https://discourse.julialang.org/u/Andres_Gonzalez)\
**Post date:** [April 10, 2020, 7:14pm UTC](https://discourse.julialang.org/t/numerical-problems/37364/6 "2020-04-10T19:14:09Z")

</div>

logOld is a Float64 number close to -617.0

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [April 10, 2020, 7:15pm UTC](https://discourse.julialang.org/t/numerical-problems/37364/7 "2020-04-10T19:15:36Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Andres\_Gonzalez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andres_gonzalez/32/13993_2.png) [@Andres\_Gonzalez](https://discourse.julialang.org/u/Andres_Gonzalez)\
**Post date:** [April 10, 2020, 7:17pm UTC](https://discourse.julialang.org/t/numerical-problems/37364/8 "2020-04-10T19:17:00Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Andres\_Gonzalez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andres_gonzalez/32/13993_2.png) [@Andres\_Gonzalez](https://discourse.julialang.org/u/Andres_Gonzalez)\
**Post date:** [April 10, 2020, 7:25pm UTC](https://discourse.julialang.org/t/numerical-problems/37364/9 "2020-04-10T19:25:30Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [April 10, 2020, 7:30pm UTC](https://discourse.julialang.org/t/numerical-problems/37364/10 "2020-04-10T19:30:47Z")

</div>

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?

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

```

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 11, 2020, 7:36am UTC](https://discourse.julialang.org/t/numerical-problems/37364/11 "2020-04-11T07:36:33Z")

</div>

@Andres_Gonzalez, please [quote your code](https://discourse.julialang.org/t/psa-how-to-quote-code-with-backticks/7530), 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

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

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 13, 2020, 4:17am UTC](https://discourse.julialang.org/t/numerical-problems/37364/12 "2020-04-13T04:17:19Z")

</div>

Is

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

```

really equivalent to

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

```

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

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 13, 2020, 4:29am UTC](https://discourse.julialang.org/t/numerical-problems/37364/13 "2020-04-13T04:29:42Z")

</div>

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