I've added a picture of what I want to compute. In the nomenclature of the picture, I want to compute a realization of y(x) given the known distributions and constants.
Let's say y(x), random variable, represents the height of a line, with x being the usual x-axis.
The probability distribution function is:
p(y)=K * exp[−(y^2/(2 * s^2))) i.e. a Gaussian and K and s are known constants
And the autocorrelation of y(x) relative to y(x') is given as:
C = s^2 * exp[−x^2/(2 * R^2)] --> C is the autocorrelation function between y(x) and y(x'), which since the process is stationary will only depend upon the difference between x and x' (i.e. x-x'). Here I replace (x-x') by x.
where R is also a known constant
What I want to do is to compute a realization of y(x). How should I go about doing this?
Given the constraints of the known distributions and constants (s, R, K) I would like to be able to compute a realization of y.
Suppose you want to compute a realization at times $t_1<t_2<\dots<t_N$. This amounts to sampling from the distribution of an $N$ dimensional Gaussian vector $X$ with mean $0$ and covariance matrix $R$ where $R_{ij}=C(t_i-t_j)$. Such a vector is given by $X=R^{1/2} Y$, where $Y$ is an $N$ dimensional vector of independent standard normal variables and $R^{1/2}$ denotes the matrix square root of $R$. $R^{1/2}$ can be computed using an eigenvalue decomposition of $R$: if $R=P D P^{-1}$ for a diagonal $D$ then $R^{1/2}=P D^{1/2} P^{-1}$, where the matrix square root of a diagonal matrix is obtained by taking square roots of the diagonal entries. This decomposition always exists because $R$ is symmetric, and the matrix square root is always real because $R$ is always nonnegative definite.
If you want to sample $Y$ "by hand", a standard algorithm for the job is called the Box-Muller algorithm. But one usually just uses a built-in library function, usually called randn, for this task.