\section{Introduction}

The Monte Carlo (MC) method is one of  the widely used numerical methods for simulating probability distributions. However, sometimes it is not possible to sample from a given target distribution.
Markov chain Monte Carlo (MCMC)  methods have been developed to address this problem. Instead of sampling independent points directly, MCMC samples
from a Markov chain whose limiting distribution is the target distribution.
 MCMC has widened the applications of MC in many different fields \cite{Caflisch1998,MeynTweedie1993}. Another deficiency of MC {algorithms} is its slow convergence rate.
Quasi-Monte Carlo (QMC) algorithms  on the other hand {perform} better in improving the convergence rate of Monte Carlo which partially depends on generating  samples with small discrepancy. For a survey of QMC  we refer to \cite{JFS2013,Josef}. Putting the QMC idea into MCMC is a good way to  improve the convergence rate and widen practical applications.  Recently many  results in this direction have been achieved \cite{ChenThesis2011,CDO2011,TribbleThesis2007,TribbleOwen2005}. With this paper we add another result in this direction by using a deterministic driver {sequence} in an acceptance-rejection algorithm. We prove discrepancy bounds of order $N^{-1/s}$, where the dimension of the state space is $s-1$ and $N$ is the number of samples. The discrepancy here is a generalization of the concept of the Kolmogorov-Smirnov test between the empirical distribution of the samples and  the target distribution to higher dimension.
A more detailed description of our results will be provided below.

\subsection{Previous work on MCMC and QMC}
In the following we describe some previous research on  QMC and MC. L'Ecuyer studied the convergence behavior of randomized quasi-Monte Carlo for discrete-time  Markov chains in \cite{Pierre2008}, known as array-RQMC. The general idea is to obtain a better approximation of the target distribution than with the plain   Monte Carlo method by using randomized QMC. In a different direction, Tribble~\cite{TribbleThesis2007} and Tribble and Owen~\cite{TribbleOwen2005}
  established {a} condition under  which low discrepancy sequences can be used for consistent MCMC estimation for finite state spaces.  It has been shown that replacing an IID sequence by a completely uniformly distributed sequence also implies a consistent estimation in finite state spaces.   A construction of  weakly completely uniformly distributed sequences is also proposed in \cite{TribbleOwen2005}. As a sequel to the work of Tribble, Chen in   his thesis \cite{ChenThesis2011} and Chen, Dick and Owen~\cite{CDO2011} demonstrated that Markov chain quasi-Monte Carlo (MCQMC) algorithms using a completely uniformly distributed sequence as
  driver sequence gives a consistent result under certain assumptions on the update function and Markov chain. Further, Chen~\cite{ChenThesis2011} also showed that MCQMC can achieve
    a  convergence rate of $O(N^{-1+\delta})$ for any $\delta>0$ under certain conditions,  but he only showed the existence of a driver sequence.

In our recent work \cite{DRZ2013}, done with  Rudolf,  we prove upper bounds on the discrepancy under the assumptions that the Markov chain is uniformly ergodic
and the driver sequence is deterministic rather than independent uniformly distributed random variables. In particular, we show the existence of driver sequences for which the discrepancy of the Markov chain from the target distribution with respect
to certain test sets converges with (almost) the usual Monte Carlo rate of $N^{-1/2}$. A drawback of this result is that we
 are currently  not able to give an explicit construction of a driver sequence for  which our discrepancy bounds hold for uniformly ergodic Markov chains. Garber and Choppin in \cite{GC2014} adapted  low discrepancy point sets instead of random numbers in sequential Monte Carlo (SMC). They proposed  a new algorithm named sequential quasi-Monte Carlo (SQMC).
  They constructed consistency and stochastic bounds based on randomized QMC point set for this algorithm. It is an open problem to obtain  deterministic bounds for SQMC.
More literature review about applying QMC to MCMC problems can be found in \cite[Section 1]{CDO2011}.

{\subsection{Acceptance-rejection algorithms}}
We now give a description of the algorithms in this paper.
Let $\psi: D\to{\mathbb{R}}_{+}=[0,\infty)$ be our target density function, where $D\subseteq \mathbb{R}^{s-1}$ and $\mathbb{R}_{+}=[0,\infty)$. We consider the cases where $D=[0,1]^{s-1}$ or $D=\mathbb{R}^{s-1}$.
Assume that it is not possible to sample directly from the target distribution.
One possible solution to obtain samples from $\psi$  is to choose a proposal density $H$ from which we can sample and then use an acceptance-rejection algorithm. Assume there exists a constant $L<\infty$ such that $\psi(\boldsymbol{x}) < L H(\boldsymbol{x})$ for all $\boldsymbol{x}\in D$. The  following algorithm can be used to obtain samples with distribution $\psi$.

\begin{algorithm}\label{Basic_AR}(Random acceptance-rejection (RAR) algorithm). Given a target density $\psi:D\to{\mathbb{R}}_{+}$  and a proposal density $H: D\to \mathbb{R}_+$. Assume that there exists a constant $L<\infty$ such that  $\psi(\boldsymbol{x}) < LH(\boldsymbol{x})$ for all $\boldsymbol{x}$ in the domain $D$.
We introduce another random variable $u$ having uniform distribution in the unit interval, i.e. $u\thicksim U([0,1])$. Then the acceptance-rejection algorithm is given by
\begin{enumerate}
  \item Draw $X\thicksim H$ and $u\thicksim U([0,1])$.
  \item Accept $Y=X$ as a sample of $\psi$ if $u\leq \frac{\psi(\boldsymbol{X})}{LH(\boldsymbol{X})}$,
otherwise go back to step 1.
\end{enumerate}
\end{algorithm}
See also \cite{BHL13,Devroye,HLD04} for a discussion of related algorithms. For a discussion on how to select proposal densities see for instance \cite{BHL13} and the references therein.
 The  acceptance-rejection sampler works to sample from an unknown density  based on a proposal density.

Acceptance-rejection sampling and importance sampling \cite[Section~3]{RobertCasella2004} are quite similar ideas. Both of them distort a sample from a  distribution in order to sample from another one.
However, there is a difference in the  selection of the constant $L>\psi(\boldsymbol{x})/H(\boldsymbol{x})$ for  $\boldsymbol{x}$ in the domain $D$.
The acceptance-rejection method does not work when $\sup_{\boldsymbol{x}\in D} \psi(\boldsymbol{x})/H(\boldsymbol{x})=\infty$, while importance
sampling is still available  \cite{OwenbookDraft}. In this paper, we only use the acceptance-rejection sampler to get samples of a given target density,
since we are interested in obtaining discrepancy bounds for MCQMC. More information on general strategies for generating nonuniform random variables can be found in the monographs \cite{Devroye, HLD04}.

 \begin{figure}\label{SobolVSPseudorandom}
 \includegraphics[trim=1.5cm 1cm 1.5cm 1cm, clip=true, totalheight=0.42\textheight, angle=360]{HistogramsSobolVSPseudorandom22.eps}\\
  \caption{Different driver sequences: deterministic  and pseudo-random points (the total number of points in each case is $2^9$) and histograms are w.r.t. samples we accepted.}
\end{figure}

The acceptance-rejection algorithm with deterministic driver sequence is one special class of MCQMC.
From the superior distribution properties in terms of the discrepancy of the Sobol sequence  \cite{Sobol1967} one could expect an improvement in the discrepancy  of the samples obtained from the acceptance-rejection algorithm based on the Sobol sequence.  In one dimension the discrepancy we study is the Kolmogorov-Smirnov test   between the target distribution and the empirical distribution of the sample points.  For a given point set in the $s$-dimensional unit cube, the star discrepancy
  measures the difference between the  proportion of points in a subinterval of $[0,1]^s$ and the Lebesgue  measure of this subinterval. We defer the precise definition
  of discrepancy to Section~\ref{SecBackground}.

In this paper, we replace the  IID initial samples with an explicit construction of the driver sequence by using $(t,m,s)$-nets \cite{Josef,Niederriter1992} (obtained from the Sobol sequence).
Figure \ref{SobolVSPseudorandom} shows a  comparison between different driver sequences: deterministic points (Sobol points) and pseudo random uniform points (they both have $2^9$ points). The acceptance-rejection sampler works by only accepting  those points under the target density curve.  The difference of driver sequences  will affect the samples we obtain by the acceptance-rejection algorithm, hence the distribution properties of the points which were  accepted will be influenced. The right two figures in Fig.\ref{SobolVSPseudorandom} show the histograms of the points which we accepted in both cases. Note that the deterministic  samples  better estimate the density function.
 Our interest in this paper is in entirely deterministic methods. However, one could also use randomized quasi-Monte Carlo point sets \cite{Owen1, Owen2, Owen3} and  study a randomized setting.

\subsection{Previous work on deterministic acceptance-rejection algorithm }
 The deterministic acceptance-rejection  algorithm has also been discussed by Moskowitz and Caflisch~\cite{MC96} and Wang~\cite{Wang1999, Wang00}. Therein a smoothing technique was introduced to improve the numerical performance of the acceptance-rejection algorithm. Wang~\cite{Wang00} gave  a heuristic argument to indicate a convergence rate of order $N^{-\frac{s+2}{2(s+1)}}$. This argument assumes that the points in elementary intervals are uniformly distributed. Thus this reasoning is not fully deterministic. { Our lower bound on discrepancy (Theorem~\ref{lowerPsiN}) indicates that this reasoning does not apply  in our case.}
 The numerical experiments in \cite{Wang00} also indicate an improvement using a well chosen deterministic driver sequence
 (in this case the so-called Halton sequence \cite{Halton}) compared to a random driver sequence. Recently, Nguyen and \"{O}kten in \cite{NO2013} presented a consistency result of an acceptance-rejection algorithm for low-discrepancy sequences. This algorithm yielded good numerical performances on standard deviation and efficiency. However, proving an explicit convergence rate of the discrepancy for this algorithm is still  an open problem. See also \cite{MC95,NO2013} for numerical experiments
 using quasi-Monte Carlo point sets for the related problem of integrating indicator functions.

  It is worth noticing
 that all results given in previous work are empirical evidence and the discrepancy of samples is not directly investigated.
 Our work focuses on discrepancy properties of points produced by totally deterministic acceptance-rejection methods.  We also
 prove discrepancy bounds on deterministic acceptance-rejection algorithms, including an upper bound and a lower bound.  {The combination with the reduced acceptance-rejection sampler provides further evidence of the good
 performance of the deterministic method.} Our algorithm here may also be combined with similar algorithms like the acceptance-complement method, see for instance \cite[Section~II.5]{Devroye}.

Before  presenting the theoretical background, we briefly describe deterministic algorithms and some numerical results which show a convergence rate comparison  using Monte Carlo and quasi-Monte Carlo methods.
\section{Our results}
\subsection{Construction of driver sequence}\label{SubSecConstruction}
In this paper, we use low discrepancy point  sets given by $(t,s)$-sequences (see Definition~\ref{DeftmsNet} and Definition~\ref{DeftsSequence} below) in base $b$ as driver sequences. The first $b^m$ points of a $(t,s)$-sequence are a so called $(t,m,s)$-net in base $b$.
Explicit constructions of $(t,s)$-sequences in base $2$ have been found by Sobol~\cite{Sobol1967}, in prime base $b\ge s$ by Faure~\cite{Faure1982} and in prime-power base $b$ by Niederreiter~\cite{Niederriter1988}.  In all these  constructions $t$ depends only on $s$ but not on $m$. In practice, since digital nets (based on Sobol points) are included in the statistics toolbox of Matlab, this method is very easy to implement. People seeking  more discussion of construction methods can  also consult \cite[Chapters 4\&8]{Josef}.

