rm(list=ls())

####  VERY SIMPLE

x=seq(0,20,length=100)
k=0.67   ## 1.8
b=30     ## 35
y=exp(k*x)/(k*(b+exp(k*x)))


pdf("Plot.pdf")
plot(x,y,type="l",col="purple",lwd=6,xlim=c(0,11),xaxt="n",yaxt="n",
                  ylab="",xlab="")
title(xlab="x", line=2, cex.lab=2)
title(ylab=expression(pi(x,beta)),line=1, cex.lab=2)
dev.off()

## k
## [1] 0.6
## b
## [1] 30


b0=3          ### .2,3,.01, .5
b1=.3         ### 2, 10,20, .3

y=exp(b0+b1*x)/(1+exp(b0+b1*x))
plot(x,y,type="l",col="red",lwd=2)