#' Returns the jth selected person in month m when the random number drawn is i.
#'
#' @param i An integer
#' @param j An integer
#' @param m An integer
#' @return The \code{j}th selected person in month \code{m} when the random number drawn is \code{i}.
#' @examples
#' samplerule(1,1:160,1)
samplerule<-function(i,j,m){((5*(m-1)+
5*(i-1)+
5000*((j-1)%/%5)+((j-1)%%5)+5*((j-1)%/%100)+40*((j-1)%/%400)))%%100000+1}
#' Return atomic rotation group 1 for seed i (1000 possible values)
#'
#' @param i An integer between 1 and 1000
#' @return The atomic rotation group 1 for seed i
#' @examples
#' sampleruleS(1,1:160,1)
#' identical(sampleruleS(i0=1),samplerule(i=1,j=1:100,m=1))
#' identical(sampleruleS(i0=16),samplerule(i=1,1:100,m=16))
#'
#' (checkrule<-function(){
#' i=sample(1:1000,1);m=sample(1:85,1);h=sample(1:8,1);
#' return(list(i=i,m=m,h=h,check=identical(sampleruleS(i0=i+(m-1)+c(0:3,12:15)[9-h]),samplerule(i,(8-h)*100+(1:100),m))))})()
#' identical((function(i,m){c(sapply((m+i-2)+c(1:4,13:16),sampleruleS))})(1,1),samplerule(1,1:800,1))
sampleruleS<-function(i0){((5*(i0-1)+5000*((0:99)%/%5)+((0:99)%%5)+5*((0:99)%/%100)))%%100000+1}
sampleruleH<-
function(i,h,m){
j=c(outer((8-h)*100,(1:100),"+"));
(((m-1)+
(i-1)+
1000*((j-1))+((j-1)%/%20)+8*((j-1)%/%80)))%%20000+1}
sampleruleR<-function(i,h,m){sampleruleS(i0=i+(m-1)+c(0:3,12:15)[9-h])}
misH<-function(i,m,misi,syntheticcpspopsHA){plyr::aaply(syntheticcpspopsHA[,sampleruleH(i,(8-misi)*20+(1:20),m),,m],c(1,3),sum)}
mis<-function(i,m,misi,syntheticcpspopsA){plyr::aaply(syntheticcpspopsA[,samplerule(i,(8-misi)*100+(1:100),m),,m],c(1,3),sum)}
misS<-function(i,m,misi,syntheticcpspopsS){plyr::aaply(syntheticcpspopsA[,samplerule(i,(8-misi)*100+(1:100),m),,m],c(1,3),sum)}
#system.time(misestimates<-plyr::maply(expand.grid(i=1:3,m=1:85,misi=1:8),mis,syntheticcpspopsA=syntheticcpspopsA,.parallel=TRUE))
#system.time(misestimates<-plyr::maply(expand.grid(i=1:3,m=1:85,misi=1:8),mis,syntheticcpspopsA=syntheticcpspopsA))
#system.time(misestimates<-plyr::maply(expand.grid(i=1:20,m=1:85,misi=1:8),misH,syntheticcpspopsHA=syntheticcpspopsHA,.parallel=TRUE))
#system.time(misestimates<-plyr::maply(expand.grid(i=1:20,m=1:85,misi=1:8),misH,syntheticcpspopsHA=syntheticcpspopsHA))
#mis(1,2,1,syntheticcpspopsHA)
roundv<-function(n,x){y<-floor(n*x);y[length(x)]<-n-sum(y[-length(x)]);y}
rectif<-function(varr,prob_1){
N<-length(varr)
eff_0=unlist(table(varr))[names(prob_1)]
eff_1<-roundv(N,prob_1)
eff<-eff_1-eff_0
qui<-names(eff)[eff<0]
remps<-names(eff)[eff>0]
varr[sample((1:N)[is.element(varr,qui)],sum(eff[eff>0]))]<-rep(remps,eff[eff>0])
return(as.factor(varr))}
changevar<-function(varr,prob2=NULL,prob,iterative.synthetic.model="Independent"){
var2<-varr
zz=unique(varr)
for(i in zz){
n<-sum(varr==i)
if(iterative.synthetic.model=="2dorder"){
var2[varr==i]<-rep(zz,roundv(n,prob2[i,]))[order(runif(n))]}
if(iterative.synthetic.model=="2dorder-indexdependent"){
var2[varr==i]<-rep(zz,roundv(n,prob2[i,]))}}
rectif(var2,prob)
}
syntheticcpsdataset<-
function(Totals,
crossTotals=NULL,
cluster.size=5,
cluster.all.sample.rate=1/10,
cluster.single.sample.size=20,
seed=1,
iterative.synthetic.models=c("Independent","2dorder","2dorder-indexdependent")){
if(is.null(Totals)){
Totals=plyr::aaply(1:85,1,function(x){y=runif(3);y/sum(y)})}
L=nrow(Totals)
nb.samples<-L+15
cluster.single.sample.rate <-cluster.all.sample.rate/nb.samples
N<-cluster.size*
cluster.single.sample.size/
cluster.all.sample.rate*
(nb.samples)
nb.clusters<-N/cluster.size
sample.rate<-8*cluster.single.sample.size/nb.clusters
set.seed(seed)
prob <- Totals/apply(Totals,1,sum)
prob2<-plyr::aaply(crossTotals,1:2,function(x){matrix(x/sum(x))})
dimnames(prob2)<-dimnames(crossTotals)
Hmisc::label(prob2)<-"P(PumlrR(t+1)=j|PumlrR(t)=i)"
creevar<-function(nmod,prob=NULL){
modal <- nmod
if( is.integer(nmod)&length(nmod)==1){modal <- 1:nmod}
as.factor(sample(modal,N,replace=TRUE, prob=prob))}
creevarpumlr<-function(prob,N=100000){
return(rep(names(prob),roundv(N,prob))[order(runif(N))])}
Popu<-data.frame(
hrlongid=rep(1:(N/cluster.size),each=cluster.size),
pulineno=rep(1:cluster.size,N/cluster.size),
pwsswgt=rep(1/sample.rate,N),
pumlrR =creevarpumlr(prob[1,],N),
pumlrRlag =creevarpumlr(prob[1,],N))
LL<-lapply(iterative.synthetic.models,function(j){
list.tablespop<-list(Popu)
Popu2=Popu
for (i in (2:L)){
Popu2<-data.frame(
hrlongid=Popu$hrlongid,
pulineno=Popu$pulineno,
pwsswgt=Popu$pwsswgt,
pumlrRlag =Popu$pumlrR,
pumlrR =changevar(Popu2$pumlrR,prob2[i-1,,],prob[i,],iterative.synthetic.model=j))
list.tablespop<-c(list.tablespop,list(Popu2))}
names(list.tablespop)<-dimnames(Totals)$month
Hmisc::label(list.tablespop,paste0('Synthetic "CPS populations" datasets, ',j," method"))
list.tablespop})
names(LL)<-iterative.synthetic.models
Hmisc::label(LL)<-'Synthetic "CPS populations" datasets'
LL
}
pumlrR <- function(i){replace(x<-as.character(-1+2*is.element(i,c("1","2"))+is.element(i,as.character(3:4))),x=="-1","_1")}
#This function modify the dataframes of the list list.tables
#and returns a list of modified tables
Chargetables <- function(list.tables){
L<-lapply(list.tables,
function(dff){
names(dff)<-tolower(names(dff))
colnames<-intersect(tolower(c("GESTFIPS","PUMLR","HRMIS","puhhmem",
'prpertyp',"pumlrlag",
'pesex' ,
'prwtrace',
'pehspnon',
'PRBLNONB',
'GESTREC')),names(dff))
colnamesnum<-setdiff(tolower(names(dff)),colnames)
dff[colnamesnum] <- lapply(dff[colnamesnum] ,as.numeric)
dff[colnames] <- lapply(dff[colnames] ,as.factor)
dff<-transform(dff,pwsswgt=pwsswgt/10000)
dff<-transform(dff,pwcmpwgt=pwcmpwgt/10000)
if(is.element("hwniwgt",colnames)){dff<-transform(dff,hwniwgt=hwniwgt/10000)}
dff$pumlrR <- as.factor(pumlrR(dff$pumlr))
if(is.element("pumlrlag",names(dff))){dff$pumlrRlag <- as.factor(pumlrR(dff$pumlrlag))
}
#assign(tables.entree[[i]],dff, envir = .GlobalEnv)
return(dff)})
names(L)<-names(list.tables)
return(L)}
#This function modify the dataframes of the list list.tablespop
#and returns a list of modified tables
syntheticcpspopsbf<-function(syntheticcpspops){
plyr::llply(syntheticcpspops,function(l){
plyr::llply(l,function(ll){
ll[["pumlrR"]][sample((1:nrow(ll))[ll$"pumlrR"=="1"],ceiling(0.2*sum(ll$"pumlrR"=="1")))]<-"0"
ll
})
})}
Chargetablespop <- function(list.tablespop){
L<-mclapply(list.tablespop,
function(dff){
names(dff)<-tolower(names(dff))
colnames<-intersect(tolower(c("GESTFIPS","PUMLR","puhhmem",
'prpertyp',"pumlrlag",
'pesex' ,
'prwtrace',
'pehspnon',
'PRBLNONB',
'GESTREC')),names(dff))
colnamesnum<-setdiff(tolower(names(dff)),colnames)
dff[colnamesnum] <- lapply(dff[colnamesnum] ,as.numeric)
#dff[colnames] <- lapply(dff[colnames] ,as.character, stringsAsFactors=FALSE)
dff[colnames] <- lapply(dff[colnames] ,as.factor)
dff<-transform(dff,pwsswgt=pwsswgt/10000)
dff<-transform(dff,pwcmpwgt=pwcmpwgt/10000)
dff<-transform(dff,hwniwgt=hwniwgt/10000)
dff$pumlrR <- as.factor(pumlrR(dff$pumlr))
dff$pumlrRlag <- as.factor(pumlrR(dff$pumlrlag))
#assign(paste(tables.entree[[i]],"_pop",sep=""),dff, envir = .GlobalEnv)
return(dff)})
names(L)<-names(list.tablespop)
return(L)}
ChargetablespopA<-function(list.tablespop){
plyr::laply(list.tablespop,
function(l){
model.matrix(~0+pumlrR,l)})}
syntheticccpspopHAf<-function(list.tablespop){
list.tablespopA<-
array(unlist(lapply(list.tablespop,function(l){
apply(array(
cbind(pumlrR_n0=l$pumlrR=="0",
pumlrR_n1=l$pumlrR=="1",
pumlrR_n_1=l$pumlrR=="_1")*1,c(5,20000,3)),2:3,sum)})),c(20000,3,length(list.tablespop)))
dimnames(list.tablespopA)[2:3]<-list(c("0","1","_1"),names(list.tablespop))
names(dimnames(list.tablespopA))<-c("h","y","m")
list.tablespopA}
#sapply(1:3,list.tablespopAf)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.