In the following we describe the algorithm  and present some numerical results.
Since in general it is computationally too expensive to compute the supremum in the definition of the star-discrepancy exactly,  we use a so-called $\delta$-cover to estimate this supremum. An introduction to $\delta$-covers is provided in the appendix.
In the numerical discussion, the driver sequence is  generated by a $(t,m,s)$-net in base $2$. Specifically, we always use a Sobol sequence \cite{Sobol1967} to generate $(t,m,s)$-nets for our experiments.
\subsection{Deterministic algorithm for target densities defined on $[0,1]^s$}\label{SubsecCube}
 We consider now the case where the target density is defined on $[0,1]^{s-1}$. The following algorithm is a deterministic version of Algorithm~\ref{Basic_AR}.  For the proofs later, we need the technical assumption that the target density is pseudo-convex. The definition of pseudo-convexity is discussed in Section~\ref{SecBackground}.
\begin{algorithm}\label{algorithmAR}(Deterministic acceptance-rejection (DAR) algorithm in $[0,1]^s$).  Let the target density $\psi:[0,1]^{s-1}\to \mathbb{R}_+$,  where $s\geq 2$, be pseudo-convex. Assume that there exists a constant $L<\infty$ such that $\psi(\boldsymbol{x}) \leq L$ for  all $\boldsymbol{x}\in[0,1]^{s-1}$. Let $A=\{\boldsymbol{x}\in[0,1]^s:\psi(x_1,\ldots,x_{s-1})\geq L x_s\}$. Suppose we aim to obtain approximately $N$ samples from $\psi$.
\begin{itemize}
  \item [i)~] Let $M=b^m\ge\left \lceil N/(\int_{[0,1]^{s-1}}\psi(\boldsymbol{x})/L d\boldsymbol{x})\right\rceil$, where $m\in\mathbb{N}$ is the smallest integer satisfying this inequality.  Generate a $(t,m,s)$-net $Q_{m,s}=\{\boldsymbol{x}_0,\boldsymbol{x}_1,\ldots,\boldsymbol{x}_{b^m-1}\}$ in base $b$.
  \item [ii)~]Use the acceptance-rejection method for the points $Q_{m,s}$ with respect to the density $\psi$, i.e. we accept the point $\boldsymbol{x}_n$ if $\boldsymbol{x}_n \in A$, otherwise reject. Let $P_N^{(s)}=A\cap Q_{m,s}=\{\boldsymbol{z}_0,\ldots, \boldsymbol{z}_{N-1}\}$ be the sample set we accept.
  \item [iii)~]  Project the points  $P_N^{(s)}$ onto  the first $(s-1)$ coordinates.
  Let $P_N^{(s-1)}=\{\boldsymbol{y}_0,\ldots, \boldsymbol{y}_{N-1}\}\subseteq [0,1]^{s-1}$ be the projections of the points  $P_N^{(s)}$.
 \item [iv)~] Return the point set $P_N^{(s-1)}$.
\end{itemize}
\end{algorithm}

The following example shows a better convergence rate when using a low-discrepancy  driver sequence  rather than random point set. In each example the reported discrepancy for  the AR algorithm using a random diver sequence is the average of $10$ independent runs, which is throughout all the numerical experiments.

\begin{example}\label{Example24Dim}
In this example we consider a non-product target density in $[0,1]^4$.
Let target density  $\psi$ be
\begin{equation*}
    \psi(x_1,x_2,x_3,x_4)=\frac{1}{4}(e^{-x_1}+e^{-x_2}+e^{-x_3}+e^{-x_4}), ~~ (x_1,x_2,x_3,x_4) \in[0,1]^4.
\end{equation*}
\end{example}
Figure~\ref{fig:Example24Dim} shows the discrepancy by using deterministic points and pseudo-random points as driver sequence.
For the RAR algorithm, we observe  a convergence rate of $N^{-0.482}$, whereas the DAR algorithm shows a convergence rate of the discrepancy of order $N^{-0.659}$.

 \begin{figure}[!htp]\begin{center}

  \includegraphics[trim=0.5cm 0.25cm 0.5cm 0.5cm, clip=true, totalheight=0.4\textheight, angle=360]{P120140311.eps}\\
  \caption{Convergence order of the star discrepancy of Example~\ref{Example24Dim}.}\label{fig:Example24Dim} \end{center}
\end{figure}
\subsection{Deterministic algorithm for target density defined in real state space}\label{SubsecR}
Now we extend the domain of the target density $\psi$ to  $\mathbb{R}^{s-1}$ with $s\ge 2$.
 Assume that there is  a proposal density function ${H}: \mathbb{R}^{s-1}\to \mathbb{R}_+$ such  that there exists a  constant ${L}<\infty$ such that $\psi(\boldsymbol{z})\le {L} {H}(\boldsymbol{z})$ holds for all $\boldsymbol{z}\in \mathbb{R}^{s-1}$.

 The inverse  Rosenblatt transformation  is used to generate samples from the proposal density
in the real state space $\mathbb{R}^{s-1}$.
 Let ${F}$ be the joint CDF of ${H}$ and ${F}_j(z_j|z_1,\ldots,z_{j-1})$ be the conditional CDF of the proposal density for $j = 1,\ldots,s-1$.
The  transformation ${T}$ is used to
generate points in $\mathbb{R}^{s-1}\times \mathbb{R}_+$ from the unit cube $[0,1]^{s}$, such that the projection of points onto the first $s-1$ coordinates has distribution $H$. More precisely,
let $T: [0,1]^{s}\to \mathbb{R}^s$ be the transformation  given by
\begin{equation}\label{InverseRosenblattTrans}\begin{array}{ll}
&\boldsymbol{z}=T(\boldsymbol{u})\\[0.25cm]
&=\left\{\begin{array}{lll}
  &z_1=F_1^{-1}(u_1), \\[0.25cm]
  &z_j=F_j^{-1}(u_j|u_1,\ldots,u_{j-1}), ~2\leq j\leq s-1,\\[0.25cm]
  &z_{s}= u_{s}H(z_1,\ldots,z_{s-1}).\\[0.25cm]
  \end{array}\right.\end{array}
\end{equation}
The first $s-1$ coordinates  are produced by  the inverse  Rosenblatt transformation which converts the points from the unit cube $[0,1]^{s-1}$ into $\mathbb{R}^{s-1}$.  The $s$th coordinate is uniformly distributed on the line
 $$\{(1-v)(z_{1},\ldots z_{s-1},0)+v(z_{1},\ldots z_{(s-1)},H(z_1,\ldots,z_{s-1}), 0\le v\le1\}$$ if
 $u_s$ is uniformly distributed in $[0,1]$. More details with respect to the Rosenblatt transformation and extensions can be found in \cite{Chentsov1967,Nohu2007,Rosenblatt1952}.

\begin{algorithm}\label{algorithmRd} (Deterministic acceptance-rejection algorithm in $\mathbb{R}^s$). Let an unnormalized target density function $\psi:\mathbb{R}^{s-1}\to \mathbb{R}_+$, where $s\ge 2$, be given. Let ${H}$ be a  proposal density $H:\mathbb{R}^{s-1}\to \mathbb{R}_+$, such that there exists a constant $L<\infty$ such that
$\psi(\boldsymbol{z})\leq L H(\boldsymbol{z})$ for all $\boldsymbol{z}\in \mathbb{R}^{s-1}$.
Let $A=\{\boldsymbol{z}\in\mathbb{R}^s:\psi(z_{1},\ldots,z_{s-1})\ge L H(z_{1},\ldots,z_{s-1}) z_s\}$. {Suppose we aim to obtain approximately $N$ samples from $\psi$.}
\begin{itemize}
                   \item [i)~]  Let $M= b^m\ge\left\lceil N\int_{[0,1]^{s-1}}H(\boldsymbol{x})d\boldsymbol{x}/(\int_{[0,1]^{s-1}}\psi(\boldsymbol{x})/L d\boldsymbol{x})\right\rceil$, where $m\in \mathbb{N}$ is the smallest integer satisfying this inequality. Generate  a $(t,m,s)$-net $Q_{m,s}=\{\boldsymbol{x}_{0},\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{M-1}\}$ in base $b$.
                    \item [ii)~]  Transform the  points into $\mathbb{R}^{s-1}\times \mathbb{R}_+$ from $[0,1]^s$ using the  transformation $T$ given in \eqref{InverseRosenblattTrans} to obtain $\{T(\boldsymbol{x}_0),\ldots,T(\boldsymbol{x}_{M-1})\}$.
                    \item [iii)~] Take the acceptance-rejection method for the sample $T(\boldsymbol{x}_{n})$ with respect to $H$ and $\psi$ in $\mathbb{R}^{s-1}\times\mathbb{R}_+$, i.e. accept the point $T(\boldsymbol{x}_{n})$ if $T(\boldsymbol{x}_{n})\in A$, otherwise reject.
                    Let $P_N^{(s)}=A\cap T(Q_{m,s})=\{\boldsymbol{z}_0,\ldots, \boldsymbol{z}_{N-1} \}$.
                    \item [iv)~]  Project the points $P_N^{(s)}$ we accepted onto  the first $(s-1)$-dimensional space. Denote the first $s-1$ coordinates of the points we accept by the acceptance-rejection method by $P_N^{(s-1)}=\{\boldsymbol{y}_0,\ldots, \boldsymbol{y}_{N-1}\}\subseteq \mathbb{R}^{s-1}$.
                     \item [v)~] Return the point set $P_N^{(s-1)}$.
\end{itemize}
\end{algorithm}
 We provide an example to demonstrate the performance of Algorithm~\ref{algorithmRd}.
