Nothing
# generates m pairs of X and Y and independent m X samples (called Xprime) and n Y samples (called Yprime)
sim.partially.matched=function(m, n.x, n.y, distr=c("normal","logistic","student","mixnormal","gamma","lognormal","beta","uniform","hybrid1","hybrid2","doublexp"), params, seed){
set.seed(seed)
N=m+n.x+n.y
distr=match.arg(distr)
if (distr %in% c("normal","student","logistic","lognormal","doublexp")) {
loc.2=params["loc.2"]
scale.2=params["scale.2"]
rho=params["rho"]; if (is.na(rho)) rho=0
if (distr=="normal") {
dat.all=mvtnorm::rmvnorm(n = N, mean=c(0,loc.2), sigma=matrix(c(1,scale.2*rho,scale.2*rho,scale.2^2),nrow=2))
} else if (distr=="lognormal") {
dat.all=exp(mvtnorm::rmvnorm(n = N, mean=c(0,0), sigma=matrix(c(1,scale.2*rho,scale.2*rho,scale.2^2),nrow=2)) )
dat.all[,2]=dat.all[,2]+loc.2
} else if (distr=="student") {
dat.all=mvtnorm::rmvt (n = N, delta=c(0,loc.2), sigma=matrix(c(1,scale.2*rho,scale.2*rho,scale.2^2),nrow=2), df=4)
} else if (distr=="logistic") {
dat.all=rbilogistic(N, 0, loc.2, scale.1=1, scale.2=scale.2, rho)
} else if (distr=="doublexp") { # double exponential also known as Laplace
dat.all=rbidoublexp(N, 0, loc.2, scale.1=1, scale.2=scale.2, rho)
}
} else if (distr=="mixnormal") {
p.1=params["p.1"]
p.2=params["p.2"]
sd.n=params["sd.n"]; if (is.na(sd.n)) sd.n=1
sd1=params["sd1"]; if (is.na(sd1)) sd1=0.5
sd2=params["sd2"]; if (is.na(sd2)) sd2=0.5
dat.all=rnorm(N, sd=sd.n) + cbind(rmixnorm (N, mix.p=p.1, mu1=0, mu2=2, sd1=sd1, sd2=sd2), rmixnorm (N, mix.p=p.2, mu1=0, mu2=2, sd1=sd1, sd2=sd2))
} else if (distr=="gamma") {
loc.2=params["loc.2"]
shape.1=params["shape.1"]
shape.2=params["shape.2"]
rate.1=params["rate.1"]
rate.2=params["rate.2"]
rho=params["rho"]
dat.all=rbigamma(n=N, shape.1, shape.2, rate.1, rate.2, rho=rho)
dat.all[,2]=dat.all[,2]+loc.2
} else if (distr=="beta") {
shape1.1=params["shape1.1"]
shape2.1=params["shape2.1"]
shape1.2=params["shape1.2"]
shape2.2=params["shape2.2"]
dat.all=cbind(rbeta(N,shape1.1,shape2.1),rbeta(N,shape1.2,shape2.2))
} else if (distr=="uniform") {
loc.2=params["loc.2"]
dat.all=cbind(runif(N),runif(N))
dat.all[,2]=dat.all[,2]+loc.2
} else if (distr=="hybrid1") {
# X gamma Y logistic
loc.2=params["loc.2"]
toggle=params["toggle"]
if (toggle==0) {
dat.all=cbind(rgamma(N,shape=2,rate=0.5), rlogis(N)+loc.2)
} else if (toggle==1) {
dat.all=cbind(rlogis(N)+loc.2, rgamma(N,shape=2,rate=0.5))
}
} else if (distr=="hybrid2") {
# X lognormal Y logistic
loc.2=params["loc.2"]
toggle=params["toggle"]
if (toggle==0) {
dat.all=cbind(exp(rnorm(N)), rlogis(N)+loc.2)
} else if (toggle==1) {
dat.all=cbind(rlogis(N)+loc.2, exp(rnorm(N)))
}
}
index.1=if(n.x==0) NULL else m+1:n.x
index.2=if(n.y==0) NULL else m+n.x+1:n.y
list(X=dat.all[1:m,1], Y=dat.all[1:m,2],
Xprime=dat.all[index.1,1], Ymissing=dat.all[index.1,2],
Xmissing=dat.all[index.2,1], Yprime=dat.all[index.2,2])
}
r2sample=function(m, n, distr=c("normal","logistic","student","mixnormal"), params, seed){
distr=match.arg(distr)
dat=sim.partially.matched(m=m, n.x=0, n.y=n, distr, params, seed)[c("X","Yprime")]
names(dat)=c("X","Y")
dat
}
sim.paired.with.replicates=function(m, meanRatio, sdRatio, within.sd, type, hyp, distr, seed){
set.seed(seed)
ratios=round(rnorm(n=m,mean=meanRatio,sd=sdRatio))
tmp=sim.partially.matched(m=m,n.x=0,n.y=0,distr=distr,params=c(loc.2=0,rho=0.5,scale.2=1),seed=seed)
mother=list()
if(type==1) {
# infant seq and mother sequences are exchangeable
infant=numeric(m)
for (i in 1:m){
tmp.1=tmp$Y[i]+rnorm(n=ratios[i]+1, sd=within.sd)
if (hyp==0) {
# null: choose a random one to be infant: 1 is a random one
sel=1
} else if (hyp==1) {
# filtering model: choose smaller ones to be infant with higher prob, to be consistent with the direction of 1-sided test
tmp.1=sort(tmp.1)
sel=rbern(1, 1/(1:length(tmp.1))^1.05, generalized=TRUE)
} else if (hyp==2) {
# location shift model
sel=1
tmp.1[sel]=tmp.1[sel]-.5
} else stop("wrong hyp")
infant[i]=tmp.1[sel]
mother=c(mother, list(tmp.1[-sel]))
}
} else if (type==2) {
# null only
# infant seq and mother sequences are not exchangeable
infant=tmp$X
for (i in 1:m){
tmp.1=tmp$Y[i]+rnorm(n=ratios[i], sd=within.sd)
if (hyp==0) {
infant[i]=infant[i]+rnorm(n=1, sd=within.sd)
mother=c(mother, list(tmp.1))
} else stop("for type 2, only hyp 0 is implemented")
}
} else stop ("wrong type")
list(X=infant, Y=mother)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.