How to generate uniform distributed samples with given auto-correlation function As I mentioned in the question title, I want to generate specific uniformly distributed samples.
I need them to model a real world scenario. For my real data, I estimated a function, which approximates the auto-correlation function by a e-function. I also modeled a distribution model for my real data.
I can generate samples, which follow my distribution model with the inverse transform sampling.
If i think correct, instead of "standard" uniform distributed samples, which i feed into the inverse transform sampling, i need to adapt those uniform distributed samples, so that their auto-correlation function follows my estimated e-function.
Now the question is, if i am correct, how to adapt my uniform distributed samples, that they correspond to my requirements.
It would be great if, somebody could help.
EDIT 1:
Here is a plot of the auto-correlation function of my real data (different scenarios) and the modeled function (green). The ACF of my uniform distributed samples should follow the green (respectively blue) line.

EDIT 2:
Thank you so much so far for the answers. To make my problem a bit clearer, i added a illustration of the problem. For a simulation, i need to generate samples, which follow my distribution model (with inverse transform sampling a easy task) AND have the same auto-correlation behavior like the real data.
My idea was, that i need to adapt the uniformly distributed samples, which i feed into the inverse transform sampling, so that their ACF follows my model.
The problem is, that i don't know if this idea is correct and that i don't know how to adapt the uniformly distributed samples.

 A: Letting $U_t=\Phi(X_t)$ where $X_t$ is a zero-mean and unit-variance stationary Gaussian process with autocorrelation function $\rho_h$ and $\Phi$ is the standard normal cdf, it follows that each $U_t$ is marginally uniformly distibuted.  It follows that the relation between $\rho_h$ and the autocovariance function of $U_t$ is