\begin{example}\label{Example3} Let the target density function  be given by
\begin{equation*}
  \psi(x_1,x_2)=\left\{\begin{array}{ll}
\displaystyle\frac{4}{\pi}e^{-(x_1+x_2)}(x_1x_2)^{1/2},& \quad x_1,x_2>0,\\
    0,& \quad \mbox{ else.}
    \end{array}\right.
\end{equation*}
\end{example}
The proposal density function  $H$,  which we  use to do the  acceptance-rejection to generate samples of $\psi(x_1,x_2)$,  is chosen as
 \begin{equation*}
   H(x_1,x_2)=\left\{ \begin{array}{ll}
  \displaystyle\frac{1}{4}, &0\leq x_1,x_2\leq 1,\\[0.2cm]
  \displaystyle\frac{1}{4x_2^2}, &0\leq x_1\leq 1,x_2 >1,\\[0.2cm]
    \displaystyle\frac{1}{4x_1^2}, &x_1 >1,0\leq x_2\leq 1,\\ [0.2cm]
  \displaystyle\frac{1}{4x_1^2x_2^2}, &x_1,x_2 >1,\\[0.2cm]
  0, &\quad \mbox{else.}\\
      \end{array}\right.
 \end{equation*}
For this choice of $H$, we use transform $T$ defined in Equation~\eqref{InverseRosenblattTrans} to obtain samples from $H$.
The sample $(x_{j,1},x_{j,2})$ is given by the following transformation
\begin{equation*}\begin{array}{ll}
x_{j,1}=&\left\{\begin{array}{lll}
  2u_{j,1}, &0\le u_{j,1}\le 1/2,\\
   1/2(1-u_{j,1}),& 1/2<u_{j,1}\le 1,\\
    \end{array}\right.\end{array}
\end{equation*}
\begin{equation*}\begin{array}{ll}
x_{j,2}=&\left\{\begin{array}{lll}
  2u_{j,2}, &0\le u_{j,2}\le 1/2,\\
   1/2(1-u_{j,2}),& 1/2<u_{j,2}\le 1.\\
    \end{array}\right.\end{array}
\end{equation*}
 Note that $(u_{j,1},u_{j,2})$ is the driver sequence given by a $(t,m,2)$-net in base $b$.

The order of the star discrepancy is demonstrated in  Figure~\ref{OrderGamma3-2} where $N$ is the number of {accepted} samples. The numerical experiments  show that the star discrepancy converges  at a rate of $N^{-0.720}$ for this example using quasi-Monte Carlo samples as proposal.  The RAR algorithm converges with order $N^{-0.390}$. Again, the DAR sampler outperforms the RAR sampler.

\begin{figure}[!htp]\begin{center}

 \includegraphics[trim=0.5cm 0.25cm 0.5cm 0.5cm, clip=true, totalheight=0.4\textheight, angle=360]{P2Gamma20140310.eps}\\
 \caption{ Convergence order of the star discrepancy of Example~\ref{Example3}.}\label{OrderGamma3-2}\end{center}
\end{figure}

\subsection{A deterministic reduced acceptance-rejection sampler}\label{SubsecSum}
In this subsection we consider  an extension of the DAR sampler.  The random version of this reduced method was recently introduced by Barekat and Caflisch in \cite{BC2013}.  For a target density function $\psi$, we carefully select $H$ such that for $\psi-H$ and $H$ the inverse CDF can be computed. For the case $\psi(\boldsymbol{x})> H(\boldsymbol{x})$, we write $\psi=(\psi-H)+H$ and get samples according to $\psi-H$ and $H$  respectively.

 Figure~\ref{FishProblem} illustrates this method. The sample sets of $\psi$ can be divided into three subsets, $R_{1,1}, R_{1,2}$ and $R_{2,2}$, where $R_{1,2}$ and $R_{2,2}$ can be directly generated by using the inverse CDF of $H$ and $\psi-H$ in a certain range. The acceptance-rejection method is only used to obtain $R_{1,1}$.  Compared with the ordinary acceptance-rejection sampler, one obvious merit of this method is that we do not require $\psi(\boldsymbol{x})\le H(\boldsymbol{x})$ in the whole domain. Also, this method might give better convergence rates {since $R_{1,2}$ and $R_{2,2}$ are obtained via inversion and therefore have low discrepancy}.  Algorithm~\ref{Alg-RAR} gives a simple version of the improved method.  More discussion of a general version is available in Section~\ref{SecImprovedAlg}.

 \begin{figure}[!htp]\begin{center}

 \includegraphics[trim=0.5cm 0.25cm 0.5cm 0.5cm, clip=true, totalheight=0.4\textheight, angle=360]{FishE3April23.eps}\\
 \caption{The solid line represents the target density $\psi$ and the dashed line is the proposal density $H$. }\label{FishProblem}\end{center}
\end{figure}

\begin{algorithm}\label{Alg-RAR}{(Deterministic reduced acceptance-rejection (DRAR) algorithm)} Let $\psi:[0,1]\to \mathbb{R}_+$ be a target density.  Choose a proposal density $H$ such that $\psi-H$ and $H$  can be sampled directly.
Let\begin{align*}
    &\mathcal{S}:=\{x\in [0,1]: \psi(x)< H(x)\}\nonumber \\
    \mbox{ and }&\\
    &\mathcal{L}:=\{x\in [0,1]: \psi(x)\ge H(x)\}.\nonumber \\
\end{align*} Assume that $\int_{\mathcal{S}}\psi(x)dx/\int_{[0,1]}\psi(x)dx,\int_{\mathcal{L}}H(x)dx/\int_{[0,1]}\psi(x)dx$ and $\int_{\mathcal{L}}(\psi-H)(x)dx/\int_{[0,1]}\psi(x)dx$ can be calculated or estimated.
Let $F_{H,\mathcal{S}}^{-1},F_{H,\mathcal{L}}^{-1}$ be the inverse  CDF of the proposal density $H$ in the domain $\mathcal{S}$ and $\mathcal{L}$ respectively and $F_{\psi-H,\mathcal{L}}^{-1}$ be the inverse CDF with respect to $\psi-H$ in $\mathcal{L}$.
Suppose we aim to generate approximately $N$ samples from $\psi$. Let
$$N_1=\left\lceil N\frac{ \int_{\mathcal{S}}\psi(x)dx}{\int_{[0,1]}\psi(x)dx}\right\rceil,  N_2=\left\lceil N\frac{\int_{\mathcal{L}}H(x)dx}{\int_{[0,1]}\psi(x)dx}\right\rceil
\mbox{ and } N_3=\left\lceil N\frac{\int_{\mathcal{L}}(\psi-H)(x)dx}{\int_{[0,1]}\psi(x)dx}\right\rceil.$$
{\begin{itemize}
 \item [i)~] Let $\{\boldsymbol{x}_0,\boldsymbol{x}_1, \boldsymbol{x}_2,\ldots\}\subset [0,1]^2$ be a $(t,2)$-sequence in base $b$.

 \item [ii)~]   Use the acceptance-rejection method with the target density $\psi$ and the proposal density $H$ on the domain $\mathcal{S}$ using  $\{\boldsymbol{x}_0,\boldsymbol{x}_1, \ldots,\boldsymbol{x}_{M-1}\}$ as driver sequence. Choose $M$ such that $N_1$ points are accepted by the DAR algorithm. Compute $F^{-1}_{H,\mathcal{S}}(\boldsymbol{x}_n)$ for $n=0,1,2,\ldots, M-1$. Let $\boldsymbol{z}_0,\boldsymbol{z}_1, \ldots,\boldsymbol{z}_{N_1-1}$ be the accepted points. Label the point set as $R_{1,1}$.
 \item [iii)~]  Compute the points $F_{H,\mathcal{L}}^{-1}(\boldsymbol{x}_n)$ for $n=0,1,\ldots, N_2-1$. Let $R_{1,2}=\{F_{H,\mathcal{L}}^{-1}(\boldsymbol{x}_n):0\le n< N_2\}$.
  \item [iv)~]  Compute the points $F_{\psi-H,\mathcal{L}}^{-1}(\boldsymbol{x}_n)$ for $n=0,1,\ldots, N_3-1$. Let $R_{2,2}=\{F_{\psi-H,\mathcal{L}}^{-1}(\boldsymbol{x}_n):0\le n< N_3\}$.
 \item [v)~]  Project the points in $R_N=R_{1,1}\cup R_{1,2}\cup R_{2,2}$ onto  the first  coordinate. Return the point set $R_N^{(1)}$.
    \end{itemize}}
\end{algorithm}
 Since the inverse transform is a measure-preserving transformation, it can preserve the uniformities of the driver sequence. Thus $R_{1,2}$ and $R_{2,2}$ are low discrepancy point sets. The following example verifies the efficiency  of the DRAR algorithm. A theoretical result about the discrepancy properties of samples obtained by this class of algorithms is provided in  Theorem~\ref{triangleDisc}.

\begin{example}\label{ExampleSumCase} Let $\psi(x)=\sin(4x)+x^2$  be a density function defined on $[0,1]$. Instead of seeking  a proposal density $H$ such that  $\psi(x)\le H(x)$, we notice that inversion  can be implied to $\sin(4x)$ and $x^2$ independently. However, it can not work for their  sum. Choose $H(x)=x^2$. We only do deterministic acceptance-rejection with respect to the target density $\psi$ and proposal density $H$ in the subinterval $\mathcal{S}=(\pi/4,1]$.  In the remaining range $\mathcal{L}=[0,\pi/4]$,  we apply the inverse transformation on $H$ and $\psi-H$ to obtain samples based on a deterministic driver sequence.\end{example}
The discrepancy of the point set generated by Algorithm~\ref{Alg-RAR}  converges at the rate of  $N^{-0.929}$, which is significantly better than the $N^{-0.501}$ convergence rate of a random driver sequence, see Figure~\ref{figDiscrepSumcase}.

\begin{figure}
\includegraphics[trim=0.5cm 0.25cm 0.5cm 0.5cm, clip=true, totalheight=0.4\textheight, angle=360]{FishProblem20140423.eps}\\
   \caption{Convergence order of the star discrepancy of Example~\ref{ExampleSumCase}.}\label{figDiscrepSumcase}
\end{figure}

\section{Background on  discrepancy theory and $(t,m,s)$-nets}\label{SecBackground}
In this section we  first establish some notation and  some useful definitions and then obtain theoretical results.
First we introduce the definition of $(t,m,s)$-nets in base $b$ (see  \cite{Josef}) which we use as the driver sequence throughout the paper. The following fundamental definitions of elementary interval and fair sets are used to define a $(t,m,s)$-net and $(t,s)$-sequence in base $b$.
\begin{definition} ({b-adic elementary interval}). Let $b\geq 2$ be an integer. An $s$-dimensional $b$-adic elementary interval is an interval of the form
\begin{equation*}
    \prod_{i=1}^s\left[\displaystyle\frac{a_i}{b^{d_i}},\displaystyle\frac{a_i+1}{b^{d_i}}\right)
\end{equation*}
with integers $0\leq a_i<b^{d_i}$ and $d_i\geq 0$ for all $1\leq i\leq s$. If $d_1,\ldots,d_s$ are such that $d_1+\cdots+d_s=k$, then we say that the elementary interval is of order $k$
\end{definition}

 \begin{definition}({fair sets}).~For a given set $P_N=\{\boldsymbol{x}_0,\boldsymbol{x}_1,\ldots,\boldsymbol{x}_{N-1}\}$ consisting of $N$ points in $[0,1)^s$, we say for a subset $J$ of $[0,1)^s$ to be fair with respect to $P_N$, if
\begin{equation*}
   \displaystyle\frac{1}{N}\sum_{n=0}^{N-1}1_J(\boldsymbol{x}_n)=\lambda(J),
\end{equation*}
where $1_J(\boldsymbol{x}_n)$ is the indicator function of the set $J$ and $\lambda$ is the Lebesgue measure. \end{definition}

\begin{definition}\label{DeftmsNet} ({$(t,m,s)$-nets in base b}).~For a given dimension $s\geq 1$, an integer base $b\geq 2$, a positive integer $m$ and an integer $t$ with $0\leq t\leq m$,
a point set $Q_{m,s}$ of $b^m$ points in $[0,1)^s$ is called a  $(t,m,s)$-nets in base $b$ if the point set $Q_{m,s}$ is fair with respect to all b-adic s-dimensional elementary intervals of order at most $m-t$.\end{definition}

