Nothing
# Fonctions internes :
#' @param obj_fct une liste contenant tous les objets de l'environnement courant de
#' calcul (current evaluation environment) dans la fonction appelante
######################################################################################################
#' Preparation de stratumID pour la sortie
#'
#' Ajoute a straumIDnoc les observations de la strate certain + format facteur.
#'
#' @return stratumID tel que definit dans le template stratumID
#'
get_stratumID <- function(obj_fct)
{
# Pour tirer de obj_fct les variables dont on a besoin ici
certain <- obj_fct$certain
stratumIDnoc <- obj_fct$stratumIDnoc
N <- obj_fct$N
L <- obj_fct$L
# Insertion des "certain" au bon endroit dans stratumIDnoc
if(is.null(certain)){
stratumID <- stratumIDnoc
} else {
stratumID <- rep("certain", N)
stratumID[-certain] <- stratumIDnoc
}
# Pour faire de ma sortie un facteur
lev <- if(is.null(certain)) 1:L else c(1:L, "certain")
stratumID <- factor(stratumID, levels=lev)
stratumID ## un seul objet retourne
}
# Creee et verifiee le 26 septembre 2012.
######################################################################################################
#' Wrapper R pour la fonction C get_EYsC
#'
get_EYs <- function(xs, Ns, ps, obj_fct)
{
resuC <- .C(get_EYs_C, as.double(xs), as.integer(Ns), as.integer(obj_fct$nmodel),
as.double(obj_fct$beta), as.double(obj_fct$sig2), as.double(ps), as.double(obj_fct$gamma),
as.double(obj_fct$epsilon), as.double(obj_fct$EX),
EYs = as.double(0), EXs = as.double(0), phis = as.double(0), NAOK = TRUE)
resuC$EYs ## un seul objet retourne
}
# Creee et verifiee avec valeurs par defaut le 26 septembre 2012.
######################################################################################################
#' Wrapper R pour la fonction C strata_bh_opti_C
#'
strata_bh_opti <- function(bhfull, takeallin, takeall.adjust, dotests, obj_fct)
{
L <- obj_fct$L
Nnoc <- obj_fct$Nnoc
minNh <- if (!dotests) 0 else obj_fct$minNh
resuC <- .C(strata_bh_opti_C, as.double(obj_fct$xnoc), as.integer(Nnoc), as.double(bhfull),
as.integer(L), as.integer(obj_fct$takenone), as.integer(takeallin), as.integer(obj_fct$Nc),
as.double(obj_fct$EYc), as.double(obj_fct$q1), as.double(obj_fct$q2), as.double(obj_fct$q3),
as.integer(obj_fct$nmodel), as.double(obj_fct$beta), as.double(obj_fct$sig2), as.double(obj_fct$ph),
as.double(obj_fct$gamma), as.double(obj_fct$epsilon), as.double(obj_fct$EX), as.double(obj_fct$EX2),
as.integer(obj_fct$findn), as.integer(obj_fct$n), as.double(obj_fct$CV), as.double(obj_fct$rhL),
as.double(obj_fct$bias.penalty), as.integer(takeall.adjust), as.integer(dotests), as.integer(minNh),
## Pour les elements en sortie
NhOK = as.integer(0), nhOK = as.integer(0),
phih = as.double(rep(0,L)), psih = as.double(rep(0,L)), gammah = as.double(rep(0,L)),
ah = as.double(rep(0,L)), U2 = as.double(0), U = as.double(0), V = as.double(0),
stratumIDnoc = as.integer(rep(0,Nnoc)), Nh = as.integer(rep(0,L)),
EYh = as.double(rep(0, L)), VYh = as.double(rep(0, L)), TY = as.double(0), TAY = as.double(0),
nhnonint = as.double(rep(0, L)), takeallout = as.integer(0), nh = as.double(rep(0, L)),
opti.nhnonint = as.double(0), opti.nh = as.double(0), NAOK = TRUE)
resuC[names(resuC) != ""] ## une liste retournee
}
# Creee et verifiee avec valeurs par defaut le 28 septembre 2012.
######################################################################################################
#' Wrapper R pour la fonction C get_RRMSE_C
#'
get_RRMSE <- function(nhcalcul, obj_fct)
{
resuC <- .C(get_RRMSE_C, as.double(obj_fct$bias.penalty), as.double(obj_fct$TY), as.double(obj_fct$TAY),
as.integer(obj_fct$Nh), as.double(obj_fct$VYh), as.double(nhcalcul), as.double(obj_fct$rhL),
as.integer(obj_fct$L), as.integer(obj_fct$takenone),
RRMSE=as.double(0), NAOK = TRUE)
RRMSE <- resuC$RRMSE ## un seul objet retourne
}
# Creee le 1e octobre 2012
######################################################################################################
#' Wrapper R pour la fonction C get_n_C
#'
get_n <- function(nhcalcul, obj_fct)
{
resuC <- .C(get_n_C, as.double(nhcalcul), as.integer(obj_fct$L), as.integer(obj_fct$Nc),
n=as.double(0), NAOK = TRUE)
resuC$n ## un seul objet retourne
}
######################################################################################################
#' Fonction pour calculer trois stats uniquemnet utiles pour la sortie, pas utile pour des calculs :
#' RMSE (de la moyenne et non de la somme), relativebias et propbiasMSE
#'
#' @param RRMSE Dans le cas d'un n cible, il s'agit du critere d'optimiasation calcule sur les nh entier
#' (optinh), dans le cas d'un CV cible, RRMSE = NA (il doit etre calcule).
#' Pour la definition des autres parametres, voir le code C
#'
#' @return RMSE tel que definit dans le template RMSEtakenone
#' @return RRMSE tel que definit dans le template RMSEtakenone
#' @return relativebias tel que definit dans le template RMSEtakenone
#' @return propbiasMSE tel que definit dans le template RMSEtakenone
#'
get_stat_out <- function(obj_fct)
{
# Pour tirer de obj_fct les variables dont on a besoin ici
RRMSE <- obj_fct$RRMSE
nh <- obj_fct$nh
bias.penalty <- obj_fct$bias.penalty
TY <- obj_fct$TY
TAY <- obj_fct$TAY
Nh <- obj_fct$Nh
VYh <- obj_fct$VYh
rhL <- obj_fct$rhL
L <- obj_fct$L
takenone <- obj_fct$takenone
N <- obj_fct$N
# Calculs
# Dans le cas d'un CV cible, le RRMSE n'a pas ete calcule prealablement.
if (is.na(RRMSE)) RRMSE <- get_RRMSE(nhcalcul=nh, obj_fct = as.list(environment()))
RMSE <- RRMSE * TY / N
relativebias <- if (takenone == 0) 0 else (bias.penalty*TAY)/TY
propbiasMSE <- if (takenone == 0) 0 else ((bias.penalty*TAY)^2)/((N*RMSE)^2)
# Sortie des resultats
list(RMSE=RMSE, RRMSE=RRMSE, relativebias=relativebias, propbiasMSE=propbiasMSE) ## une liste retournee
# Note : Ici, dans le package, l'estimateur est : moyenne de Y (Ey) meme si
# dans l'article l'estimateur est : somme de Y (Ty),
# Lien entre les deux : Ey = Ty/N (voir definition de mean dans la sortie d'une fonction)
# RMSE de Ey = (RMSE de Ty)/N car Var(Ey) = Var(Ty/N) = Var(Ty)/N^2 et
# biais(Ey) = Ty/N pop - Ty/N ech = Tay/N = biais(Ty)/N
# Definition RMSE : sqrt(biais^2+variance) donc RMSE de Ey = sqrt((biais(Ty)^2+Var(Ty))/N^2) =
# sqrt((biais(Ty)^2+Var(Ty)))/N = (RMSE de Ty)/N
# Le facteur bias.penalty devant la biais est le meme pour les deux RMSE
# RRMSE de Ey = RRMSE de Ty
# car RRMSE de Ey = RMSE de Ey / Ey = ((RMSE de Ty)/N)/(Ty/N) = RMSE de Ty / Ty = RRMSE de Ty
}
# Creee le 1e octobre 2012
######################################################################################################
#' Fonction pour preparerla sortie de strata.geo ou strata cumrootf a partir de la
#' sortie de strata.bh.internal.
#'
#' @param out la sortie de strata.bh.internal
#' Pour la definition des autres parametres, voir le code C
#'
#' @return out la sortie adaptee
#'
out_rule <- function(out, bhfull, nclassh = NULL)
{
if (!is.null(nclassh)) out <- c(list(nclassh = nclassh), out) # Pour cumrootf seulement
out <- c(list(bh = as.vector(bhfull[-c(1, length(bhfull))])), out) # Pour ajouter les bornes (de longueur L-1) a la sortie
out[c("relativebias", "propbiasMSE")] <- NULL # Pour enlever les elements de la sortie en lien avec la strate takenone
names(out)[names(out)=="RMSE"] <- "stderr"
names(out)[names(out)=="RRMSE"] <- "CV"
class(out)<-"strata"
out ## une liste est retournee
}
# Creee le 2 octobre 2012
######################################################################################################
#' Fonction qui permet de d'initialiser quelques statistiques supplementaires utiles aux calculs
#'
init_stat <- function(obj_fct)
{
# Pour tirer de obj_fct les variables dont on a besoin ici
x <- obj_fct$x; N <- obj_fct$N
# Variables relatives a la strate certain
certain <- obj_fct$certain; Nc <- obj_fct$Nc;
# Variables relatives au model
nmodel <- obj_fct$nmodel; beta <- obj_fct$beta; sig2 <- obj_fct$sig2; pcertain <- obj_fct$pcertain;
gamma <- obj_fct$gamma; epsilon <- obj_fct$epsilon;
# Calcul de stat sur les observations
EX <- mean(x) ## Moyenne de toutes les observations (utile seulement pour le modele "random")
EX2 <- sum(x^2)/N ## Moyenne de toutes les observations au carre (utile seulement pour le modele "random")
# Calculs pour obtenir les informations relatives a la certainty stratum
EYc <- if (is.null(certain)) NA else {
get_EYs(xs=x[certain], Ns=Nc, ps=pcertain, obj_fct = as.list(environment()))
}
# Sortie
list(EX=EX, EX2=EX2, EYc=EYc)
}
# Creee le 5 octobre
######################################################################################################
#' Fonction qui permet de changer l'echelle des bornes :
#' elle permet de passer de l'echelle des donnees a l'echelle de la position dans le vecteur x1noc.
#'
#' @param bh un vecteur de bornes de longueur L-1 (b0 et bL non inclus), sur l'echelle des donnees
#' Pour la definition des autres parametres, voir le code C
#'
#' @return pbh vecteur de longueur L-1 representant des bornes de strates, mais sur l'echelle des rangs des donnees :
#' chaque element de pbh est un entier representant la position dans le vecteur x1noc d'une borne.
#'
bh2pbh <- function(bh, x1noc){
pbh <- vector(length=length(bh))
for (i in 1:length(bh)) {
dif <- x1noc - bh[i]
pbh[i] <- sum(dif<0) + 1
}
pbh
}
# Creee le 17 octobre 2012
######################################################################################################
#' Fonction pour determiner les bornes initiales robustes
#'
#' Parametres dont on a besoin : N1noc,Ls,takenone,takeall,minNh,wtx1noc
#'
robust_initpbh <- function(obj_fct)
{
# Pour tirer de obj_fct les variables dont on a besoin ici
N1noc <- obj_fct$N1noc
Ls <- obj_fct$Ls
takenone <- obj_fct$takenone
takeall <- obj_fct$takeall
minNh <- obj_fct$minNh
wtx1noc <- obj_fct$wtx1noc
# Autres initialisation de variables
ipbh <- vector(length=Ls-1)
wtx1noc_copy <- wtx1noc
### strates takeall s'il y en a ##
# On les veux les plus petites possibles, tout en respectant minNh, car elles peuvent causer un n negatif.
# On ne se soucie pas de pouvoir calculer une variance dans ces strates (N1noch>1 ou au moins 2 unites differentes)
# car cette variance peut etre nulle sans causer de problemes dans les calculs subsequents.
if (takeall!=0) {
for (i in 1:takeall) {
pos <- sum(cumsum(rev(wtx1noc_copy))<minNh)
ipbh[Ls-i] <- N1noc-pos
N1noc <- N1noc - pos - 1
wtx1noc_copy <- wtx1noc_copy[1:N1noc]
}
}
### strates takesome ###
# On veut une repartition la plus egale possible en terme des N1noch (nombres d'unites differentes dans les strates),
# ce qui donne des strates avec des plus petites variances que de prendre une repartition la plus egale possible
# en terme des Nh. Si les N1noch ne peuvent etre egaux, on ajoute une unite unique dans les premieres strates
# jusqu'a ce que sum(Nh)=N. Il faut aussi s'assurer que minNh soit respecte. S'il ne l'est pas, les strates avec N1noch
# environ egaux sont modifiees de la facon suivante :
# On veut augmenter Nh dans la ou les strates avec un Nh<checkNh. Pour ce faire, on essaye tous les deplacements possibles
# d'une valeur (representant possiblement plus d'une unite) a partir de n'importe quelle strate vers la strate avec le plus
# petit Nh. Si une ou des solutions verifient checkNh, on arrete la et on conserve celle qui donne la plus petite variance des Nh
# car on les veut aussi les plus egaux possible. Si aucune solution ne verifie checkNh, on reprend la procedure a partir de la
# solution qui donne la plus petite variance des Nh. En cas d'egalite des variances, on favorise des petites strates pour les
# grandes unites.
if (Ls-takeall>1) {
N1nochobj <- N1noc/(Ls-takeall)
N1noch <- rep(floor(N1nochobj),Ls-takeall)
plus <- N1noc-sum(N1noch)
if (plus>0) N1noch[1:plus] <- N1noch[1:plus]+1
pos <- cumsum(N1noch[-(Ls-takeall)])+1
# Il faut maintenant s'assurer que minNh est respecte
Nh <- vector(length=Ls-takeall)
for (i in 1:(Ls-takeall)) Nh[i] <- sum(wtx1noc[ifelse(i>1,pos[i-1],1):ifelse(i<Ls-takeall,pos[i]-1,N1noc)])
checkNh <- all(Nh>=minNh)
iter <- 0
while (!checkNh && iter<2*minNh*Ls) {
change <- matrix(rep(c(pos,NA,NA),Ls-takeall-1),byrow=TRUE,nrow=Ls-takeall-1,ncol=Ls-takeall+1)
# id colonnes de change : pos (longueur Ls-takeall-1) + checkNh + varNh
posmin <- which.min(Nh)
idrow <- 1
for(j in rev((1:(Ls-takeall))[-posmin])) { # ici rev permet de favoriser des petites strates pour les grandes unites
if (posmin<j) {
change[idrow,posmin:(j-1)] <- pos[posmin:(j-1)] + 1
} else {
change[idrow,j:(posmin-1)] <- pos[j:(posmin-1)] - 1
}
for (i in 1:(Ls-takeall)) Nh[i] <- sum(wtx1noc[ifelse(i>1,change[idrow,i-1],1):ifelse(i<Ls-takeall,change[idrow,i]-1,N1noc)])
change[idrow,Ls-takeall] <- all(Nh>=minNh)
change[idrow,Ls-takeall+1] <- var(Nh)
idrow <- idrow + 1
}
checkNh <- sum(change[,Ls-takeall])>0
if (checkNh) change <- change[as.logical(change[,Ls-takeall]),,drop=FALSE]
keeprow <- which.min(change[,Ls-takeall+1])
pos <- change[keeprow,1:(Ls-takeall-1)]
iter <- iter + 1 # pour eviter les boucles infinies imprevues
}
ipbh[1:(Ls-takeall-1)] <- pos
}
### strate takenone s'il y en a ###
# strate takenone initiale nulle car elle est la principale cause d'un n negatif.
if (1==takenone) ipbh <- c(1, ipbh)
### Sortie
ipbh
}
# Creee le 17 octobre 2012 a partir de initbh.robust
######################################################################################################
#' Fonction pour passer des bornes pleines aux bornes de longueur L-1
bhfull2bh <- function(bhfull) { as.vector(bhfull[-c(1, length(bhfull))])}
# Creee le 17 octobre 2012
#' Fonction pour passer des bornes de longueur L-1 aux bornes pleines
bh2bhfull <- function(bh, x) { as.vector(c(min(x), bh, max(x) + 1)) }
# Creee le 17 octobre 2012
#' Fonction pour ajouter une borne pour la strate takenone a des bornes initiales, au besoin
add_b1_takenone <- function(ibhfull, Ls, takenone, x)
{
if ((length(ibhfull)==Ls+1) && (takenone==1)){
# The initial boundary of the take-none stratum is set to the first percentile of x.
b1 <- quantile(x, probs=0.01)
# If this first percentile is equal to the minimum value of x, this initial boundary would
# lead to an empty take-none stratum. In that case, the initial boundary of the take-none
# stratum is rather set to the second smallest value of x
if (b1==min(x)) b1 <- unique(x)[2]
# On insere la nouvelle borne b1 (position 2) et on ordonne le nouveau vecteur car il est
# possible (mais c'est une situation extreme) que le nouveau b1 soit superieur a l'ancien.
sort(c(ibhfull[1], b1, ibhfull[-1]))
} else {
ibhfull # S'il y a deja L - 1 bornes, on n'a rien a faire
}
}
# Creee le 18 octobre 2012
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.