scaled_exp(x,m,beta) = exp( (m-x)/beta)
gumbel(x,m,beta) = scaled_exp(x,m,beta)* exp(-scaled_exp(x,m,beta)) /beta
lognormal(x,mu,sigma) = exp(-(log(x)-mu)**2/(2*sigma**2))/(x*sigma*sqrt(2*pi))

Gumbel_P(x,m,beta) = 1-exp(-scaled_exp(x,m,beta))
# Note: Gumbel_P may have rounding error problems as exp(-epsilon) is
#	very close to 1.  
# We can approximate it in these cases better using exp(-epsilon) approx 1-epsilon.
# (with double precision arithmetic, we have enough accuracy it
# doesn't matter which way we do it for this problem)
Gumbel_P_approx(x,m,beta)= scaled_exp(x,m,beta)


lognormal_P(x,mu,sigma) = 0.5*erfc( (log(x)-mu)/(sigma*sqrt(2)))

set ylabel 'probability'
set xlabel "ORF length (codons)'




set title 'Uzilov, model1'

unset logscale x
set logscale y
m=60; beta=20;
fit  gumbel(x,m,beta) 'uzilov-model1.hist' using 1:(1.e-4*$2) via m,beta
plot [10:400] gumbel(x,m,beta), 'uzilov-model1.hist' using 1:(1.e-4*$2)

print Gumbel_P(388, m, beta), Gumbel_P_approx(388, m, beta)

pause -1 "press return to continue"

mu=log(50); sigma=0.3;
fit lognormal(x,mu,sigma) 'uzilov-model1.hist' using 1:(1.e-4*$2) via mu,sigma

set logscale xy
plot [10:400] lognormal(x,mu,sigma), 'uzilov-model1.hist' using 1:(1.e-4*$2)

print lognormal_P(388, mu, sigma)

pause -1 "press return to continue"




set title 'Uzilov, model2'

unset logscale x
set logscale y
m=60; beta=20;
fit  gumbel(x,m,beta) 'uzilov-model2.hist' using 1:(1.e-4*$2) via m,beta
plot [10:400] gumbel(x,m,beta), 'uzilov-model2.hist' using 1:(1.e-4*$2)

print Gumbel_P(388, m, beta), Gumbel_P_approx(388, m, beta)

pause -1 "press return to continue"

mu=log(60); sigma=4;
fit lognormal(x,mu,sigma) 'uzilov-model2.hist' using 1:(1.e-4*$2) via mu,sigma

set logscale xy
plot [10:400] lognormal(x,mu,sigma), 'uzilov-model2.hist' using 1:(1.e-4*$2)

print lognormal_P(388, mu, sigma)

pause -1 "press return to continue"