\begin{definition}\label{DeftsSequence}($(t,s)$-sequence).~For a given dimension $s\geq 1$, an integer base $b\geq 2$ and a positive integer $t$, a sequence $\{\boldsymbol{x}_0, \boldsymbol{x}_1,\ldots\}$ of points in $[0,1)^s$ is called a  $(t,s)$-sequence in base $b$ if for all integers $m\ge t$ and $k\ge 0$, the point set consisting of the points $\boldsymbol{x}_{kb^m},\ldots,\boldsymbol{x}_{kb^m+b^m-1}$ forms a $(t,m,s)$-net in base $b$.
\end{definition}
The concept of discrepancy is introduced in \cite{Niederriter1974} to measure the deviation of
 a sequence from the uniform distribution.  Now we give the definition of the so-called
 star discrepancy which enables us to distinguish  the quality of point sets with respect to the uniform distribution.
\begin{definition}\label{stardiscrepancy} ({star discrepancy}).~Let  $P_N=\{\boldsymbol{x}_0, \boldsymbol{x}_1,\ldots,\boldsymbol{x}_{N-1}\}$ be a point set in  $[0,1)^s$. The star discrepancy $D_N^*$ is defined by
\begin{equation*}
    D_N^*(P_N)=\sup_{J\subset [0,1)^s}\left| \displaystyle\frac{1}{N}\sum_{n=0}^{N-1}1_J(\boldsymbol{x}_n)-\lambda(J)\right|,
\end{equation*}
where the supremum is taken over all  $J=\prod_{i=1}^s [0,\beta_i)\subseteq [0,1)^s$.\end{definition}
See Figure~\ref{Fig:StarDiscrepancy} for an illustration of the concept of discrepancy in the unit square.
\begin{figure}[!htp]\begin{center}

 \includegraphics[trim=1cm 0.25cm 0.5cm 0.5cm, clip=true, totalheight=0.4\textheight, angle=360]{StarDiscrepancy.eps}\\
 \caption{The discrepancy measures the difference between the proportion of points in each rectangle $J$ and the Lebesgue measure of $J$. The star discrepancy
  is defined by the supremum of discrepancies over all rectangles $J$. If we change $J$ to be  convex sets, we obtain the so-called isotropic discrepancy.}\label{Fig:StarDiscrepancy}
 \end{center}
\end{figure}

If we extend the supremum in Definition \ref{stardiscrepancy} over all convex sets in $[0,1]^s$, we get another interesting discrepancy, the so-called isotropic discrepancy.
 It is  another measure of the distribution properties of point sets with respect to convex sets.
\begin{definition}  ({isotropic discrepancy}). \label{isotropic_dis} ~Let  $P_N=\{\boldsymbol{x}_0,\boldsymbol{x}_1,\ldots,\boldsymbol{x}_{N-1}\}$  be a point set in  $[0,1)^s$. The isotropic discrepancy $J_N$ is defined to be
\begin{equation*}
    J_N(P_N)=\sup_{J\subset\mathcal{C}}\left| \displaystyle\frac{1}{N}\sum_{n=0}^{N-1}1_J(\boldsymbol{x}_n)-\lambda(J)\right|,
\end{equation*}
where $\mathcal{C}$ is the family of all convex subsets of $[0,1)^s$.\end{definition}

For  further reading about the definition and properties of discrepancy, we refer for instance to \cite{Josef,Niederriter1974}.

For our purposes here we need the definition of pseudo-convex sets which we introduce in the following (see also \cite[Definition 2]{CJJ2011} and Figure~\ref{FigPseudoconvex} for an example).

\begin{definition}\label{pseudoConvexSet}(pseudo-convex set). Let $A$ be an open subset of $[0,1]^s$ such that there exists a collection of $p$ convex  subsets $A_1,\ldots,A_p$ of $[0,1]^s$
satisfying\begin{enumerate}
            \item $A_i\cap A_j=\emptyset$ for $i\neq j$;
            \item $A\subseteq ( A_1\cup\cdots \cup A_p)$ of $[0,1]^s$ ;
            \item either $A_j$ is a convex part of $A$ ($A_j\subseteq A$ for $j=1,\ldots,q$) or the complement of $A$ with respect to $A_j$, $A_j'=A_j\backslash A$ is convex.
          \end{enumerate}
Then $A$ is called a pseudo-convex set and $A_1,\ldots, A_p$ is an admissible convex covering for $A$ with $p$ parts and with $q$ convex parts of $A$.\end{definition}

\begin{figure}\label{FigPseudoconvex}
\includegraphics[trim=0.5cm 0.25cm 0.5cm 0.5cm, clip=true, totalheight=0.4\textheight, angle=360]{Pseudoconvex2.eps}\\
   \caption{Shows a  pseudo-convex set in the unit square given by the area below graph of the density function and its admissible convex covering.

  Let A be given by the graph under the curve. Then $A_i, i=1,\ldots,5$, is an admissible convex covering of $A$, where $A_1, A_2$ and $A_3$ are convex parts in $A$ but $A_4$ and $A_5$ are rectangles shadowed region covering the remaining part of $A$. The regions $A_4\backslash A$ and $A_5\backslash A$ are convex.}
\end{figure}

\begin{remark}\label{remark-Pseudoconvex}For convenience, we call a nonnegative function pseudo-convex if and only if  the region below its graph is a pseudo-convex  set.\end{remark}

 Next we present a bound on the isotropic discrepancy of points generated by  $(t,m,s)$-nets. A detailed proof is given in Appendix~\ref{AppendixUpperBound}.

\begin{lemma}\label{IsoDisBound} Let the point set $Q_{m,s}=\{\boldsymbol{x}_0,\boldsymbol{x}_1,\ldots, \boldsymbol{x}_{M-1}\}\subseteq[0,1]^s$ be a  $(t,m,s)$-net in base $b$ where $M=b^m$.
 For the  isotropic discrepancy of $Q_{m,s}$  we have
\begin{equation*}
J_M(Q_{m,s})\leq 2sb^{t/s}M^{-1/s}.
\end{equation*}
\end{lemma}

A slightly weaker result than Lemma~\ref{IsoDisBound} can also be obtained from \cite[Korollar~3]{NW75}.
\begin{lemma}\label{StarIsotropic} For any point set $P_N$ in $[0,1]^s$ we have
\begin{equation*}
J_N(P_N) \le 2 s \left(\frac{4s}{s-1}\right)^{(s-1)/s} ( D_N^*(P_N) )^{1/s}.
\end{equation*}
\end{lemma}
Further it is known from \cite{Kr06} that the star discrepancy of a $(t,m,s)$-net $Q_{m,s}$ in base $b$, where $M = b^m$, satisfies
\begin{equation*}
D_{M}^*(Q_{m,s}) \le M^{-1} b^{t} (\log M)^{s-1} \frac{b^s}{(b+1) 2^s (s-1)! (\log b)^{s-1}} + C_s M^{-1} b^{t} (\log M)^{s-2},
\end{equation*}
for some constant $C_s > 0$. These two inequalities therefore yield a convergence rate of order $M^{-1/s} (\log M)^{1-1/s}$.

The following lemma will be used to get a discrepancy bound  for a point set on a pseudo-convex set. It is an extension of \cite[Lemma 5]{CJJ2011} to the unit cube.
\begin{lemma} \label{lemmaConvex}
Let $A$ be a pseudo-convex subset of $[0,1]^s$ with admissible convex covering of $p$ parts with $q$ convex parts of $A$. Then for any point set $P_N=\{\boldsymbol{x}_0,\boldsymbol{x}_1,\ldots,\boldsymbol{x}_{N-1}\}\subseteq [0,1]^s$ we have
\begin{equation*}
    \left|\displaystyle\frac{1}{N}\sum_{n=0}^{N-1}1_A(\boldsymbol{x}_n)-\lambda(A)\right|\leq (2p-q) J_N(P_N).
\end{equation*}
\end{lemma}

\section{Discrepancy investigation of deterministic acceptance-rejection sampler}\label{TheoricalResults}

The first result we get is a discrepancy bound with respect to the target density of samples generated by the acceptance-rejection algorithm with deterministic driver sequences.
The star discrepancy of points generated by the acceptance-rejection algorithm with respect to the target density  converges  at the rate of $N^{-1/s}$, where $N$ is the number of accepted samples. See Theorem~\ref{Maintheorem} for details.
The proof  uses a bound on the discrepancy of our driver sequence with respect to convex sets (which is called isotropic discrepancy, see Definition~\ref{isotropic_dis} for details).

\subsection{Upper bound }\label{UpboundCube}
Let an unnormalized  density function $\psi:[0,1]^{s-1}\to \mathbb{R}_+$ be pseudo-convex, and $\int_{[0,1]^{s-1}}\psi (\boldsymbol{z})d\boldsymbol{z} >0$, but not necessarily $1$.
 Assume that there exists a constant $L<\infty$ such that $\psi(\boldsymbol{x}) \leq L$ for  all $\boldsymbol{x}\in[0,1]^{s-1}$. Let the subset under the graph of $\psi/L$  be defined as
\begin{equation}\label{SubsetA}
  A=\{\boldsymbol{x}\in[0,1]^s:\psi(x_1,\ldots,x_{s-1})\geq L x_s\},
\end{equation}
 which is pseudo-convex in $[0,1]^s$ as $\psi$ is a  pseudo-convex function. Assume that there is an admissible convex covering of $A$ with $p$ parts and with $q$ convex parts of $A$.
  Without loss of generality, let $A_1,\ldots,A_q$ be the convex subsets of $A$ and $A_{q+1},\ldots,A_p$, such that  $A_j'=A_j\backslash A$ is convex for  $q+1\leq j\leq p$.

 The definition of the star discrepancy of a point set $\{\boldsymbol{y}_0,\boldsymbol{y}_1,\ldots, \boldsymbol{y}_{N-1}\}$ with respect to a
 density function $\psi$ is given as follows.

\begin{definition} Let $\psi:[0,1]^{s-1}\to \mathbb{R}_+$ be an unnormalized  target density. Let $\{\boldsymbol{y}_0,\boldsymbol{y}_1,\ldots, \boldsymbol{y}_{N-1}\}$ be a
 point set in $[0,1]^{s-1}$. The star discrepancy of $\{\boldsymbol{y}_0,\boldsymbol{y}_1,\ldots, \boldsymbol{y}_{N-1}\}$ with respect to the  density $\psi$ is defined by
\begin{equation*}\label{DefDiscrepancyWRTdensity}
   D_{N,\psi}^*(P_N^{(s-1)})=\sup_{\boldsymbol{t}\in[0,1]^{s-1}} \left|\displaystyle\frac{1}{N}\sum_{n=0}^{N-1}1_{[\boldsymbol{0},\boldsymbol{t})}({\boldsymbol{y}_n})
    -\displaystyle\frac{1}{C}\int_{[\boldsymbol{0},\boldsymbol{t})}\psi({\boldsymbol{z}})d\boldsymbol{z}\right|,
