# NOT RUN {
## simulate some fake data
G=100 ## for demonstration only. Normally, G should be much larger
sdncp=1.3
n1=n2=5
df=n1+n2-2
set.seed(54457704)
x=runif(G,1,G)
f=function(x)sin(x*pi/1000)+1
Pi.i=1/(1+exp(f(x)))
Z.i=rbinom(G,1,1-Pi.i)
t0.i=rt(G,df)
ncp.i=rnorm(G,0,sdncp)
t1.i=rt(G,df, ncp.i)
t.i=ifelse(Z.i==0,t0.i,t1.i)
## fit model
(plfit=penLik.EMNewton(t.i, x, df, spar=10^seq(0,8,length=30),plotit=FALSE))
(plfit0=scaledTMix.null(t.i, df))
# }
# NOT RUN {
plot(plfit)
plot(t.i, plfit$lfdr, pch='.')
lines(sort(t.i), plfit0$lfdr[order(t.i)], col=2, lwd=3)
# }
Run the code above in your browser using DataLab