Load script
Nlocus=1000
# Effect of s
###############
meanshap <- 0.01
shapesel <- 2
set.seed(1)
#shap <- c(0.0002,0.0005,0.00075,0.001, 0.005, 0.01)
shap <- sort(rgamma(Nlocus, shape=shapesel, scale=meanshap/shapesel))
mu <- 0.0001
hsim <- 0.25
c <- 0.5
q <- (mu*exp(-shap/2))/(1-exp(-(shap/2)* (1+hsim*c)))
shap <- shap[q<1]
q <- q[q<1]
hist(shap)

plot(shap, q, type="l", las=1)

plot(shap, q*(1-q), type="l", las=1, xlab="Haploid selection coefficient", xlim=c(0, 0.05), ylim=c(0, 0.25))

(cest <- mean(c*shap)/mean(shap))
## [1] 0.5
# Effect of h
###############
shapeh <- 2
h <- 0.25
a <- shapeh
b <- a*(1-h)/h
hsim <- sort(rbeta(n=Nlocus, shape1=a, shape2=b))
mu <- 0.0001
c <- 0.5
shap=0.01
q <- (mu*exp(-shap/2))/(1-exp(-(shap/2)* (1+hsim*c)))
plot(hsim, q, type="l", las=1, ylim=c(0, 0.02))

hist(hsim)

plot(hsim, q*(1-q), type="l", las=1, ylim=c(0, 0.02), xlab="Level of dominance")

(cest <- mean(c*shap)/mean(shap))
## [1] 0.5
# Effect of c
###############
shapec <- 2
c <- 0.5
a <- shapec
b <- a*(1-c)/c
csim <- sort(rbeta(n=Nlocus, shape1=a, shape2=b))
mu <- 0.0001
h <- 0.25
shap<-0.01
q <- (mu*exp(-shap/2))/(1-exp(-(shap/2)* (1+h*csim)))
hist(csim)

plot(hsim, q, type="l", las=1, ylim=c(0, 0.02))

plot(csim, q*(1-q), type="l", las=1, ylim=c(0, 0.02), xlab="Ratio of fitness effect in diploids vs haploids")

(cest <- mean(csim*shap)/mean(shap))
## [1] 0.5078183