\end{equation*}
where $C=\int_{[0,1]^{s-1}}\psi(\boldsymbol{z})d\boldsymbol{z}$ and $[\boldsymbol{0},\boldsymbol{t})=\prod_{j=1}^{s-1}[0,t_j)$.
\end{definition}
\begin{remark}
Note that $\frac{1}{C}\psi$ is a probability density function on $[0,1]^{s-1}$. Thus  the discrepancy  in Definition \ref{DefDiscrepancyWRTdensity}
measures the difference between the distribution $\frac{1}{C}\psi$ and the empirical distribution of the sample points with respect to the test
 sets $[\boldsymbol{0}, \boldsymbol{t})$ for $\boldsymbol{t}\in [0,1]^{s-1}$.
\end{remark}

\begin{theorem} \label{Maintheorem}
Let the unnormalized  density function $\psi:[0,1]^{s-1}\to \mathbb{R}_+$, with $s\geq 2$, be pseudo-convex.  Assume that there is an admissible convex covering of
 $A$ given by Equation~\eqref{SubsetA}
 with $p$ parts and with $q$ convex parts of $A$.
Then the  discrepancy of the point set $\{\boldsymbol{y}_0,\boldsymbol{y}_2,\ldots, \boldsymbol{y}_{N-1}\} \subseteq [0,1]^{s-1}$ generated by
Algorithm \ref{algorithmAR} using a $(t,s)$-sequence in base $b$, for large enough $N$, satisfies
\begin{equation*}
   D_{N,\psi}^*(P_N^{(s-1)})
    \le 8C^{-1}L s b^{t/s}(2p-q) N^{-1/s},
  \end{equation*}
  where $C=\int_{[0,1]^{s-1}}\psi({\boldsymbol{z}})d\boldsymbol{z}$ and $\psi (\boldsymbol{x})\le L$ for all $\boldsymbol{x}\in[0,1]^{s-1}$.
\end{theorem}

We postpone the proof of this theorem to Appendix~\ref{AppendixUpperBound}.

\subsection{Lower bound}\label{Lowerbound}
In this section, we provide a lower bound on the star discrepancy with respect to a  convex density function. The general idea is to find, for a given driver point set, a density function satisfying a certain convergence rate.
{
\begin{theorem}\label{lowerPsiN} Let $P_M$ be an arbitrary  point set   in $[0,1]^s$. Then there exists a concave density function $\psi$ defined in $[0,1]^{s-1}$ such that, for $N$ samples generated by the acceptance-rejection algorithm with respect to $P_M$ and $\psi$, we have
\begin{equation*}
 D_{N,\psi}^*(P_N)\ge c_sN^{-\frac{2}{s+1}},
\end{equation*}
where $c_s >0$ is independent of $N$ and $P_M$ but depends on $s$.
\end{theorem}}

A detailed proof is provided in Appendix~\ref{AppendixLowerbound}. We would like to point out that
the lower bound also limits the convergence rates which we can obtain in our current approach via convex sets.
Note that a concave function is also pseudo-convex as defined in Remark~\ref{remark-Pseudoconvex}.

Additionally,  note that \cite{Beck} (in dimension $s=2$) and \cite{Stu77} (for dimension $s > 2$)
showed the existence of points with discrepancy with respect to convex sets bounded from above by $N^{-2/(s+1)} (\log N)^{c(s)}$ (where $c(s)$ is a function of only $s$). This would yield an improvement of our results from $N^{-1/s}$ to $N^{-2/(s+1)}(\log N)^{c(s)}$, however, those constructions are not explicit and can therefore not be used in computation.

\subsection{Generalization to real state space}\label{BoundsGeneralization}
We consider now the case where the target density is defined on $\mathbb{R}^{s-1}$ with $s\ge 2$.  The aim is to show a discrepancy bound on samples generated by the deterministic acceptance-rejection method.
The discrepancy  with respect to a given density function $\psi:\mathbb{R}^{s-1}\to \mathbb{R}_+$ is defined as follows.
\begin{definition}
Let $P_N=\{\boldsymbol{z}_0,\boldsymbol{z}_1,\ldots,\boldsymbol{z}_{N-1}\}$ be a point set in $\mathbb{R}^{s-1}$. Let $\psi:\mathbb{R}^{s-1}\to \mathbb{R}_+ $ be an unnormalized probability density function. Then the star discrepancy
 $D_{N,\psi}^*(P_N)$ is defined by
\begin{equation*}
    D_{N,\psi}^*(P_N)=\sup_{\boldsymbol{t}\in \mathbb{R}^{s-1}}\left|\frac{1}{N}\sum_{n=0}^{N-1}1_{(-\boldsymbol{\infty},\boldsymbol{t}]}(\boldsymbol{z}_n)-\frac{1}{C}\int_{(-\boldsymbol{\infty},\boldsymbol{t}]} \psi(\boldsymbol{z})d\boldsymbol{z}\right|,
   \end{equation*}
where $C=\int_{\mathbb{R}^{s-1}} \psi(\boldsymbol{z})d\boldsymbol{z}$ and $(-\boldsymbol{\infty},\boldsymbol{t}]=\prod_{j=1}^{s-1}(-\infty, t_j]$ for $\boldsymbol{t}=(t_1, \ldots,t_{s-1})$.
\end{definition}

We use the transformation $T$ given in Equation~\eqref{InverseRosenblattTrans} to generate samples of ${H}$. For the sake of investigating discrepancy, the following result is helpful. The lemma shows that the transformations $T$ and its inversion $T^{-1}$  are both measure-preserving. For the proofs later, we assume that the proposal density $H$ is a product measure, i.e. $H=\prod_{j=1}^{s-1}H_j$, where $H_j$ is the marginal density with respect to $z_j$. In our numerical examples, the proposal density is not necessarily of product type.

\begin{lemma} The transformation $T$ from the $s$-dimensional unit cube to $\mathbb{R}^{s-1}\times \mathbb{R}_+$ given in (\ref{InverseRosenblattTrans})
 is measure-preserving, i.e. $\mbox{Volume}(T(D))= \mbox{Volume}(D)$ holds for any measurable set $D\subseteq [0,1]^s$. This is true for $T^{-1}$ as well.
\end{lemma}

To prove a bound on the discrepancy of the samples generated by Algorithm~\ref{algorithmRd},
the following assumption is needed.
\begin{assumption}\label{AssumptionGeneral} Let $\psi$ be the target density and ${H}$ be a product measure proposal density function, which is chosen such that its inverse CDF can be computed.
Let $A=\{\boldsymbol{z}\in\mathbb{R}^s:\psi(z_{1},\ldots,z_{s-1})\ge Lz_s H(z_{1},\ldots,z_{s-1})\}$ and the transformation $T^{-1}$  is defined as the inversion of transform $T$. Then we assume that $T^{-1}(A)$ is pseudo-convex.
\end{assumption}

 As the mappings $T$ and $T^{-1}$ are measure preserving, and since there are the same number of samples in an arbitrary subset $D\subseteq [0,1]^{s}$ and the corresponding subset $T(D)\subseteq \mathbb{R}^{s-1}\times \mathbb{R}_+$, we can consider the discrepancy in the unit cube instead of that in $\mathbb{R}^{s-1}\times \mathbb{R}_+$. Following by  similar proof arguments as for Theorem \ref{Maintheorem} and Theorem~\ref{lowerPsiN}, we obtain the same discrepancy bounds including an upper bound and a lower bound for the general density $\psi$ defined in the real state space $\mathbb{R}^{s-1}$.

\begin{theorem}\label{GneralMaintheorem}  Let the unnormalized target density  $\psi:\mathbb{R}^{s-1}\to \mathbb{R}_+$ and the proposal density $H:\mathbb{R}^{s-1}\to\mathbb{R}_+$ satisfy Assumption \ref{AssumptionGeneral}. Then the  discrepancy of the point set $P_N^{(s-1)}=\{\boldsymbol{y}_0, \boldsymbol{y}_1,\ldots, \boldsymbol{y}_{N-1}\} \subseteq \mathbb{R}^{s-1}$ generated by Algorithm \ref{algorithmRd} satisfies
\begin{equation*}
D_{N,\psi}^*(P_N^{(s-1)})
    \le 8LC^{-1} s b^{t/s}(2p-q) N^{-1/s},
\end{equation*}
for $N$ large enough, where $C=\int_{\mathbb{R}^{s-1}} \psi(\boldsymbol{z})d\boldsymbol{z}$ and $L$ is such that $\psi(\boldsymbol{x})\le LH(\boldsymbol{x})$ for all $\boldsymbol{x}\in \mathbb{R}^{s-1}$.
\end{theorem}

\begin{theorem}\label{lowerRealState} Let $H$ be a product density function defined on $\mathbb{R}^{s-1}$. Let $T$ be the transformation given in Equation~\eqref{InverseRosenblattTrans} associated  to $H$. Let $P_M$ be an arbitrary  point set in $[0,1]^s$, then $T(P_M)$ is a point set in $\mathbb{R}^{s-1}$. Then there exists an unnormalization density function $\psi$ defined in $\mathbb{R}^{s-1}$ satisfying the assumption in Theorem~\ref{GneralMaintheorem} such that  the star discrepancy of the points generated by the acceptance-rejection sampler with respect to $\psi$ and $H$ satisfies
\begin{equation*}
  D_{N,\psi}^*(P_N^{(s-1)})\ge c_sN^{-\frac{2}{s+1}},
\end{equation*}
where $c_s$ is independent of $N$ and $P_M$, but only dependent on $s$.
\end{theorem}

