# NOT RUN {
library(DIME)
# generate simulated datasets with underlying exponential-normal components
N1 <- 1500; N2 <- 500; K <- 4; rmu <- c(-2.25,1.50); rsigma <- c(1,1);
rpi <- c(.05,.45,.45,.05); rbeta <- c(12,10);
set.seed(1234);
chr1 <- c(-rgamma(ceiling(rpi[1]*N1),shape = 1,scale = rbeta[1]),
rnorm(ceiling(rpi[2]*N1),rmu[1],rsigma[1]),
rnorm(ceiling(rpi[3]*N1),rmu[2],rsigma[2]),
rgamma(ceiling(rpi[4]*N1),shape = 1,scale = rbeta[2]));
chr2 <- c(-rgamma(ceiling(rpi[1]*N2),shape = 1,scale = rbeta[1]),
rnorm(ceiling(rpi[2]*N2),rmu[1],rsigma[1]),
rnorm(ceiling(rpi[3]*N2),rmu[2],rsigma[2]),
rgamma(ceiling(rpi[4]*N2),shape = 1,scale = rbeta[2]));
chr3 <- c(-rgamma(ceiling(rpi[1]*N2),shape = 1,scale = rbeta[1]),
rnorm(ceiling(rpi[2]*N2),rmu[1],rsigma[1]),
rnorm(ceiling(rpi[3]*N2),rmu[2],rsigma[2]),
rgamma(ceiling(rpi[4]*N2),shape = 1,scale = rbeta[2]));
# analyzing only chromosome 1 and chromosome 3
data <- list(chr1,chr3);
# Fitting a GNG model only
bestGng <- gng.fit(data,K=2);
# Plot the estimated GNG model imposed on the histogram of the observed data
hist(unlist(data),freq=FALSE,breaks=100,xlim=c(-20,20))
gng.plot.mix(bestGng,resolution=1000,new.plot=FALSE,col="red");
# }
Run the code above in your browser using DataLab