\begin{align}
\gamma_h&=\operatorname{Cov}(U_t,U_{t+h})
\\&=E(U_tU_{t+h})-E(U_t)E(U_{t+h})
\\&=E(\Phi(X_t)\Phi(X_{t+h}))-1/4
\\&=P(Z_1\le X_t\cap Z_2\le X_{t+h})-1/4
\\&=P\left(\frac{Z_1-X_1}{\sqrt{2}}\le0\cap \frac{Z_2-X_{t+h}}{\sqrt{2}}\le 0\right)-1/4
\\&=\Phi_2\left(\begin{bmatrix}0\\0\end{bmatrix};\frac{\rho_h}2\right)-1/4.
\end{align}
Here $Z_1$ and $Z_2$ are independent standard normal random variables, and $\Phi_2$ is the cdf of the standard bivariate normal distribution with correlation $\rho_h/2$.
Using this result, the above relation simplifies to
$$
\gamma_h = \frac1{2\pi}\operatorname{arcsin}\frac{\rho_h}2.
$$
Solving this equation, we find that the autocorrelation of $X_t$ needed to achieve a target autocovariance $\gamma_h$ of $U_t$ is
$$
\rho_h = 2\sin(2\pi \gamma_h).
$$
Even if the target autocovariance function $\gamma_h$ is positive semi-definite, the above construction may not be feasible.  For example, if the target autocovariance function is that of a MA(1)-process with a unit root,
$$
\gamma_h=\begin{cases}
1/12 &\text{for }h=0 \\
1/24 &\text{for }h=1 \\
0 &\text{for }h>1
\end{cases},
$$
this would imply that
$$
\rho_h=\begin{cases}
1 &\text{for }h=0 \\
0.518 &\text{for }h=1 \\
0 &\text{for }h>1
\end{cases}
$$
which is not a positive semi-definite autocorrelation function.
A: Here's a pragmatic and easy approach with room to expand and establish proofs, with the focus on the main problem: how do you generate a correlated uniform sample?
Let $U_1, U_2$ be uniformly distributed and independent on the unit interval. Let $V_1 = U_1$. For a desired correlation $d$, let $G_{1,2}$ be yet another uniformly distributed random variable.
Let $$V_2 = \left\{ \begin{array}{ccc} U_1 & \text{if} & G_{1,2} < d \\
U_2 & \text{if} &G_{1,2} \ge d\end{array} \right.$$
I claim that:

*

*One can show that $V_2$ is indeed uniformly distributed, with a correlation of $d$ with $V_1$

*By way of induction, one can set an arbitrary sequence of uniform random variables and create a new sample having a desired covariance structure.

As usual, this site is always most compelled by code, so I can show a couple cases of the spherical AR-1 auto correlation, where the correlation between observations with a lag of 1 is set to $d$, but otherwise it's relatively straightforward to use any structure you want.
do.one <- function(n,N,d) {
  u <- matrix(runif(n*N), N, n)
  g <- matrix(runif((n-1)*N/2), N, n-1)

  v <- u
  for( i in 2:n) {
    v[g[, i-1] < d,i] <- v[g[, i-1] < d, i-1]
  }  
    
  acf <- sapply(2:n, function(i) cor(v[,i], v[, 1]))
  acf
}

set.seed(123)
ds <- c(0,0.25, 0.5, 0.95)
acfs <- sapply(ds, do.one, n=10, N=1000)
matplot(acfs,type='l', xlab='Lags')
legend('topright', title = 'AR-1 spherical correlation', legend = ds, lty=1:4, col=1:4)


A: If I understand correctly your problem, you need to simulate random variables which marginally follow an uniform distribution while the joint distribution is a multivariate uniform distribution with some correlation between your marginals.
Formally:
$X_i\sim U(0,1), \forall_i, i = 1,2,...,K$ , where $X\sim U(a_i, b_i, Cor)$.
In the method I am firstly proposing, your uniform random variables will be uniformly distributed between 0 and 1. However, I also provide a "Normal to Anything" kind of approach.
One trick is to use generate a Multivariate normal distribution, specifying a correlation matrix. Let's say $\Sigma$ is your auto-correlation structure for your $K$ random variables.
$$X \sim MVN(\mu, \Sigma)$$
$\mu$ is a vector of means of size $K$, while $\Sigma$ is a $K*K$ matrix
Then you have to transform your quantiles to probabilities from any normal distribution you want, for one $X$:
$$CDF_{Normal}(X_{i}, \mu_i, \Sigma_{ii}) \sim U(0,1)$$
Now if you assume that your variables are not uniformly distributed between a common support $[a,b], \forall K$ but a variable-dependent support instead. You just need to convert your probabilities obtained previously with an inverse uniform distribution specifying a and b.  For a random variable,
$$\theta^{-1}(CDF_{Normal}(X_{i}, \mu_i, \Sigma_{ii}), a_i, b_i)$$
It is worth noting that you will respect globally your correlation structure, while under-estimate it a little bit.
I can provide script to illustrate the method.
Hope this helps.
A: You can try implying AR(p) process coefficients $\phi_i$ from the given ACF $r(p)$. You could apply Yule Walker equations:

*

*form a vector $r$ of ACF for lags $p$: $1, r_1, r_2,\dots, r_p$

*construct a correlation matrix $R$ as described in the link above from $r_p$, e.g. the third row would be $(r_2,r_1,1,p_1,\dots,r_{p-2})$

*calculate $\phi=R^{-1}r$
Use these coefficients to produce autocorrelated samples
A: The following spells out the details of the approach proposed in the other answer by @AdamO and in its comments by @LucaCiti.
For $i=1,2,\dots,\infty$, let $|\phi_i|$ denote the probability that $U_t$ takes a value identical to either $U_{t-i}$ or $1 - U_{t-i}$ and let these two possibilities be determined by the sign of $\phi_i$.  Let the remaining fraction
$$
\phi_0=1-\sum_{i=1}^\infty |\phi_i|
$$
denote the probability that $U_t$ takes a uniformly distributed value independent of the history of the process.  Clearly, we must have
$$
0\le \phi_0\le 1. \tag{1}
$$
and
$$
-1\le \phi_i\le 1 \tag{2}
$$
for $i=1,2,\dots,\infty$.
Letting $V_t=U_t-\frac12$ denote the mean-centered process, and using the law of total expectation, we have
\begin{align}
E(V_t|V_{t-1},V_{t-2},\dots)
&=|\phi_1|\operatorname{sgn}(\phi_1) V_{t-1} + |\phi_2|\operatorname{sgn}(\phi_2) V_{t-2} + \dots
\\&=\phi_1 V_{t-1} + \phi_2 V_{t-2} + \dots.
\end{align}
Thus it is immedeately clear that $\phi_1,\phi_2,\dots$ are the coefficients in the  $\operatorname{AR}(\infty)$ representation of the model.
Unlike ordinary ARMA models, constraint (1) and (2) implies that not all positive semi-definite autocovariance functions are possible via this construction, however. For example, if the target autocovariance function is that of a MA(1) model with MA polynomial $1-\theta B$, the infinite AR polynomial would equal
$$
\frac1{1-\theta B}=1+\theta B+\theta^2 B^2+\dots,
$$
and we would have
$$
\sum_{i=1}^\infty|\phi_i|=\sum_{i=1}^\infty |\theta^i|=\sum_{i=1}^\infty |\theta|^i = \frac{|\theta|}{1-|\theta|}.
$$
Combined with (1) this limits possible values of $\theta$ to
$$
|\theta|\le \frac12
$$
and the correlation at lag 1 to
$$
-\frac25\le \rho_1=\frac{\theta}{1+\theta^2}\le\frac25
$$
In contrast, via the copula described in my other answer the correlation at lag 1 is limited to $|\rho_1|<0.4825837$ only .  Semipositive definiteness in itself limits the same correlation to $|\rho_1|\le 1/2$.
A: Uniform distributed samples are the set of samples in which every element is distributed uniformly, if we place further constraints on this, it ceases to be uniform distributed samples. But anyway, if you meanе something else, we can easily upgrade my answer.
Let us just to generate samples with given autocorrelation function. Our idea is put all needed constraints on samples and let scipy.optimize to do everything for us.
While the key concept is simple, optimization and numerical problems can be arbitrary hard in real application and it can be necessary to adjust scipy solvers parameters here and/or optimize some computations.
The code implements this idea:
import numpy as np
import scipy.optimize as sco
import matplotlib.pyplot as plt

def cov(x,y):
    return np.sum(x*y)/len(x)

def cor(x,y):
    return cov(x,y)/(np.std(x)*np.std(y))

#define our autocorrelation
def autoMoment(x, moment, n):
    out = []
    for i in range(n):
        out = [moment(x,x)]
        for i in range(1,n):
            out.append(moment(x[i:],x[:-i]))
    return np.array(out)

#define generation of function for scipy.optimize
def generateWithSpecificAutomoment(moment, values): 
    def f(x):
        out = autoMoment(x, moment, len(values))
        return np.sum((np.array(out) - np.array(values))**2)
    return f

desiredAutocorr = [1,0.3,0.1,-0.23,0.03,0.07,-0.03]

initial = np.random.randn(50) #starting point for the solver
solution = sco.minimize(generateWithSpecificAutomoment(cor, desiredAutocorr), initial)

print(solution.success)
out = solution.x #samples obtained after optimization

#check autocorrelation on generated samples
realAutocorr = autoMoment(out, cor, len(desiredAutocorr))

#compare desired and real result
plt.title("Autocorrelations")
plt.plot(desiredAutocorr, linewidth=7, label="desired")
plt.plot(realAutocorr, label="real")
plt.legend()
plt.show()

Generated picture:

A: This problem can be approached by first generating samples from the desired distribution and then reordering them to match the desired autocorrelation function.
The R code below demonstrates an approach based on this answer that can be modified for any desired ACF and distribution. The example generates $n=10^6$ samples from $\text{Gamma}(0.9,1)$ whose ACF follows 50 random samples from $\text{Beta}(1,3)$, sorted descending.
The process is as follows.

*

*Generate the desired number of samples, $X=\{x_1,...,x_n\}$, from the target distribution. Set $\alpha_0$ equal to the desired ACF. Initialize the target ACF, $\alpha$, as the desired ACF.

*Find a set of weights that, when passed to filter along with $n$ random normal variates, results in a series, $Y$, with $ACF=\alpha$ (see the answer linked above).

*Reorder $X$ so that its rank ordering matches that of $Y$. If $X$ is normally distributed, the resulting series should have the desired ACF; however, the more $X$ deviates from normality, the more the ACF will deviate from $\alpha_0$ (the example below has a target distribution of $\text{Gamma}(0.9,1)$, which is very "non-normal"). Update the target ACF, $\alpha$, according to $\alpha'=\frac{\alpha}{2}\Big(\frac{\alpha_0}{ACF}+1\Big)$ and repeat steps 1-3 until the ACF of the reordered $X$ converges.

The function that performs the reordering (it works only for positive values for alpha):
acf.reorder <- function(x, alpha) {
  tol <- 1e-5
  maxIter <- 10L
  n <- length(x)
  xx <- sort(x)
  y <- rnorm(n)
  w0 <- w <- alpha1 <- alpha
  m <- length(alpha)
  m1 <- 1:(m - 1)
  tol10 <- tol/10
  iter <- 0L
  x <- xx[rank(filter(y, w, circular = TRUE))]
  SSE0 <- Inf
  
  while ((SSE <- sum((acf(x, lag.max = m - 1, plot = FALSE)$acf[1:m] - alpha)^2)) > tol) {
    if (SSE < SSE0) {
      SSE0 <- SSE
      w <- w0
      if ((iter <- iter + 1L) == maxIter) break
    } else break
    
    w1 <- w0
    a <- 0
    sse0 <- Inf
    
    while (max(abs(alpha1 - a)) > tol10) {
      a <- c(1, sapply(m1, function(i) sum(head(w1, -i)*tail(w1, -i)))/sum(w1^2))

      if ((sse <- sum((a - alpha1)^2)) < sse0) {
        sse0 <- sse
        w0 <- w1
      } else {
        # w0 failed to converge; try optim
        w1 <- exp(
          optim(
            log(w0),
            function(ww) {
              ww <- exp(ww)
              sum((c(1, sapply(m1, function(i) sum(head(ww, -i)*tail(ww, -i)))/sum(ww^2)) - alpha1)^2)
            },
            method = "L-BFGS-B"
          )$par
        )
        
        a <- c(1, sapply(m1, function(i) sum(head(w1, -i)*tail(w1, -i)))/sum(w1^2))
        sum((a - alpha1)^2) < sse0
        if (sum((a - alpha1)^2) < sse0) w0 <- w1
        break
      }
      
      w1 <- (w1*alpha1/a + w1)/2
    }
    
    x <- xx[rank(filter(y, w0, circular = TRUE))]
    alpha1 <- (alpha1*alpha/acf(x, lag.max = m - 1, plot = FALSE)$acf[1:m] + alpha1)/2
  }
  
  xx[rank(filter(y, w, circular = TRUE))]
}

Generate samples from the target distribution and specify the desired ACF:
set.seed(1960841256)
x <- rgamma(1e6, 0.9, 1)
alpha <- c(1, sort(rbeta(50, 1, 3), TRUE))

Reorder x and plot its ACF against alpha:
x <- acf.reorder(x, alpha)
acf(x, lag.max = length(alpha) - 1)
lines(seq_along(alpha) - 1, alpha, col = "green")


The resulting ACF is a good match with the target, and since x has simply been reordered, it is known to have the desired distribution.