\section{{Discrepancy properties of the deterministic reduced acceptance-rejection sampler}}\label{SecImprovedAlg}
Algorithm~\ref{Alg-RAR} can be extended to a  more general case.
Consider the target density $\psi(\boldsymbol{x})=\sum_{\ell=1}^{k} H_\ell(\boldsymbol{x}), \boldsymbol{x}\in D\subset \mathbb{R}^s$. If it is possible to  sample from $H_{\ell}(\boldsymbol{x})$ individually and the expectations of $H_{\ell}$  can be calculated or estimated with low cost, then we can use an embedding deterministic reduced acceptance-rejection sampler in each step.
{Let
\begin{align}\label{SjLj}
    &\mathcal{S}_{\ell}=\{\boldsymbol{x}\in D: \psi_{k-{\ell}+1}(\boldsymbol{x})< H_{\ell}(\boldsymbol{x})\},\nonumber \\
    \mbox{and}&\nonumber \\
    &\mathcal{L}_{\ell}=\{\boldsymbol{x}\in D: \psi_{k-{\ell}+1}(\boldsymbol{x})\ge H_{\ell}(\boldsymbol{x})\},\nonumber \\
\end{align} where $\psi_{k-{\ell}+1}(\boldsymbol{x})=\sum_{i={\ell}}^{k}H_i(\boldsymbol{x})$ for ${\ell}=1,\ldots,k-1$, and, in particular, $\psi_k$ is the target density.

Suppose we aim to sample $N$ points from the target density $\psi$. The sample set can be divided into two types, namely, points generated from the sets $\mathcal{S}_{\ell}$'s and $\mathcal{L}_{\ell}$'s respectively.
We apply a deterministic acceptance-rejection method given in Algorithm~\ref{algorithmRd} in each $\mathcal{S}_{\ell}$ with respect to $\psi_{k-{\ell}+1}$ and $H_{\ell}$. Note that we get  $\lceil N\int_{\mathcal{S}_{\ell}}\psi_{k-{\ell}+1}(\boldsymbol{x})d\boldsymbol{x}/\int_D\psi_k (\boldsymbol{x}) d\boldsymbol{x}\rceil$ points from $\mathcal{S}_{\ell}$ for ${\ell}=1,\ldots,k-1$.
For sampling from $\mathcal{L}_{\ell}$, the remaining samples come from applying the inverse transformation of $H_{\ell}$ in $\mathcal{L}_{\ell}$. Then we obtain additional   $\lceil N\int_{\mathcal{L}_{\ell}}H_{\ell}(\boldsymbol{x}) d\boldsymbol{x}/\int_D\psi_k(\boldsymbol{x}) d\boldsymbol{x}\rceil$ points from $\mathcal{L}_{\ell}$ for ${\ell}=1,\ldots,k$. We conduct the  procedure inductively until we get samples from $H_k(\boldsymbol{x})$. We assume that $\int_{\mathcal{S}_{\ell}}\psi_{k-{\ell}+1}(\boldsymbol{x})d\boldsymbol{x}/\int_D\psi_k (\boldsymbol{x}) d\boldsymbol{x}$ and $\int_{\mathcal{L}_{\ell}}H_{\ell}(\boldsymbol{x}) d\boldsymbol{x}/\int_D\psi_k(\boldsymbol{x}) d\boldsymbol{x}$ can be calculated or estimated.

The following algorithm is an extension of the  DRAR algorithm, which summarizes the embedding idea.
\begin{algorithm}\label{Alg-RAR-2} Let $\psi(\boldsymbol{x})=\sum_{{\ell}=1}^{k} H_{\ell}(\boldsymbol{x}), \boldsymbol{x}\in D\subset \mathbb{R}^{s-1}$, be a target density we aim to sample from. Define $\psi_{k-{\ell}+1}(\boldsymbol{x})=\sum_{i={\ell}}^{k}H_i(\boldsymbol{x})$ for $j=1,\ldots,k-1$. Denote $\mathcal{S}_{\ell}$ and $\mathcal{L}_{\ell}$ like in Equation~\eqref{SjLj} and assume that $\int_{\mathcal{S}_{\ell}}\psi_{k-{\ell}+1}(\boldsymbol{x})d\boldsymbol{x}/\int_D\psi (\boldsymbol{x}) d\boldsymbol{x}$ and $\int_{\mathcal{L}_{\ell}}H_{\ell}(\boldsymbol{x})d\boldsymbol{x}/\int_D\psi (\boldsymbol{x}) d\boldsymbol{x}$ can be calculated or estimated.
Further assume that we can sample from $H_{\ell}$ individually by applying the  transformation $T_{H_{\ell},\mathcal{S}_{\ell}}$ and $T_{H_{\ell},\mathcal{L}_{\ell}}$ given in Equation~\eqref{InverseRosenblattTrans} in $\mathcal{S}_{\ell}$ and $\mathcal{L}_{\ell}$ respectively. Suppose we aim to generate $N$ samples from $\psi$.  Let  $$N_{1,{\ell}}=\left\lceil N\frac{\int_{\mathcal{S}_{\ell}}\psi_{k-{\ell}+1}(\boldsymbol{x})d\boldsymbol{x}}{\int_D\psi (\boldsymbol{x}) d\boldsymbol{x} }\right\rceil \mbox{ and } N_{2,j}=\left\lceil N\frac{\int_{\mathcal{L}_{\ell}}H_{\ell}(\boldsymbol{x})d\boldsymbol{x}}{\int_D\psi (\boldsymbol{x}) d\boldsymbol{x}} \right\rceil.$$
For  ${\ell}$ from 1 to $k$ do:
\begin{itemize}
 \item [i)~]Let $\{\boldsymbol{x}_0,\boldsymbol{x}_1, \boldsymbol{x}_2,\ldots\}\subset [0,1]^s$ be a $(t,s)$-sequence in base $b$.

 \item [ii)~] Compute $T_{H_{\ell},\mathcal{S}_{\ell}}(\boldsymbol{x}_n)$ for $n=0,1,2,\ldots$. Use the acceptance-rejection method with respect to  $\psi_{k-{\ell}+1}$ and  $H_{\ell}$ on the domain $\mathcal{S}_{\ell}$ using
 $\{\boldsymbol{x}_0,\boldsymbol{x}_1, \ldots,\boldsymbol{x}_{M-1}\}$ as driver sequence. Choose $M$ such that $N_{1,{\ell}}$ points are accepted by the DAR algorithm. Let $R_{1,{\ell}}=\{\boldsymbol{z}_0,\boldsymbol{z}_1, \ldots,\boldsymbol{z}_{N_{1,{\ell}}-1}\}$ be the accepted points.

 \item [iii)~]  Compute $T_{H_{\ell},\mathcal{L}_{\ell}}(\boldsymbol{x}_n)$ for $ n=0,1,\ldots,N_{2,{\ell}}-1$. Let $R_{2,{\ell}}=\{T_{H_{\ell},\mathcal{L}_{\ell}}(\boldsymbol{x}_n): 0\le n<N_{2,{\ell}}\}$.

 \end{itemize}
 Let $R_N^{(s)}=\bigcup_{{\ell}=1}^{k} (R_{1,{\ell}}\cup R_{2,{\ell}})$ and let $R_N^{(s-1)}$ denote the projection of $R_N^{(s)}$ onto the first $s-1$ coordinates. Return the set $R_N^{(s-1)}$.
\end{algorithm}
Now we consider the discrepancy properties of sample points produced by this algorithm.  Note that the sample set of $\psi=\sum_{{\ell}=1}^{k}H_{\ell}$ can be decomposed into several subsets with different star discrepancy.
\begin{theorem}\label{triangleDisc}
For a given target density $\psi(\boldsymbol{x})=\sum_{{\ell}=1}^{k} H_{\ell}(\boldsymbol{x}), \boldsymbol{x}\in D\subset \mathbb{R}^{s-1}$, let $\psi_{k-{\ell}+1}(\boldsymbol{x})=\sum_{i={\ell}}^{k}H_i(\boldsymbol{x})$.  Let $\mathcal{S}_{\ell}$ and $\mathcal{L}_{\ell}$ be given by \eqref{SjLj}.
Let $R^{(s-1)}_N$ be the sample set generated by Algorithm~\ref{Alg-RAR-2}, where $$N_{1,{\ell}}=\left\lceil N\int_{\mathcal{S}_{\ell}}\psi_{k-{\ell}+1}(\boldsymbol{x})d\boldsymbol{x}/\int_D\psi_k (\boldsymbol{x}) d\boldsymbol{x}\right\rceil,$$ which is the number of points generated from $\mathcal{S}_{\ell}$, and
$$N_{2,{\ell}}=\left\lceil N\int_{\mathcal{L}_{\ell}}H_{\ell}(\boldsymbol{x}) d\boldsymbol{x}/\int_D\psi_k(\boldsymbol{x}) d\boldsymbol{x}\right\rceil,$$ which is the number of points generated from $\mathcal{L}_{\ell}$ for ${\ell}=1,\ldots,k$. Assume that $N_{1,{\ell}}$ and $N_{2,{\ell}}$ can be calculated or estimated for the  given target density $\psi$ and $N$. Then we have
\begin{equation*}
    D_{N,\psi}^*(R^{(s-1)}_N)\le \sum_{{\ell}=1}^{k-1}\frac{N_{1,{\ell}}}{N}D^*_{\mathcal{S}_{\ell},\psi_{k-{\ell}+1}} +\sum_{{\ell}=1}^{k}\frac{N_{2,{\ell}}}{N}D^*_{\mathcal{L}_{\ell},H_{\ell}}+\frac{1}{N},
\end{equation*} where $D^*_{\mathcal{S}_{\ell},\psi_{k-{\ell}+1}}$ and $D^*_{\mathcal{L}_{\ell},H_{\ell}}$ is  the discrepancy of the samples in $\mathcal{S}_{\ell}$ and $\mathcal{L}_{\ell}$ respectively.
\end{theorem}
The proof of Theorem~\ref{triangleDisc} is given in Appendix~\ref{SumCaseDiscrepancy}. Note that this  method  achieves an improved acceptance rate of points since we are only  rejecting points in a certain range. For the remaining domain, we  get samples by applying the inverse transform. To be more exact, all point sets from $\mathcal{L}_{\ell}$ have low discrepancy since the inverse transformation is directly applied with respect to $H_{\ell}$ for ${\ell}=1,\ldots,k$. Now we consider the star discrepancy of points generated from $\mathcal{S}_{\ell}$.

