samp = data.frame(g=numeric(0), e=numeric(0), l=numeric(0), side=logical(0),p=numeric(0))
K = 0.001
t = qnorm (1-K)
herit = 0.5
sdg = sqrt(herit)
sde = sqrt(1-herit)
g = rnorm (1) * sdg
e = rnorm (1) * sde
l = g + e
side = l >t
p.weight = dnorm(g/sdg)*dnorm(e/sde)*(side*(1-K) + (!side)*K)


for (i in 1:10^6){
# next two lines are a draw from the candidate distribution
	g.new = g + rnorm(1)*0.1
	e.new = e + rnorm(1)*0.1

	l.new = g.new + e.new
	side.new = l.new>t

# Ratio of probabilities of new and old sample
	p.new = dnorm(g.new/sdg)*dnorm(e.new/sde)*(side.new*(1-K) + (!side.new)*K)

	
	if (runif(1)< p.new/p.weight){
		g = g.new
		e = e.new
		l=l.new
		side=side.new
		p.weight = p.new
	}
	if (i%% 500 == 0){
		cat (i,g,e,l,side,p.weight,"\n")
		samp = rbind(samp, data.frame(g,e,l,side,p.weight))
	}
}