How to create 2D Gaussian noise

Share icon

Thu, 05/12/2011 - 08:14 pm

I want to create a 2D variable (Y1, Y2) whose probability density function is
1/sqrt(2*pi*sigma1*sigma2)*exp[(Y1-mu1)^2/2/sigma1^2]*exp[(Y2-mu2)^2/2/sigma2^2]

Do you have any idea how to do it? I know gnoise can create 1D Gaussian noise, but it seems Igor does not have any function creating 2D Gaussian noise.

Thanks!
Arthur
I think you need to know how to refer to rows and columns of a matrix in a wave assignment statement. Here is a tutorial. Execute the commands one-at-a-time and read the comments and notice the effect on the table and image plot.
Make/N=(3,3) mat
Edit mat
NewImage mat
mat = p         // Set each element to its row number
mat = q         // Set each element to its column number
mat = p + 10*q      // Set each element to a combination of its row and column number
mat = x + 10*y      // Initially x and y are the same as p and q
SetScale/P x -1, 1, "", mat     // Change meaning of x for wave mat. x starts at -1 and steps by 1.
SetScale/P y -1, 1, "", mat     // Change meaning of y for wave mat. y starts at -1 and steps by 1.
mat = x + 10*y              // Now x and y are different from p and q
Redimension/N=(100,100) mat
SetScale x -3, 3, "", mat       // Make x go from -3 to 3
SetScale y -3, 3, "", mat       // Make y go from -3 to 3
mat = gauss(x, 0, 1, y, 0, 1)


I believe your Y1 and Y2 are equivalent to scaled row and columns indices which in Igor wave assignment statements are represented by x and y.

For more on wave assignment statements, execute:
DisplayHelpTopic "Waveform Arithmetic and Assignment"

and
DisplayHelpTopic "Multidimensional Waves"

May 12, 2011 at 10:10 pm - Permalink

Hello Arthur,

In general, if you don't have a specific noise function you could use the Rejection method (see e.g., Numerical Recipes 3rd Ed. Chapter 7 or http://en.wikipedia.org/wiki/Rejection_sampling) which only requires that you have an expression for your distribution and access to enoise().

The expression that you provide corresponds to a 2D Gaussian with zero correlation between the two dimensions. As such you could obtain such distribution from a product of two distributions (one for each dimension). The extra complication is that gnoise() is centered at zero so you need to add an offset. For example, to get a Gaussian distribution centered around mu1 with standard deviation sigma1 you would use mu1+gnoise(sigma1).

It is not clear to me what you mean by a "2D variable (Y1, Y2)". I expect that the answer to that would be to create two 1D waves such that the first corresponds to Y1 and the second to Y2:
Make/O/N=(someNumber) Y1,Y2
Y1=mu1+gnoise(sigma1)
Y2=mu2+gnoise(sigma2)

In this case the normalized joint-histogram of the two waves will have the desired distribution.


A.G.
WaveMetrics, Inc.

May 13, 2011 at 10:16 am - Permalink

AG beat me to it. I second his remarks, and illustrate them using
function gerkin(mu1, sigma1, mu2, sigma2)
    variable mu1, sigma1, mu2, sigma2
    make /o/n=100000 data1=mu1+gnoise(sigma1), data2=mu2+gnoise(sigma2)
    make /o/n=(50,50) myHist
    setscale x,mu1-3*sigma1, mu1+3*sigma1,myHist//set histogram limits to mu +/- 3 sigma
    setscale y,mu2-3*sigma2,mu2+3*sigma2, myHist
    JointHistogram(data1,data2,myHist)
end


•gerkin(1, 1.5, 2, 5)


See see http://www.igorexchange.com/node/1373 for the JointHistogram function. If the 2D joint Gaussian variables are correlated, the procedure is more complicated, but doable.
Hist2D.png

May 13, 2011 at 11:00 am - Permalink

Here is an example for two correlated Gaussian variables I had worked out a while ago. For convenience, they have zero mean, but the extension is obvious. It also uses the Wavemetrics #include procedures for BiVariate histograms. The calculation relies on using the orthogonalized Gaussian random variables, which are uncorrelated. See any standard text on statistics, stochastic processes, or signal processing (e.g. Van Trees).
#pragma rtGlobals=1     // Use modern global access method.
#include<Bivariate Histogram>
#include<Bivariate Histogram 2>