The following result from \cite{Kr06} gives an improved upper bound on the star discrepancy on the first $M$ terms of a $(t,s)$-sequence in  base $b$ with $s\ge 2$.
\begin{lemma}\label{ImprovedMsequence}The star discrepancy of the first $M$ terms of a $(t,s)$-sequence $P_M$ in  base $b$ with $s\ge 2$ satisfies
\begin{equation*}
D_{M}^*(P_M) \le M^{-1} b^{t} (\log M)^{s} \frac{b^s(b-1)}{(b+1) 2^{s+1} (s)! (\log b)^{s}} + C_s M^{-1} b^{t} (\log M)^{s-1},
\end{equation*}
for some constant $C_s > 0$ only depending on $s$.
\end{lemma}
With the help of Lemma~\ref{StarIsotropic}, we obtain a bound on the isotropic discrepancy of the first $M$ points of a $(t,s)$-sequence.
\begin{lemma}Let the point set $P_M=\{\boldsymbol{x}_0,\boldsymbol{x}_1,\ldots,\boldsymbol{x}_{M-1}\}\subset[0,1]^s$ be the first $M$ terms of a $(t,s)$-sequence. For the isotropic discrepancy of $P_{M}$ we have
\begin{equation*}
  J_M(P_M)\le 2 s \left(\frac{4s}{s-1}\right)^{(s-1)/s} \Big( \frac{b^{t/s}  \big(\frac{b^s(b-1)}{(b+1) 2^{s+1} (s)! (\log b)^{s}} \big)^{1/s}\log M}{M^{1/s} } + \frac{C'_s  b^{t/s} (\log M)^{(s-1)/s}}{M^{1/s}} \Big),
\end{equation*} for some constant $C'_s>0$ depending only on $s$.
\end{lemma}

Hence, for the star discrepancy of $R_{1,{\ell}}$ for $1\le {\ell}< k$, using a $(t,s)$-sequence as a diver sequence in the DAR algorithm we have  a convergence rate of order $N_{1,{\ell}}^{-1/s}\log N_{1,{\ell}}$. We omit  a detailed proof since similar arguments as for proving Theorem~\ref{Maintheorem} can be used. The following corollary holds by substituting the proper upper bounds and $N_{1,{\ell}}, N_{2,{\ell}}$ in terms of $N$.
\begin{corollary} Suppose that  the target density $\psi(\boldsymbol{x})=\sum_{{\ell}=1}^{k} H_{\ell}(\boldsymbol{x}), \boldsymbol{x}\in D\subset \mathbb{R}^{s-1}$ satisfies all assumptions stated in  Theorem~\ref{triangleDisc}. Let $R^{(s-1)}_N$ be the sample set generated by Algorithm~\ref{Alg-RAR-2}. Then  we have
\begin{equation*}
   D_{N,\psi}^*(R^{(s-1)}_N)\le \sum_{j=1}^{k-1}\frac{C_{\mathcal{S}_{\ell},\psi_{k-{\ell}+1}}\alpha_j^{1-1/s}\log(\alpha_{\ell}N) }{N^{1/s}} +\sum_{{\ell}=1}^{k}\frac{C_{\mathcal{L}_{\ell},H_{\ell}}(\log\beta_{\ell}N)^{s-1}}{N}+\frac{1}{N},
\end{equation*} where
 $$\alpha_{\ell}=\frac{\int_{\mathcal{S}_{\ell}}\psi_{k-{\ell}+1}(\boldsymbol{x}) d\boldsymbol{x}}{\int_{D}\psi(\boldsymbol{x}) d\boldsymbol{x}}
  \quad \mbox{and} \quad \beta_{\ell}=\frac{\int_{\mathcal{L}_{\ell}}H_{\ell}(\boldsymbol{x}) d\boldsymbol{x}}{\int_{D}\psi(\boldsymbol{x}) d\boldsymbol{x}},$$
  and $C_{\mathcal{S}_{\ell},\psi_{k-{\ell}+1}}$ and $C_{\mathcal{L}_{\ell},H_{\ell}}$ are constants associated with $\mathcal{S}_{\ell},\psi_{k-{\ell}+1}$ and $\mathcal{L}_{\ell},H_{\ell}$ respectively.
\end{corollary}

\section{Conclusion and outlook}\label{Conclusion}
As is well known, the integration error using a Monte Carlo method converges at the rate of $N^{-1/2}$. The acceptance-rejection sampler with
 a deterministic driver sequence, which is a simple class of MCQMC methods,  performs much better in our numerical
 experiments than the theoretical result $N^{-1/s}$ of Theorem~\ref{Maintheorem} or Theorem~\ref{GneralMaintheorem},
  for a density defined on $[0,1]^{s-1}$ or $\mathbb{R}^{s-1}$, would imply. The three examples even demonstrate that it is
   possible to achieve a better convergence rate than  with standard  Monte Carlo using a well-chosen deterministic driver sequence.

The main drawback of the acceptance-rejection sampler is that it might reject many points in high dimension.
Some methods to improve the  acceptance rate of points are included in \cite{HLD04}. For a special class of density functions given by a finite sum, we propose an embedding deterministic reduced acceptance-rejection algorithm. This algorithm  produces a better star discrepancy convergence rate in our numerical example.


\section{{{Introduction}}}

The idea to use noncommutative coordinates is quite old - it goes
back to Heisenberg and was firstly formalized by Snyder in
\cite{snyder}. Recently, however,  there were  found new formal
arguments based mainly on Quantum Gravity \cite{2w}, \cite{2a} and
String Theory models \cite{recent}, \cite{string1}, indicating that
space-time at Planck scale  should be noncommutative, i.e. it should
have a quantum nature.

Presently, it is well known, that in accordance with the
Hopf-algebraic classification of all deformations of relativistic
and nonrelativistic symmetries, one can distinguish three
types of quantum spaces \cite{class1}, \cite{class2} (for details see also \cite{nnh}):\\
\\
{ \bf i)} Canonical ($\theta^{\mu\nu}$-deformed) type of quantum space \cite{oeckl}-\cite{dasz1}
\begin{equation}
[\;\hat{ x}_{\mu},\hat{ x}_{\nu}\;] = i\theta_{\mu\nu}\;, \label{noncomm}
\end{equation}
\\
{ \bf ii)} Lie-algebraic modification of classical space-time \cite{dasz1}-\cite{lie1}
\begin{equation}
[\;\hat{ x}_{\mu},\hat{ x}_{\nu}\;] = i\theta_{\mu\nu}^{\rho}\hat{ x}_{\rho}\;,
\label{noncomm1}
\end{equation}
and\\
\\
{ \bf iii)} Quadratic deformation of Minkowski and Galilei  spaces \cite{dasz1}, \cite{lie1}-\cite{paolo}
\begin{equation}
[\;\hat{ x}_{\mu},\hat{ x}_{\nu}\;] = i\theta_{\mu\nu}^{\rho\tau}\hat{
x}_{\rho}\hat{ x}_{\tau}\;, \label{noncomm2}
\end{equation}
with coefficients $\theta_{\mu\nu}$, $\theta_{\mu\nu}^{\rho}$ and  $\theta_{\mu\nu}^{\rho\tau}$ being constants.\\
\\
Besides, it has been demonstrated in \cite{nnh}, that in the case of
so-called N-enlarged Newton-Hooke Hopf algebras
$\,{\mathcal U}^{(N)}_0({ NH}_{\pm})$ the twist deformation
provides the new  space-time noncommutativity of the
form\footnote{$x_0 = ct$.},\footnote{ The discussed space-times have been  defined as the quantum
representation spaces, so-called Hopf modules (see e.g. \cite{oeckl}, \cite{chi}), for quantum N-enlarged
Newton-Hooke Hopf algebras.}
\begin{equation}
{ \bf 4)}\;\;\;\;\;\;\;\;\;[\;t,{ x}_{i}\;] = 0\;\;\;,\;\;\; [\;{ x}_{i},{ x}_{j}\;] =
if_{\pm}\left(\frac{t}{\tau}\right)\theta_{ij}(x
\;, \label{nhspace}
\end{equation}
with time-dependent  functions
$$f_+\left(\frac{t}{\tau}\right) =
f\left(\sinh\left(\frac{t}{\tau}\right),\cosh\left(\frac{t}{\tau}\right)\right)\;\;\;,\;\;\;
f_-\left(\frac{t}{\tau}\right) =
f\left(\sin\left(\frac{t}{\tau}\right),\cos\left(\frac{t}{\tau}\right)\right)\;,$$
$\theta_{ij}(x) \sim \theta_{ij} = {\rm const}$ or
$\theta_{ij}(x) \sim \theta_{ij}^{k}x_k$ and  $\tau$ denoting the time scale parameter
 -  the cosmological constant. Moreover,  the  different relations  between all mentioned above quantum spaces ({\bf 1)}, { \bf 2)}, { \bf 3)}
and { \bf 4)}) have been summarized in paper \cite{nnh}.

It should be noted that the  described  above classification can be supplemented by the proper deformations of so-called $l$-conformal Galilei Hopf structure \;${\cal U}^{(l)}_0(\mathcal{G})$ provided in \cite{czar}. In general, a conformal extension of the Galilei Hopf algebra is parametrized by a positive
half integer $l$, which justifies the term $l$-conformal Galilei Hopf structure. Particularly, the instance of
$l=\frac{1}{2}$ (well-known in the literature as the Schroedinger algebra) has been the focus of most studies (for a review see e.g. \cite{przeglad}).
Besides, motivated by current investigation of the nonrelativistic version of the AdS/CFT correspondence, the interest in conformal
Galilei Hopf algebras with $l>\frac{1}{2}$ is growing rapidly in the last time  \cite{pocz}-\cite{kon}.

In this article  we investigate  the Abelian twist deformations of \;${\cal U}^{(l)}_0(\mathcal{G})$ Hopf structure
which provide (as we shall see for a moment)  six types of space-time noncommutativity (see formulas (\ref{rspacetime1})-(\ref{rspacetime6})).
 The motivations for such a kind of  studies are manyfold.  First of all, we construct explicitly the deformation of  quite general (known)
Hopf algebra at nonrelativistic level.
Secondly, we get the completely new quantum (twist-deformed) space-times associated with \;${\cal U}^{(l)}_{\alpha}(\mathcal{G})$ Hopf
structure. Such a result seems to be quite interesting due to the fact, that it extends in natural way the mentioned above classification
of the quantum spaces. Finally, it should be noted  that the obtained results permit to consider the classical as well as quantum particle models defined on new  noncommutative nonrelativistic space-times (\ref{rspacetime1})-(\ref{rspacetime6}).

The paper is organized as follows. In first section we recall  the basic facts concerning $l$-conformal Galilei Hopf algebra \;${\cal U}^{(l)}_0(\mathcal{G})$. The
second section is devoted to its twist deformations and to the derivation of corresponding quantum space-times. The  final remarks are provided in the last section.

\section{{{$l$-conformal Galilei Hopf algebra \;${\cal U}^{(l)}_0(\mathcal{G})$}}}

In this section we recall  basic facts associated with the $l$-conformal Galilei Hopf algebra \;${\cal U}^{(l)}_0(\mathcal{G})$ provided in article \cite{czar}.
Hence, the $l$-conformal Galilei Hopf algebra includes the generators of time translations, dilatations, special conformal transformations, spatial rotations,
spatial translations, Galilei boosts
and accelerations. Denoting the mentioned generators by $H$, $D$, $K$, $M_{ij}$ and $G^{(n)}_i$, respectively, where $i=1,\dots,d$ is a spatial index and $n=0,1,\dots, 2l$,  one can write the following algebraic
\begin{eqnarray}
&&\left[\;M_{ij},M_{kl}\;\right] = -i(\delta_{ik} M_{jl}+\delta_{jl} M_{ik} - \delta_{il} M_{jk}-\delta_{jk} M_{il})\;, \label{algebra1}\\
&~~&  \cr
&&\left[\;M_{ij},G^{(n)}_k\;\right] = -i (\delta_{ik} G^{(n)}_j-\delta_{jk} G^{(n)}_i)\;, \label{algebra2}\\
&~~&  \cr
&&\left[\;K,G^{(n)}_i\;\right] = i (n-2l) G^{(n+1)}_i\;\;\;,\;\;\;
[D,G^{(n)}_i] = i (n-l) G^{(n)}_i\;, \label{algebra3}\\
&~~&  \cr
&&\left[\;H,K\;\right] = 2 i D\;\;\;,\;\;\; \left[\;D,K\;\right] = i K\;,\label{algebra4}\\
&~~&  \cr
&&\left[\;H,D\;\right] = i H\;\;\;,\;\;\; \left[\;H,G^{(n)}_i\;\right]\;=\;i n G^{(n-1)}_i\;,\label{algebra5}
\end{eqnarray}
as well as coalgebraic
\begin{eqnarray}
\Delta_0(a) = a\otimes 1 +1\otimes a\;\;\;,\;\;\;S_{0}(a) =-a\;,\label{cop}
\end{eqnarray}
sectors. It should be also observed that operators $H$, $D$ and $K$ form $so(2,1)$ subalgebra, which is the conformal algebra in one dimension. Besides, one can
notice that the instances of $n=0$ and $n=1$ in $G^{(n)}_i$ correspond to the spatial translations and Galilei boosts respectively, while the
higher values of index $n$ are linked to the accelerations. Finally, it is easy to check that all above generators are represented on the classical space of functions as follows \cite{anton}
\begin{eqnarray}
M_{ij}&\rhd& f(t,\overline{x}) =i\left( x_{i }{\partial_j} -x_{j
}{\partial_i} \right) f(t,\overline{x})\;,\label{a1}\\
H&\rhd& f(t,\overline{x})=i{\partial_t}f(t,\overline{x})\;, \label{a2}\\
G_i^{(n)}&\rhd& f(t,\overline{x})=it^n
{\partial_i}f(t,\overline{x})\;,  \label{a3}\\
D&\rhd& f(t,\overline{x})= i\left(t \partial_t+l x^i \partial_i\right)f(t,\overline{x})\;,\label{a4}\\
K&\rhd& f(t,\overline{x})=i\left(t^2 \partial_t+2 l t x^i \partial_i\right)f(t,\overline{x})\;.\label{a5}
\end{eqnarray}
Of course, for $D=K=0$ and $n=0,1$ we reproduce well-known (ordinary) Galilei Hopf structure $\;{\cal U}_{0}(\mathcal{G})$.

