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.
EDIT 3: Following the approach from jblood94, the solution is based on a reordering of the output samples after the inverse transform sampling. So, the underlying distribution stays the same AND the variables follow the desired acf.
I tried to translate the R code from jblood94 to matlab:
clear;
close all;
clc;
x = gamrnd(0.9, 1, 1e6, 1);
alpha = [1, sort(betarnd(1, 3, 1, 50), 'descend')];
tic
xx = acf_reorder(x, alpha);
toc
figure;
autocorr(xx, 'NumLags', length(alpha)-1);
hold on;
plot((1:length(alpha))-1, alpha);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xx = acf_reorder(x, alpha)
tol = 1e-5;
maxIter = 10;
n = length(x);
xx = sort(x);
y = randn(n,1);
[w0, w, alpha1] = deal(alpha);
m = length(alpha);
m1 = 1:(m-1);
tol10 = tol/10;
iter = 0;
x = xx(tiedrank(filter(w, 1, y)));
SSE0 = inf;
acfInit = autocorr(x, 'NumLags', m-1)';
SSE = sum((acfInit(1:m)-alpha).^2);
while SSE > tol
if(SSE < SSE0)
SSE0 = SSE;
w = w0;
iter = iter + 1;
if(iter == maxIter)
break;
end
else
break;
end
w1 = w0;
a = 0;
sse0 = inf;
while (max(abs(alpha1 - a)) > tol10)
a = [1, arrayfun(@(x) sum(head(w1, -x).*tail(w1, -x)), m1)/sum(w1.^2)];
sse = sum((a - alpha1).^2);
if (sse < sse0)
sse0 = sse;
w0 = w1;
else
% w0 failed to converge; try optim
disp('w0 failed to converge; try optim');
w1 = exp(fminsearch(@fun, log(w0), [], m1));
a = [1, arrayfun(@(x) sum(head(w1, -x).*tail(w1, -x)), m1)/sum(w1.^2)];
if (sum((a - alpha1).^2) < sse0)
w0 = w1;
break;
end
end
w1 = (w1.*alpha1./a + w1)/2;
end
x = xx(tiedrank(filter(w0, 1, y)));
acf = autocorr(x, 'NumLags', m-1);
alpha1 = (alpha1.*alpha./acf(1:m) + alpha1)/2;
SSE = sum((acf(1:m)-alpha).^2);
end
xx = xx(tiedrank(filter(w, 1, y)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function output = fun(ww, m1)
ww = exp(ww);
output = sum([1, arrayfun(@(x) sum(head(ww, -x).*tail(ww, -x)), m1)/sum(ww.^2)]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dataOut = head(dataIn, num)
dataIn = dataIn(:);
if num>0
dataOut = dataIn(1:num, :);
elseif num<0
dataOut = dataIn(1:(length(dataIn)+num));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dataOut = tail(dataIn, num)
dataIn = dataIn(:);
if num>0
dataOut = dataIn((length(dataIn) - num + 1):length(dataIn),:);
elseif num<0
dataOut = dataIn((-num+1):length(dataIn),:);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This yields SOMETIMES plots like this, and sometimes it's trapped in the inner while loop (I need to fix this). As you can see, the plot shows, that the result is not ideal.
EDIT 4:
Maybe I didn't formulate my problem correctly at first: I want to generate a time series that follows a defined distribution (in my case approximately Half-Gaussian) and whose autocorrelation function follows a given autocorrelation function (e.g. up to lag 30). So I have two properties that my time series must fulfil.
I tried the following approaches (without success):
1)
- generation of uniformly distributed samples
- invese transform sampling with desired distribution
- Reordering the samples to match the ACF
- Failed: it takes far too long for my time series (1e6 samples) + acf could not be reached.
- generation of normally distributed samples with desired autocorrelation (with the approach from jblood94)
- invese transform sampling with uniform distribution (based on this solution)
- invese transform sampling with desired distribution
- Failed: the autocorrelation of the transformed samples does not correspond to the desired Acf, which was correct before the transformation.
I think approach 2 might work. However, it must be prevented that the autocorrelation is changed by the transformation. Or the initial autocorrelation must be adjusted in such a way that it matches the desired Acf again, at least after the transformation.
Now the question is how the initial autocorrelation needs to be adjusted. Or is the approach correct? First generate normally distributed samples with autocorrelation, transform them into normally distributed samples, and then transform them into my desired distribution?
In addition, so far I have only managed to change individual lags of the normal distribution, but not that lag 1-30, for example, follow my function.
EDIT 5:
I tested approach 2 again. It seems, that the change in the autocorrelation is dependent from the acf itself. I generated samples, which follow now my distribution and my acf. I will post plots later.
This approach seems to work for me, but I can't guarantee, that it works in general.