function BIVAR(vN, rho)
    variable vN, rho // number of samples, correlation coefficient
    variable/G vNG = vN
    variable/G rhoG = rho
    WAVE BivariateHistWave
    variable sig1sq =   1 + rho//   variance for independent x1
    variable sig2sq =   1  - rho//  variance for independent x2
   
    make/O/N=(vN)   x1, x2, Y1, Y2
    x1  =   gnoise(sqrt(sig1sq))//  Gaussian random generator
    x2  =   gnoise(sqrt(sig2sq))
   
    Y1  =   ( x1 + x2 ) / sqrt(2)// joint (CORRELATED) data
    Y2  =   ( x1  - x2 ) / sqrt(2)//    obeying given pdf
    BiHistSetBins(Y1, Y2, 201, -5, 0.05, 201, -5,   0.05, "BivariateHistWave")
    BivariateHistWave = log(BivariateHistWave)
    BivariateHistWave = numtype(BivariateHistWave[p][q])==1 ? -0.1:bivariatehistwave[p][q]//change -inf's to -0.1  
end

Note that for display purposes, I have plotted the log of the histogram values. The graph recreation procedure for displaying the histogram is (be careful with possible line wraps)
Window Graph0() : Graph
    PauseUpdate; Silent 1       // building window...
    Display /W=(264,48.5,594,353.75)
    AppendImage BivariateHistWave
    ModifyImage BivariateHistWave ctab= {*,*,Grays,1}
    ModifyGraph height={Aspect,1}
    ModifyGraph grid=2
    ModifyGraph mirror=2
    ModifyGraph lblMargin(left)=6
    ModifyGraph standoff=0
    ModifyGraph lblRot(left)=-90
    Label left "R\\B2"
    Label bottom "R\\B1"
    SetAxis left -5,5
    SetAxis bottom -5,5
    TextBox/C/N=text0/F=0/H=50/A=MT/X=-4.59/Y=7.02/E "Bivariate Histogram with \\{vNG} samples, \\F'Symbol' r\\]0 =\\{rhoG}"
    AppendText "\\JC( log\\B10\\M scaling, white = 0 unscaled bin content)"
    ColorScale/C/N=text1/F=0/H=50/A=RC/X=3.41/Y=-4.22/E image=BivariateHistWave
EndMacro

A typical graph is attached, from executing
BIVAR(100000, 0.90)
BiVariateHistogram.png

May 14, 2011 at 06:17 am - Permalink

The efficient way to make correlated Gaussian random variables, with arbitrary covariance matrix, uses Cholesky decomposition. Here's a fragment we have used. It would need to be cleaned up by some helpful person to be a proper snippet, but it's a starting point:

Function MakeCorrRV(sigma,corr_range,dat_range,m,n) // generates a set of corr. noise trials
    variable sigma,corr_range,dat_range,m,n             // m = # obs. points, n = # trials
    make /d/o/n=(m,m) corrmat                       // correlation matrix
    variable corr_range1 = corr_range*m/(2*dat_range)   // convert from x-units to points
    corrmat = exp(-abs(p-q)/corr_range1)                // exponential kernel
    MatrixOp /o Cholesky=chol(corrmat)                  // Cholesky decomposition
    Make /free/n=(n,m) uncorr=gnoise(1)                 // uncorrelated data
    MatrixOp /o corrnoise = uncorr x Cholesky               // generate the n correlated RVs   
    MatrixOp /o corrtest = syncCorrelation(corrnoise)       // check that RVs have right cov
    corrnoise *= sigma                          // scale noise by sigma
    SetScale/I x -dat_range,dat_range,"", corrnoise, corrmat,corrtest
End


See Numerical Recipes or wikipedia: http://en.wikipedia.org/wiki/Cholesky_decomposition [under Monte Carlo]


John Bechhoefer
Department of Physics
Simon Fraser University
Burnaby, BC, Canada

May 20, 2011 at 07:55 pm - Permalink