\gdef\mujnazevCS{V\'{y}b\v{e}r podle kv\'{o}t}
\gdef\mujnazevEN{Sample by quota}
\gdef\mujnazevPR{\mujnazevCS}
\gdef\mujnazevDR{\mujnazevEN}
\gdef\mujauthor{Martin Anděl}

\bonushyper
\renewcommand\arraystretch{1.25}

\nazev{\mujnazevCS}

\nazev{\mujnazevEN}

\author{\mujauthor}

\Adresa{STEM/MARK, a.s., Smrčkova 2485/4, 180\,00 Praha~8}

\Email{andel@stemmark.cz}



\Abstrakt{Při realizaci marketingových průzkumů se můžeme setkat s~požadavkem na výběr jednotek z~nějakého seznamu tak, aby tento výběr splňoval předepsané kvóty, které jsou dány jen nekříženě. Tento úkol lze formulovat jako úlohu celočíselného lineárního nebo kvadratického programování a řešit pomocí optimalizačních knihoven v~prostředí R.}
\smallskip
\KlicovaSlova{Kvótní výběr, celočíselné programování}

\Abstract{In market research it is possible to solve the task about selection of  units from the list in such a way, that this sample satisfies given condition on quotas, which are given only marginaly. This task is possible to formulate as the problem of integer linear or quadratic programming and solve with optimisation packages in R.}
\smallskip
\KeyWords{Quota sampling, integer programming}


\section{Úvod}\label{sec01}

 Při výběrových šetřeních se realizují různé druhy výběrových plánů. Může to být obyčejný prostý náhodný výběr nebo některý z~pokročilejších postupů jako jsou vícestupňový nebo stratifikovaný výběr. Chceme-li vybrat vyvážený výběr (balanced sample) a máme k dispozici celou oporu výběru, mohli bychom použít Tillého metodu kostky (cube method), viz \cite{tille}, která je naprogramována např. v~prostředí R v~knihovně {\tt sampling}. Někdy ale může nastat situace, že dostupná opora výběru je jen nějaký částečný seznam, který nemusí ve sledovaných znacích proporčně odpovídat zkoumané populaci. Úlohou je z~tohoto seznamu vybrat výběr o~velikosti $n$, který vyhovuje předepsaným kvótám ve sledovaných znacích. Uvažujme situaci, kdy tyto kvóty jsou  stanoveny jen nekříženě. Chceme například dodržet předepsaný počet mužů a žen ve znaku pohlaví, předepsaný počet v~jednotlivých věkových a vzdělanostních kategoriích. S~touto úlohou se můžeme setkat v~praxi například při doplňování panelu respondentů. Panelisté jsou v~průběhu času postupně obměňováni nebo z~vlastního rozhodnutí z~panelu odejdou. Potřebujeme za tyto panelisty vybrat náhradu, a to v~obdobné struktuře jako je struktura panelistů odešlých.

Představme si, že máme k~dispozici $N=500$ kontaktů a u~každého z~nich údaj o~pohlaví, věkové kategorii a dosaženém %ho 
vzdělání. Z~databáze těchto kontaktů budeme chtít vybrat výběr velikosti $n=100$ v~požadované struktuře (tabulka~\ref{table:dtbstruct}).


\begin{table}[!ht]
\centering
\caption{Databáze a požadovaná struktura výběru}
%\centering
\begin{tabular}{|l|c|c|}
\hline 
\multirow{2}{*}{}  & \multicolumn{2}{c|}{Počty respondentů} \\  
 \cline{2-3}
  & Databáze & \makecell{Požadovaná \\ struktura \\ výběru} \\ 
\hline  
\multicolumn{1}{|l}{Pohlaví} & \multicolumn{2}{c|}{}\\ 
\hline 
\ \ Muž & 156 & \phantom{1}44 \\ 
\hline 
\ \ Žena & 344 & \phantom{1}56 \\ 
\hline 
%\multicolumn{3}{|l|}{}\\ 
\multicolumn{1}{|l}{Věková kategorie \ \ \ \ \ \ \ \ \ \ \ \ \ } & \multicolumn{2}{c|}{}\\ 
\hline 
\ \ \ 15--29 let & 188 & \phantom{10}8 \\ 
\hline 
\ \ \ 30--44 let & 177 & \phantom{1}38 \\ 
\hline 
\ \ \ 45--59 let & \phantom{5}85 & \phantom{1}37 \\ 
\hline 
\ \ \ 60+ let & \phantom{5}50 & \phantom{1}17 \\ 
\hline 
%\multicolumn{3}{|l|}{}\\ 
\multicolumn{1}{|l}{Vzdělání} & \multicolumn{2}{c|}{}\\ 
\hline 
\ \ \ Základní & \phantom{5}32 & \phantom{10}2 \\ 
\hline 
\ \ \ Bez maturity & 166 & \phantom{1}26 \\ 
\hline 
\ \ \ Maturita & 220 & \phantom{1}43 \\ 
\hline 
\ \ \ VŠ & \phantom{5}82 & \phantom{1}29 \\ 
\hline 
\multicolumn{3}{|l|}{}\\ 
\hline 
\ Celkem & 500 & 100 \\ 
\hline 
\end{tabular} 
\label{table:dtbstruct}  % is used to refer this table in the text
\end{table} 


\section{Když výběr existuje}
Označme $\matr{X}$ vstupní datovou matici 0--1 indikátorů sledovaných charakteristik v~databázi. Matice $\matr{X}$ bude mít rozměry $(N \times p)$, kde $N$ je počet kontaktů v~databázi a $p$ je sledovaný počet parametrů. V~našem případě máme $N=500$ a $p=11$. Poslední jedenáctý parametr kontroluje celkový počet vybraných jedinců, v~matici $\matr{X}$ mu odpovídá poslední sloupec obsahující všude hodnoty~$1$. Úlohu výběru se stanovenými kvótami můžeme formulovat jako úlohu binárního lineárního programování. Informace o~optimalizačních úlohách lze nalézt např. v~\cite{optimisation}. Úlohu zapíšeme ve tvaru
\begin{equation}
\begin{gathered}
 \min_{\matr{w}} {\matr{1}}^T {\matr{w}}, \\
 {\matr{X}}^T {\matr{w}} = \matr{s},
\end{gathered}
\end{equation}
kde $\matr{w}$ je hledaný vektor binárních 0--1 proměnných indikujících zahrnutí daného kontaktu do výběru a vektor $\matr{s}$ je vektor s~požadovanou strukturou výběru vyjádřenou absolutními počty kontaktů. Předepsané podmínky zahrnují v~poslední rovnici požadavek na velikost výběru, součet složek vektoru~$\matr{w}$ je tedy pevně dán a úlohu bychom mohli přeformulovat na maximalizaci stejné účelové funkce za shodných vedlejších podmínek. 

Vstupní data úlohy si připravíme v~prostředí R.


\begin{alltt}
make.X <- function()\{
  Xpocet <- matrix(
    c(1,1,1, 6,  1,1,2, 6,  1,1,3,26,  1,1,4,10,  1,2,1, 2,
      1,2,2,25,  1,2,3,17,  1,2,4, 9,  1,3,2,12,  1,3,3,16,
      1,3,4, 3,  1,4,1, 2,  1,4,2, 1,  1,4,3,12,  1,4,4, 9,
      2,1,1,16,  2,1,2,46,  2,1,3,67,  2,1,4,11,  2,2,1, 4,
      2,2,2,44,  2,2,3,50,  2,2,4,26,  2,3,2,20,  2,3,3,25,
      2,3,4, 9,  2,4,1, 2,  2,4,2,12,  2,4,3, 7,  2,4,4, 5),
    ncol=4, byrow=T)
  Xdf <- as.data.frame(Xpocet[rep(1:nrow(Xpocet),Xpocet[,4]),
                              1:3])
  Xdffac <- as.data.frame(lapply(Xdf,as.factor))
  names(Xdffac)<-c('sex','agecat','educat')

X <- with(Xdffac,\{
  cbind(
    model.matrix(~sex-1), 
    model.matrix(~agecat-1),
    model.matrix(~educat-1),
    total = 1)
 \})
\}

X <- make.X()
N <- dim(X)[1]
p <- dim(X)[2]
s <- c(44,56,8,38,37,17,2,26,43,29,100)
\end{alltt}

Funkce {\tt make.X} nejprve definuje matici {\tt Xpocet}, která obsahuje v~prvních třech sloupcích postupně číselné hodnoty kategorií znaků pohlaví, věková skupina a vzdělání, přičemž např. u~znaku pohlaví kód 1 odpovídá kategorii muž, hodnota 2 kategorii žena. Poslední čtvrtý sloupec matice {\tt Xpocet} pak udává četnost výskytu dané kombinace hodnot znaků. Z~této agregované podoby   
funkce vytvoří neagregovaná data a současně kódování jednotlivých znaků převede na 0--1 tvar pomocí příkazu {\tt model.matrix}. Současně vytvoří sloupec samých jedniček, který odpovídá poslednímu parametru pro celkový počet vybraných jedinců. Na posledním řádku kódu jsou do vektoru {\tt s} uloženy požadované počty jednotlivých kategorií ve výběru, přičemž poslední hodnota 100 udává celkový požadovaný počet jedinců ve výběru.    

Úlohu můžeme řešit pomocí knihovny {\tt Rsymphony} (viz \cite{symphony}). 
\begin{alltt}
library('Rsymphony')
result <- Rsymphony_solve_LP(obj=rep(1,N),mat=t(X),
          dir=rep('==',p),rhs=s, types=rep('B',N), max=F)
\end{alltt}

Hledané {\tt w} je uloženo v~objektu {\tt result$\$$solution}. Nemusí jít o jediné možné řešení. Pokud bychom chtěli najít více řešení, mohli bychom použít například knihovnu {\tt Rcplex} (viz \cite{cplex}), která v~jednom ze svých parametrů umožňuje nastavit počet hledaných řešení. Předpokladem ovšem je mít na počítači nainstalovaný optimalizační program CPLEX. Další možnost jak nalézt více řešení bude zmíněna v~kapitole~\ref{roi}. Při hledání dalšího řešení různého od již nalezeného {\tt w} bychom mohli postupovat také tak, že bychom přidali do úlohy další podmínku, která by zajistila různost řešení. Druhé řešení bychom dostali příkazem
\begin{alltt}
result2 <- Rsymphony_solve_LP(obj=rep(1,N),mat=rbind(t(X),w),
           dir=c(rep('==',p),'<='),rhs=c(s,99),
           types=rep('B',N), max=F)
\end{alltt}


\section{Když výběr neexistuje}\label{neexistuje}
Pokud bychom chtěli z~databáze kontaktů uvedené v~kapitole \ref{sec01} udělat výběr o~velikosti $n=250$ s~obdobnou podílovou strukturou v~jednotlivých kategoriích, zjistili bychom, že řešení neexistuje, protože ve věkové skupině 45--59 let máme v~databázi $85$ kontaktů, ale ve výběru bychom chtěli nyní mít 92 kontaktů, viz tabulka~\ref{table:sampstruct}, kde je pro úsporu místa již uveden i výběr, ke kterému dospějeme na konci kapitoly.

\begin{table}[!ht]
\centering
\caption{Výběr blízký požadované struktuře}
%\centering
\begin{tabular}{|l|c|c|c|}
\hline 
\multirow{2}{*}{}  & \multicolumn{3}{c|}{Počty respondentů} \\
\cline{2-4}
  & Výběr & \makecell{Požadovaná \\ struktura \\ výběru} & Databáze \\ 
\hline  
\multicolumn{1}{|l}{Pohlaví} & \multicolumn{3}{c|}{}\\ 
\hline 
\ \ \ Muž & 111 & 111 & 156 \\ 
\hline 
\ \ \ Žena & 139 & 139 & 344\\ 
\hline 
%\multicolumn{4}{|l|}{}\\ 
\multicolumn{1}{|l}{Věková kategorie \ \ \ \ \ \ \ \ \ \ \ \ \ } & \multicolumn{3}{c|}{}\\ 
\hline 
\ \ \ 15--29 let & \phantom{2}21 & \phantom{2}21 & 188\\ 
\hline 
\ \ \ 30--44 let & 102 & \phantom{2}95 & 177 \\ 
\hline 
\ \ \ 45--59 let & \phantom{2}85 & \phantom{2}92 & \phantom{2}85\\ 
\hline 
\ \ \ 60+ let & \phantom{2}42 & \phantom{2}42 & \phantom{2}50\\ 
\hline 
%\multicolumn{4}{|l|}{}\\ 
\multicolumn{1}{|l}{Vzdělání} & \multicolumn{3}{c|}{}\\ 
\hline 
\ \ \ Základní & \phantom{25}5 & \phantom{25}5 & \phantom{2}32\\ 
\hline 
\ \ \ Bez maturity & \phantom{2}66 & \phantom{2}66 & 166 \\ 
\hline 
\ \ \ Maturita & 107 & 107 & 220 \\ 
\hline 
\ \ \ VŠ & \phantom{2}72 & \phantom{2}72 & \phantom{2}82\\ 
\hline 
\multicolumn{4}{|l|}{}\\ 
\hline 
\ Celkem & 250 & 250 & 500\\ 
\hline 
\end{tabular} 
\label{table:sampstruct}  % is used to refer this table in the text
\end{table} 
Musíme úlohu přeformulovat a rovnosti v~podmínkách uvolnit na nerovnosti. Chtěli bychom najít takové nerovnosti, které se od rovností liší co nejméně. Úloha smíšeného binárního a celočíselného lineárního programování bude nyní vypadat takto:
\begin{equation}
\begin{gathered}
\label{eqn:ulohaneexistuje}
%\begin{split}
\min_{(\matr{w},\matr{d})} {\matr{1}}^T {(\matr{w},\matr{d})},\\
\matr{s}-\matr{d} \leq {\matr{X}}^T {\matr{w}} \leq \matr{s}+\matr{d}, \\
 0 \leq d_{i} \leq s_{i},\quad i=1\ldots p-1 \\
 0 = d_{p}
%\end{split}
\end{gathered}
\end{equation}
kde $\matr{w}$ je vektor binárních 0--1 proměnných indikujících zahrnutí daného kontaktu do výběru a $\matr{d}$ je celočíselný vektor o~kolik se může řešení lišit od požadované struktury. Poslední podmínka $0 = d_{p}$ říká, že požadovaný počet jedinců ve výběru $n=250$ má být splněn přesně. Podmínky na vektor odchylek $\matr{d}$ můžeme případně upravit, mohli bychom například některé z nich chtít splnit přesně. Alternativou může být požadavek na splnění odchylky vyjádřené podílem z vektoru $\matr{s}$. Pokud by vektor $\matr{s}$ obsahoval v~některé složce malé číslo, můžeme naopak horní mez zvýšit. Úlohu můžeme opět řešit pomocí R a knihovny {\tt Rsymphony}.
\begin{alltt}
s <- c(111,139,21,95,92,42,5,66,107,72,250)
obj <- c(rep(1,N),rep(1,p))
A <- matrix(0, nrow=2*p, ncol=p+N)
A[1:p,1:N] <- t(X)
A[1:p,(N+1):(N+p)] <- -diag(rep(1,p))
A[(p+1):(2*p),1:N] <- -t(X)
A[(p+1):(2*p),(N+1):(N+p)] <- -diag(rep(1,p))
op <- c(rep('<=',2*p))
rhs <- c(s,-s)
bounds <- list(
  lower = list(ind=c((N+1):(N+p)), val=c(rep(0,p))),
  upper = list(ind=c((N+1):(N+p)), val=c(s[1:p-1],0))  
)
result <- Rsymphony_solve_LP(obj, A, op, rhs,
          types = c(rep('B',N),rep('I',p)),
          max = F, bounds=bounds)
w <- result\$solution[1:N]
tabvysledek <- cbind(t(X)%*%w, s)
colnames(tabvysledek) <- c('vyber','s')
tabvysledek
\end{alltt}
Matice {\tt A} v~kódu slouží k~zápisu nerovností
$$
 \matr{s}-\matr{d} \leq {\matr{X}}^T {\matr{w}} \leq \matr{s}+\matr{d}
$$
ve funkci {\tt Rsymphony\_solve\_LP}. Horní část matice {\tt A} je pro horní nerovnosti $\leq$ a je tvořena dvojicí matic, v~levé části je transponovaná matice {\tt X}, což jsou koeficienty pro hledaný vektor indikující vybrané jednotlivce {\tt w}, v~pravé části je pak se záporným znaménkem jednotková diagonální matice dimenze {\tt p} pro minimalizovaný vektor odchylek $\matr{d}$. Obdobně je pak ještě v~kódu sestavena dolní polovina matice {\tt A} pro dolní nerovnosti.

Výsledný výběr je uveden v~tabulce~\ref{table:sampstruct}. Vzhledem k~nedostatečnému počtu kontaktů v~databázi u~věkové kategorie 45--59 let musíme nějak ve výběru navýšit ostatní věkové kategorie. V~našem případě došlo k~navýšení pouze u~jedné z~věkových kategorií, konkrétně u~kategorie 30--44 let. U~znaků pohlaví a vzdělání se požadovanou strukturu výběru podařilo splnit zcela přesně. Nalezené řešení není jediné možné, k~hledání více řešení bychom mohli opět použít např. program CPLEX nebo postup, který bude uveden na konci kapitoly~\ref{roi}. Také bychom mohli postupně přidávat další podmínky na různost řešení jak bylo naznačeno v předchozí kapitole, nyní ještě s ověřením, zda další řešení dosahuje stejného minima jako již nalezené první řešení.


\section{Kvadratická minimalizace}

I v~této kapitole budeme vycházet ze stejné databáze kontaktů jako v~předchozích úlohách. Navýšení sousední věkové kategorie v~tabulce~\ref{table:sampstruct} by nám mohlo vyhovovat, můžeme ale také preferovat rozložení do více kategorií. Toho dosáhneme minimalizací součtu čtverců odchylek $\sum_{i=1}^{p}d_{i}^{2}$. Pak již ale nevystačíme s~minimalizací lineární funkce, ale půjde o~složitější úlohu smíšeného binárního a celočíselného kvadratického programování. Takovou úlohu můžeme řešit pomocí optimalizačních programů jako např. CPLEX nebo Gurobi. Tyto programy jsou komerční, existuje ale možnost zkušební nebo akademické licence.

Úlohu zformulujeme v~následujícím tvaru.
\begin{equation}
\begin{gathered}
\min_{(\matr{w},\matr{d})} {(\matr{w},\matr{d})}^T \matr{I}_{(N+p,N+p)} {(\matr{w},\matr{d})}, \\
 \matr{s}-\matr{d} \leq {\matr{X}}^T {\matr{w}} \leq \matr{s}+\matr{d}, \\
 0 \leq d_{i} \leq s_{i},\quad i=1\ldots p-1, \\
 0 = d_{p} \\
\end{gathered}
\end{equation}
kde $\matr{w}$ jsou binární 0--1 proměnné a $\matr{d}$ jsou celočíselné proměnné, matice $\matr{I}_{(N+p,N+p)}$ je čtvercová matice s $1$ na diagonále a $0$ mimo diagonálu. Podmínku na velikost výběru chceme splnit přesně. Podmínky na vektor odchylek~$\matr{d}$ můžeme případně upravit, jak již bylo zmíněno u úlohy \ref{eqn:ulohaneexistuje}. Řešení budeme hledat pomocí programu CPLEX, pro který existuje knihovna {\tt Rcplex} (viz \cite{cplex}).
\begin{alltt}
library('Rcplex')
s <- c(111,139,21,95,92,42,5,66,107,72,250)
A <- matrix(0,nrow=2*p,ncol=p+N)
A[1:p,1:N] <- t(X)
A[1:p,(N+1):(N+p)] <- -diag(rep(1,p))
A[(p+1):(2*p),1:N] <- -t(X)
A[(p+1):(2*p),(N+1):(N+p)] <- -diag(rep(1,p))
rhs <- c(s,-s)
result <- Rcplex(cvec = c(rep(0,N+p)), Amat = A, bvec = rhs,
 Qmat = diag(c(rep(1,N),rep(1,p))),lb = 0,
 ub = c(rep(1,N),s[1:(p-1)],0),
 objsense = c('min'), sense = 'L',
 vtype = c(rep('B',N),rep('I',p)))
w <- result\$xopt[1:N]
tabvysledek <- cbind(t(X)\%*\%w,s)
colnames(tabvysledek) <- c('vyber','struktura')
tabvysledek
\end{alltt}

Význam a sestavení matice {\tt A} je zde obdobné jako v~kódu pro úlohu \ref{eqn:ulohaneexistuje}. Řešení pomocí programu Gurobi (viz \cite{gurobi}) by se lišilo voláním optimalizační funkce a načtením knihovny {\tt gurobi}.
\begin{alltt}
library('gurobi')
\ldots
model <- list(Q = diag(c(rep(1,N), rep(1,p))),
            obj = 0, A = A, rhs = rhs,
            sense = c(rep('<=',2*p)),
            lb = c(rep(0,N+p)), ub = c(rep(1,N),s[1:(p-1)],0),
            vtype = c(rep("B",N),rep("I",p)) )
params <- list(OutputFlag = 0)
result <- gurobi(model,params)
w <- result\$x[1:N]
\ldots
\end{alltt}

Dosažené optimální struktury výběru jsou uvedeny v~tabulce~\ref{table:sampoptim}.
Řádek Celkem je proměnná jako ostatní řádky Muž, Žena, \ldots, nejde o~typický součtový řádek, který by byl kontrolním součtem čísel uvedených
výše. Konkrétně v~naší ukázce se podařilo vybrat přesně 250 respondentů, odchylka je~0.

\begin{table}[!htb]
\centering
\caption{Optimální výběr}
%\centering
\begin{tabular}{|l|c|c|c|c|c|c|}
\hline 
\multirow{3}{*}{}  & \multicolumn{6}{c|}{Počty respondentů} \\
\cline{2-7}
 {}  & \multicolumn{2}{c|}{CPLEX} & \multicolumn{2}{c|}{Gurobi} &  
 %\makecell{\\ \\ Požad. \\ strukt. \\ výběru} 
 & 
 %\makecell{\\[2\baselineskip] Data-\\ báze} 
 \\  
 \cline{2-3} \cline{4-5}
  & \makecell{Výběr} &  \makecell{Čt.\\ odch.} 
  & \makecell{Výběr}  &  \makecell{Čt.\\ odch.} 
  & %\raisebox{-0.5\baselineskip}[0.5\baselineskip][0pt]{
  \makecell{Požad.\\ strukt. \\ výběru}
  %}
  & \makecell{Data-\\ báze} 
  \\
\hline  
\multicolumn{1}{|l}{Pohlaví} & \multicolumn{6}{c|}{}\\ 
\hline 
\ \ \ Muž & 111 & \phantom{5}0 & 111 & \phantom{5}0 & 111 & 156\\ 
\hline 
\ \ \ Žena & 139 & \phantom{5}0 & 139 & \phantom{5}0 & 139 & 344\\ 
\hline 
%\multicolumn{7}{|l|}{}\\ 
\multicolumn{1}{|l}{Věková kategorie \ } & \multicolumn{6}{c|}{}\\ 
\hline 
\ \ \ 15--29 let & \phantom{1}24 & \phantom{5}9 & \phantom{5}23 & \phantom{5}4 & \phantom{5}21 & 188 \\ 
\hline 
\ \ \ 30--44 let & \phantom{1}97 & \phantom{5}4 & \phantom{5}98 & \phantom{5}9 & \phantom{5}95 & 177 \\ 
\hline 
\ \ \ 45--59 let & \phantom{1}85 & 49 & \phantom{5}85 & 49 & \phantom{5}92 & \phantom{5}85 \\ 
\hline 
\ \ \ 60+ let & \phantom{1}44 & \phantom{5}4 & \phantom{5}44 & \phantom{5}4 & \phantom{5}42 & \phantom{5}50 \\ 
\hline 
%\multicolumn{7}{|l|}{}\\ 
\multicolumn{1}{|l}{Vzdělání} & \multicolumn{6}{c|}{}\\ 
\hline 
\ \ \ Základní & \phantom{15}5 & \phantom{5}0 & \phantom{55}5 & \phantom{5}0 & \phantom{55}5 & \phantom{5}32 \\ 
\hline 
\ \ \ Bez maturity & \phantom{1}66 & \phantom{5}0 & \phantom{5}66 & \phantom{5}0 & \phantom{5}66 & 166 \\ 
\hline 
\ \ \ Maturita & 107 & \phantom{5}0 & 107 & \phantom{5}0 & 107 & 220 \\ 
\hline 
\ \ \ VŠ & \phantom{1}72 & \phantom{5}0 & \phantom{5}72 & \phantom{5}0 & \phantom{5}72 & \phantom{5}82 \\ 
\hline 
\multicolumn{7}{|l|}{}\\ 
\hline 
\ Celkem & 250 & \phantom{4}0 & 250 & \phantom{4}0 & 250 & 500 \\ 
\hline 
\end{tabular} 
\label{table:sampoptim}  % is used to refer this table in the text
\end{table} 
%\newpage

Nalezená řešení {\tt w} nejsou jediná možná. Například funkci {\tt Rcplex} můžeme v~parametru {\tt n} říci, kolik řešení chceme hledat. Pokud parametr nenastavíme, program hledá jedno řešení. Při nastavení {\tt n = NA} program hledá všechna možná řešení, musíme však počítat s~velkou časovou náročností.  

Dále bychom mohli místo součtu čtverců odchylek minimalizovat veličinu $\chi^2=\sum_{i=1}^{p}d_{i}^{2}/s_{i}$. V~kódu bychom do parametru {\tt Qmat}, resp. {\tt Q}, zadali matici {\tt diag(c(rep(1,N),1/s))}. Výsledek by se lišil jen u~věkových kategorií (viz tabulka~\ref{table:sampchi}). I v~tomto případě existuje více optimálních řešení hledaného vektoru {\tt w}, tato řešení však dají shodné vybrané počty u všech kategorií. 

\begin{table}[!htb]
\centering
\caption{Výběr s~minimálním $\chi^2$}
%\centering
\begin{tabular}{|l|c|c|c|}
\hline 
\multirow{2}{*}{}  & \multicolumn{2}{c|}{Počty respondentů} \\  
\cline{2-3}
  & \makecell{Výběr $\min \chi^2$} & \makecell{Požadovaná \\ struktura \\ výběru} \\ 
\hline  
\multicolumn{1}{|l}{Věková kategorie \ \ \ \ \ \ \ \ \ \ \ \ \ } & \multicolumn{2}{c|}{}\\ 
\hline 
\ \ \ 15--29 let & 22 & 21 \\ 
\hline 
\ \ \ 30--44 let & 99 & 95 \\ 
\hline 
\ \ \ 45--59 let & 85 & 92 \\ 
\hline 
\ \ \ 60+ let & 44 & 42 \\ 
\hline 
\end{tabular} 
\label{table:sampchi}  % is used to refer this table in the text
\end{table} 
%\newpage

%\newpage
\section{Největší možný výběr}
Vidíme, že při $n=250$ můžeme zkonstruovat jen výběr, který je požadované struktuře v~nějakém smyslu nejbližší. Mohlo by nás ještě zajímat, jaký by byl největší možný výběr splňující danou strukturu s~předepsanou přesností. Vyjdeme z~předepsané struktury pro $n=250$, pro ostatní $n$ strukturu proporčně přepočteme, nebudou to tedy již pouze celá čísla. Z~tohoto důvodu nebudeme požadovat nulovou odchylku, ale zkusíme nastavit pro každou z~kategorií maximální možnou odchylku $0{,}9$.

\newpage % temp
Optimalizační úloha bude ve tvaru.
\begin{equation}
\begin{gathered}
\max_{\matr{w}} {\matr{1}}^T {\matr{w}}, \\ 
\matr{-b} \leq {(\matr{X}}^T - \mathrm{diag}{\matr{(s/250})}\matr{1}_{(p,\ N)}){\matr{w}} \leq \matr{b}, \\ 
\end{gathered}
\end{equation}
kde $\matr{w}$ je vektor binárních 0--1 proměnných indikujících zahrnutí daného kontaktu do výběru a $\matr{b}$ je předepsaná odchylka o~kolik se může řešení lišit od požadované struktury.

Řešení budeme hledat pomocí knihovny {\tt Rsymphony}.
\begin{alltt}
zaklad <- s[11]
r <- s/zaklad
obj <- c(rep(1,N)) 
A <- matrix(0, nrow=2*p, ncol=N)
J <- matrix(1, nrow=p, ncol=N)
dr <- diag(r) %*% J
A[1:p,1:N] <- t(X) - dr
A[(p+1):(2*p),1:N] <- t(X) - dr
op <- c(rep('>=',p),rep('<=',p))
b <- c(rep(0.9,p))
rhs <- c(rep(0,p)-b,rep(0,p)+b)
result <- Rsymphony_solve_LP(obj, A, op, rhs,
                             types = c(rep('B',N)),
                             max = T)
w <- result\$solution
tabvysledek <- as.data.frame(cbind(c(t(X)\%*\%w,sum(w)),
  c(round(s/zaklad*sum(w),1),sum(w))))
colnames(tabvysledek) <- c('vyber','zadani')
tabvysledek
\end{alltt}
Výběr je uveden v~tabulce~\ref{table:sampmax}. Našemu zadání tedy vyhovuje výběr o~velikosti $n=233$.

\begin{table}[!ht]
\centering
\caption{Maximální možný výběr}
%\centering
\begin{tabular}{|l|c|c|}
\hline 
\multirow{2}{*}{}  & \multicolumn{2}{c|}{Počty respondentů} \\  
\cline{2-3}
  & Výběr & Zadání \\ 
\hline  
\multicolumn{1}{|l}{Pohlaví} & \multicolumn{2}{c|}{}\\ 
\hline 
\ \ \ Muž & 103 & 103,5 \\ 
\hline 
\ \ \ Žena & 130 & 129,5 \\ 
\hline 
%\multicolumn{3}{|l|}{}\\ 
\multicolumn{1}{|l}{Věková kategorie \ \ \ \ \ \ \ \ \ \ \ \ \ } & \multicolumn{2}{c|}{}\\ 
\hline 
\ \ \ 15--29 let & \phantom{2}19 & \phantom{2}19,6 \\ 
\hline 
\ \ \ 30--44 let & \phantom{2}89 & \phantom{2}88,5 \\ 
\hline 
\ \ \ 45--59 let & \phantom{2}85 & \phantom{2}85,7 \\ 
\hline 
\ \ \ 60+ let & \phantom{2}40 & \phantom{2}39,1 \\ 
\hline 
%\multicolumn{3}{|l|}{}\\ 
\multicolumn{1}{|l}{Vzdělání} & \multicolumn{2}{c|}{}\\ 
\hline 
\ \ \ Základní & \phantom{22}4 & \phantom{22}4,7 \\ 
\hline 
\ \ \ Bez maturity & \phantom{2}61 & \phantom{2}61,5 \\ 
\hline 
\ \ \ Maturita & 100 & \phantom{2}99,7 \\ 
\hline 
\ \ \ VŠ & \phantom{2}68 & \phantom{2}67,1 \\ 
\hline 
\multicolumn{3}{|l|}{}\\ 
\hline 
\ Celkem & 233 & 233,0 \\ 
\hline 
\end{tabular} 
\label{table:sampmax}  % is used to refer this table in the text
\end{table} 


\section{Výběr s použitím knihovny ROI a jejích doplňků}\label{roi}
Pokud nemáme k dispozici na počítači optimalizační programy, můžeme využít k řešení úlohy NEOS server. Přístup k NEOS serveru je zdarma, ale je nutné se zaregistrovat a poskytnout emailovou adresu. Postup si ukážeme na úloze z~kapitoly~\ref{neexistuje} Využijeme k tomu obecné rozhraní {\tt ROI} (viz \cite{roineos}). Níže je pouze upravený kus kódu, do parametru {\tt email} je pak nutno napsat vlastní registrovanou emailovou adresu.     
%\newpage

\begin{alltt}
library(ROI)
library(ROI.plugin.neos)
...
lp_bound <- V_bound(li=c((N+1):(N+p)),ui=c((N+1):(N+p)),
 lb=c(rep(0,p)), ub=round(c(s[1:p-1],0),0)) 
types <- c(rep('B',N),rep('I',p))
lp <- OP(objective = obj,
         L_constraint( L = A,
                       dir = op,
                       rhs = round(rhs,0)),
         bounds=lp_bound,
         types=types,
         maximum = FALSE)
sol <- ROI_solve(lp, solver = "neos", method = "cplex",
                 email = "emailova adresa")
w <- sol\$solution[1:N]
...
\end{alltt}

V příkazu {\tt ROI\_solve} můžeme nastavit i jiné metody. Jedna z možností by mohla být nastavení {\tt method="mosek"}.

Když nyní máme optimalizační úlohu zapsánu ve tvaru pro {\tt ROI}, mohli bychom ještě zkusit hledat více řešení pomocí doplňku {\tt ROI.plugin.msbinlp}. List s více řešeními bychom po načtení knihoven dostali příkazem
\begin{alltt}
library(ROI.plugin.msbinlp)
library(ROI.plugin.symphony)
sol <- ROI_solve(lp, solver = "msbinlp", method = "symphony",
                  nsol_max = 4)
\end{alltt}

Zde jsme do parametru {\tt nsol\textunderscore max} uvedli, že chceme hledat 4 řešení. K~jednotlivým hledaným vektorům {\tt w} pak přistupujeme pomocí indexů, např. první řešení bychom dostali příkazem
\begin{alltt}
w1 <- sol[[1]]\$solution[1:N]
\end{alltt}

Můžeme použít i jiné optimalizační metody, např. {\tt glpk}.
\begin{alltt}
library(ROI.plugin.glpk)
sol <- ROI_solve(lp, solver = "msbinlp", method = "glpk",
                  nsol_max = 4)
\end{alltt}

Nastavení {\tt nsol\textunderscore max=Inf} hledá všechna řešení.




\section{Závěr}
Úlohu vytvoření výběru splňujícího předepsané marginální kvóty můžeme zformulovat jako úlohu binárního či smíšeného binárního a celočíselného lineárního nebo kvadratického programování. K~řešení těchto úloh jsme zvolili prostředí~R a jeho knihovny {\tt Rsymphony}, {\tt RCplex} a {\tt gurobi}. Pro soubory o~velikosti řádu stovek lze řešení nalézt pomocí volně dostupných knihoven, u~větších souborů v~řádu tisíců může být nutné k~řešení zvolit některý z~komerčních optimalizačních programů.

Základní zdrojový kód na vyzkoušení najde čtenář na Githubu:
\smallskip

\noindent\hfil
\url{https://github.com/MartakAnd/Sample-by-quota}


\section*{Poděkování}
Autor děkuje prof. RNDr. Jaromíru Antochovi, CSc., a dvěma anonymním recenzentům za cenné poznámky a návrhy, které text doplnily a učinily ho srozumitelnějším.

%\bigskip

\renewcommand{\refname}{Literatura}
\begin{thebibliography}{9}
\setlength\itemsep{-1pt}

\bibitem{symphony}
    An R interface to the SYMPHONY solver for mixed-integer linear programs. 
    \url{https://cran.r-project.org/web/packages/Rsymphony/}.\\

\bibitem{cplex}
    R Interface to CPLEX. \hfill \break
    \url{https://cran.r-project.org/web/packages/Rcplex/}.

\bibitem{gurobi}
    Gurobi optimisation. 
    \url{https://www.gurobi.com/}.

\bibitem{roineos}
    ROI: R Optimization Infrastructure. \hfill \break
    \url{https://cran.r-project.org/web/packages/ROI/}.

\bibitem{optimisation}
Nemhauser, George L., Wolsey, Laurence A.: \textit{Integer and combinatorial optimization.} Wiley-Interscience series in discrete mathematics and optimization. New York, N.Y.: John Wiley, 1999. ISBN 0-471-35943-2.

\bibitem{tille}
Tillé, Y.: \textit{Sampling Algorithms,} Springer, 2006.

\end{thebibliography}



