Stochastic simulation error

Hello,

I am newbie in Julia, trying to run a MCMC algorithm, but I getting some complaints about the parameters of the gamma distribution (which are supposed to be > 0), but they are, so I don’t understand this complaints. My code is below:

using Random, Distributions

# Data

M=52
NN=30
n=26
m=34
x = [49., 5., 17., 2., 39., 84., 7., 0., 35., 36., 1400., 5., 34., 15., 11.,
	2., 1., 39., 8., 101., 2., 148., 1., 68., 31., 1., 20., 118., 91.,
	427.]
dc = [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.,
	1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 0., 0.]
y = [0., 35., 50., 11., 25., 16., 36., 27., 19., 17., 7., 11., 2., 82., 24.,
	70., 15., 16., 50., 22., 45., 18., 4., 1., 40., 57., 0., 1., 20.,
	35., 82., 31., 40., 9., 66., 20., 77., 2., 26., 32., 13., 56., 2.,
	9., 4., 30., 3., 26., 4., 45., 25., 5.]
z = [15., 3., 624., 46., 127., 61., 1350., 312., 24., 10., 1024., 39., 730.,
	136., 1379., 1., 836., 60., 1140., 1153., 54., 47., 0., 43., 971.,
	868., 44., 780., 51., 710., 663., 253., 147., 51., 479., 322., 442.,
	65., 419., 362., 64., 228., 65., 264., 25., 193., 196., 63., 12.,
	103., 60., 43.]
dt = [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1.,
	0., 0., 1., 1., 1., 1., 0., 0., 1., 0., 1., 0., 0., 1., 1., 1., 0.,
	1., 0., 1., 0., 0., 1., 1., 1., 0., 1., 0., 0., 1., 1., 0., 0., 0.]

#defining vectors and matrices:
lambda=zeros(Float64,1000)
tau=zeros(Float64,1000)
p=zeros(Float64,1000)
theta=zeros(Float64,(1000, NN))
eta=zeros(Float64,(1000, M))


# setting initial values:
lambda[1]=0.1
tau[1]=0.1
p[1]=3
theta[1,:].=0.1
eta[1,:].=0.1
alpha=0.01
beta=0.01
gama=0.01
delta=0.01


# simulation
for k in 2:1000
	lambda[k]=rand(Gamma(2+NN*p[k-1]+M*p[k-1]+alpha,2+beta+sum(theta[k-1,:])+sum(eta[k-1,:])))
	tau[k]=rand(Gamma(m+gama, delta+sum(eta[k-1,:].*z)))
	p[k]=rand(Gamma(2+lambda[k]*((prod(theta[k-1,:])*prod(eta[k-1,:]))^(1/(M+NN)))^(0.999*(M+NN)),1/gamma(0.999)))
	for i in 1:NN
	theta[k,i]=rand(Gamma(p[k]+dc[i],lambda[k]+x[i]))
	end
	for j in 1:M
		eta[k,j]=rand(Gamma(p[k]+dt[j],lambda[k]+y[j]+tau[k]*z[j]))
	end
end

It complains that some of the gamma distributions within the for have parameters not satisfying the definition, but we compute the initial parameters manually they are greater than zero. I can’t see why it is not working. I would appreciate any help. thanks.

Hi @jailton, welcome to the Julia forum. Have a read of PSA: make it easier to help you and edit your post to have the code-blocks commented.