\section{Twist deformations of $l$-conformal Galilei Hopf algebra and the corresponding quantum space-times}

Let us now turn to the twist deformations of the Hopf structure described in pervious section. First of all, in accordance with Drinfeld  twist procedure
\cite{twist1}-\cite{twist3}, the algebraic sector of twisted
$l$-conformal Galilei Hopf algebra $\;{\cal U}^{(l)}_{\alpha}(\mathcal{G})$ remains
undeformed  (see (\ref{algebra1})-(\ref{algebra5})), while
the   coproducts and antipodes  transform as follows (see formula (\ref{cop}))
\begin{equation}
\Delta _{0}(a) \to \Delta _{\alpha }(a) = \mathcal{F}_{\alpha }\circ
\,\Delta _{0}(a)\,\circ \mathcal{F}_{\alpha }^{-1}\;\;\;,\;\;\;
S_{\alpha}(a) =u_{\alpha }\,S_{0}(a)\,u^{-1}_{\alpha }\;,\label{fs}
\end{equation}
with $u_{\alpha }=\sum f_{(1)}S_0(f_{(2)})$ (we use Sweedler's notation
$\mathcal{F}_{\alpha }=\sum f_{(1)}\otimes f_{(2)}$).
Besides, it should be noted, that the twist factor
$\mathcal{F}_{\alpha } \in {\cal U}^{(l)}_{\alpha}(\mathcal{G}) \otimes
{\cal U}^{(l)}_{\alpha}(\mathcal{G})$
satisfies  the classical cocycle condition
\begin{equation}
{\mathcal F}_{{\alpha }12} \cdot(\Delta_{0} \otimes 1) ~{\cal
F}_{\alpha } = {\mathcal F}_{{\alpha }23} \cdot(1\otimes \Delta_{0})
~{\mathcal F}_{{\alpha }}\;, \label{cocyclef}
\end{equation}
and the normalization condition
\begin{equation}
(\epsilon \otimes 1)~{\cal F}_{{\alpha }} = (1 \otimes
\epsilon)~{\cal F}_{{\alpha }} = 1\;, \label{normalizationhh}
\end{equation}
with ${\cal F}_{{\alpha }12} = {\cal F}_{{\alpha }}\otimes 1$ and
${\cal F}_{{\alpha }23} = 1 \otimes {\cal F}_{{\alpha }}$.

It is well known, that the twisted algebra $\;{\cal U}^{(l)}_{\alpha}(\mathcal{G})$ can be described in terms of
so-called classical $r$-matrix $r\in {\cal U}^{(l)}_{\alpha}(\mathcal{G}) \otimes {\cal U}^{(l)}_{\alpha}(\mathcal{G})$,
which satisfies the  classical Yang-Baxter equation (CYBE)
\begin{equation}
[[\;r_{\alpha},r_{\alpha}\;] ] = [\;r_{\alpha 12},r_{\alpha13} +
r_{\alpha 23}\;] + [\;r_{\alpha 13}, r_{\alpha 23}\;] = 0\;,
\label{cybe}
\end{equation}
where   symbol $[[\;\cdot,\cdot\;]]$ denotes the Schouten bracket
and for $r = \sum_{i}a_i\otimes b_i$
$$r_{ 12} = \sum_{i}a_i\otimes b_i\otimes 1\;\;,\;\;r_{ 13} = \sum_{i}a_i\otimes 1\otimes b_i\;\;,\;\;
r_{ 23} = \sum_{i}1\otimes a_i\otimes b_i\;.$$

In this article we consider six types of Abelian twist deformation of $l$-conformal Galilei Hopf algebra, described by the following
$r$-matrices\footnote{$a \wedge b = a \otimes b - b\otimes a.$}
\begin{eqnarray}
r_1 &=&  \frac{1}{2}{\alpha^{ij}_1} G_i^{(n)} \wedge
G_i^{(m)}\;\;\;\;\;\;\;\, [\;\alpha^{ij} = -\alpha^{ji}\;]\;,
\label{rmacierze01}\\
r_2 &=& \alpha_2 G_i^{(n)}
\wedge M_{kl}\;\;\;\;\;\;[\;i,k,l - {\rm fixed},\;\;i \neq
k,l\;]\;,\label{rmacierze02}\\
r_3 &=&  {\alpha_3} H \wedge M_{kl}\;,
\label{rmacierze03}\\
r_4 &=&  {\alpha_4} K \wedge M_{kl}\;,
\label{rmacierze04}\\
r_5 &=&  {\alpha_5} K \wedge
G_i^{(n)}\;\;\;\;\;\;\;\;[\;n=2l\;]\;,
\label{rmacierze05}\\
r_6 &=&  {\alpha_6} D \wedge M_{kl}\;,
\label{rmacierze06}
\end{eqnarray}
where $\alpha^{ij}_1$, $\alpha_2$, $\alpha_3$, ... $\alpha_6$ denote the deformation parameters.
Due to the Abelian character of the above carriers (all of them contain
 the mutually commuting elements of the algebra), the
corresponding twist factors can be obtained in a  standard way
\cite{twist1}-\cite{twist3}, i.e. they take the form
\begin{eqnarray}
{\cal F}_{{a}} = \exp
\left(ir_a\right)\;\;\;;\;\;\;a=1,2,\dots ,6\;.
\label{factors}
\end{eqnarray}

The corresponding quantum space-times are defined as the representation spaces (Hopf modules) for $l$-conformal Galilei Hopf algebra
\;${\cal U}_{\alpha}^{(l)}({{\mathcal{G}}})$, with action of the generators $M_{ij}$, $H$, $G_i^{(n)}$, $K$ and $D$ given by (\ref{a1})-(\ref{a5})
(see e.g. \cite{oeckl}, \cite{chi}). Besides, the $\star$-multiplication of arbitrary two functions covariant under
$\;{\cal U}^{(l)}_{\alpha}(\mathcal{G})$ is defined as follows
\begin{equation}
f(t,\overline{x})\star_a g(t,\overline{x}):=
\omega\circ\left(
 (\mathcal{F}_a)^{-1}\rhd  f(t,\overline{x})\otimes g(t,\overline{x})\right)
 \;,\label{star1}
\end{equation}
where symbol  $\mathcal{F}_a$ denotes the  twist factors (see
(\ref{factors}))   and $\omega\circ\left( a\otimes b\right) =
a\cdot b$. Consequently, we get
\begin{eqnarray}
&{\bf 1.~~~~~~}&[\,t,x_a\,]_{{\star_1}} =0\;\;\;,\;\;\;
[\,x_a,x_b\,]_{{\star}_1}
= i\alpha^{ij}t^{n+m}
 (\delta_{ai}\delta_{bj} -
\delta_{aj}\delta_{bi})\;,\label{rspacetime1}\\
&~~&  \cr
&{\bf 2.~~~~~~}&[\,t,x_a\,]_{{\star}_2} =0\;\label{rspacetime2}\\
&~~&  \cr
&&[\,x_a,x_b\,]_{{\star}_2} =2i\alpha_2 t^n
\left[\;\delta_{ia}(x_k\delta_{bl} - x_{l}\delta_{bk}) -
\delta_{ib}(x_k\delta_{al} -
x_{l}\delta_{ak})\;\right]\;,\nonumber\\
&~~&  \cr
&{\bf 3.~~~~~~}&[\,t,x_a\,]_{{\star}_3} =2i\alpha_3(x_k\delta_{la}-x_l\delta_{ka})\;\;\;,\;\;\;
[\,x_a,x_b\,]_{{\star}_3} =0\;,\label{rspacetime3}\\
&~~&  \cr
&{\bf 4.~~~~~~}&[\,t,x_a\,]_{{\star}_4} =2i\alpha_4t^2(x_k\delta_{la}-x_l\delta_{ka})\;,\label{rspacetime4}\\
&~~&  \cr
&&[\,x_a,x_b\,]_{{\star}_4} =4i\alpha_4 lt
\left[\;x_{a}(x_k\delta_{bl} - x_{l}\delta_{bk}) -
x_{b}(x_k\delta_{al} -
x_{l}\delta_{ak})\;\right]\;,\nonumber\\
&~~&  \cr
&{\bf 5.~~~~~~}&[\,t,x_a\,]_{{\star}_5} =2i\alpha_5t^{2(1+l)}\delta_{ia}\;,\label{rspacetime5}\\
&~~&  \cr
&&[\,x_a,x_b\,]_{{\star}_5} =4i\alpha_5 lt^{2l+1}
(x_{a}\delta_{bi} - x_{b}\delta_{ai}) \;,
\end{eqnarray}
and
\begin{eqnarray}
{\bf 6.~~~~~~~}[\,t,x_a\,]_{{\star}_6} &=&2i\alpha_6t(x_k\delta_{la}-x_l\delta_{ka})\;,\label{rspacetime6}\\
&~~&  \cr
[\,x_a,x_b\,]_{{\star}_6} &=&2i\alpha_6 l
\left[\;x_{a}(x_k\delta_{bl} - x_{l}\delta_{bk}) -
x_{b}(x_k\delta_{al} -
x_{l}\delta_{ak})\;\right]\;,\nonumber
\end{eqnarray}
respectively. It should be noted that three first spaces are the same as in the case of so-called twisted $N$-enlarged Galilei  algebra \cite{nnh}, while the remaining ones
correspond to the conformal sector of considered in present article Hopf structure. Obviously, for   deformation parameters $\alpha^{ij}_1$ and $\alpha_2,\dots , \alpha_6$ approaching
zero the above quantum  space-times  become classical.

\section{Final remarks}

In this article we consider six  Abelian twist-deformations of
$l$-conformal Galilei Hopf algebra $\;{\cal U}^{(l)}_{\alpha}(\mathcal{G})$. The corresponding twisted
space-times are derived as well. It should be noted, however, that  present studies can be extended in various
ways.  First of all, one can find the dual Hopf structures
$\,{\mathcal D}_{\alpha}^{(l)}(G)$ with the use of FRT
procedure \cite{frt} or by canonical quantization of the
corresponding Poisson-Lie structures \cite{poisson}. Besides, as it
was already mentioned in Introduction, one should ask about the
basic dynamical models corresponding to the $l$-conformal space-times
(\ref{rspacetime1})-(\ref{rspacetime6}).   Finally, one can also consider more
complicated (non-Abelian) twist deformations of
$l$-conformal Hopf algebras, i.e. one can find
the twisted coproducts, corresponding noncommutative space-times and
dual Hopf structures. Such  problems are now under consideration.

\section*{Acknowledgments}
The author would like to thank J. Lukierski
for valuable discussions. This paper has been financially  supported  by Polish
NCN grant No 2011/01/B/ST2/03354.