Correlated Gaussian Random Variables

Function gnoise1(var,n) // n draws of m-dim vector of corr. noise wave var; variable n // var must be sym, pos-def. matrix variable m=dimsize(var,0) MatrixOp /free cholesky=chol(var) // compute Cholesky decomposition make /d/free/n=(n,m) uncorr=gnoise(1) // uncorrelated N(0,1) data MatrixOp /o corrdat = uncorr x cholesky End

Function Demo2d(s1,s2,rho,nsamples,) variable s1,s2,rho,nsamples // correlation coeff: -1 < rho < 1 make /d/o/n=(2,2) var={{s1^2, rho*s1*s2},{rho*s1*s2,s2^2}} // make 2x2 covariance matrix gnoise1(var,nsamples); wave corrdat Display corrdat[*][1] vs corrdat[*][0] ModifyGraph height={Plan,1,left,bottom}, mode=2, lSize=2 End

andyfaff

Bech's code snippet can be used to visualise covariance matrices obtained from a Curvefitting algorithm, typically

Funcfit

.
It needs to be adapted a bit in this case, as the covariance matrix has NaN entries for the parameters that are being held. The following code snippet takes the fitted coefficients from Funcfit, along with the covariance matrix, and the holdstring (specifying which parameters are being held), and creates 2D wave called M_montecarlo.
This wave has howmany rows, and numpnts(coefs) columns.
Each of the rows contains a set of parameters which could be fed to the fit function used. Futhermore, for a given (parameter) column the mean and standard deviation of that column is equal to the fitted parameter and it's uncertainty. Moreover, different columns possess the correct covariance as specified by the input covariance matrix.


Function gen_synthMCfmCovar(coefs, M_covar, holdstring, howmany)
        //coefs from the curvefitting
        //M_covar from the curvefitting output
        //holdstring - the holdstring used for the fit.
        //how many samples you want to synthesize.

	wave coefs, m_covar
	string holdstring
	variable howmany

	variable varying = 0, ii, whichcol
	duplicate/free coefs, tempcoefs
	duplicate/free M_covar, tempM_covar

	varying = strlen(holdstring)
	for(ii = dimsize(coefs, 0)-1 ; ii >= 0 ; ii-=1)
		if(str2num(holdstring[ii]))
			deletepoints/M=0 ii, 1, tempcoefs, tempM_covar
			deletepoints/M=1 ii, 1, tempM_covar
			varying -=1
		endif
	endfor

	//create some gaussian noise
	make /free/d/n=(varying, howmany) noises = gnoise(1., 2)
	//and create the correlated noise from the covariance matrix.
	matrixop/free correlatedNoises = (chol(tempm_covar) x noises)^t

	make/n=(howmany, dimsize(coefs, 0))/d/o M_montecarlo
	Multithread M_montecarlo[][] = coefs[q]

	//now add the correlated noise back to the parameters
	whichcol = dimsize(correlatedNoises, 1) - 1
	for(ii = dimsize(coefs, 0)-1 ; ii >= 0 ; ii-=1)
		if(!str2num(holdstring[ii]))
			M_montecarlo[][ii] += correlatedNoises[p][whichCol]
			whichCol -= 1
		endif
	endfor

End

July 14, 2013 at 04:54 pm - Permalink

s.r.chinn

Very nice snippet!

I couldn't resist the temptation to try it in 3D. The new command line, demo function and Gizmo macro are:


•Demo3d(1,1,1,0.9,1e4)

Function Demo3d(s1,s2,s3,rho,nsamples,)	
	variable s1,s2,s3,rho,nsamples				       // correlation coeff:  -1 < rho < 1
	make /d/o/n=(3,3) var2={{s1^2, rho*s1*s2, rho^2*s1*s3},{rho*s1*s2,s2^2, rho*s2*s3}, {rho^2*s1*s3,rho*s2*s3,s3^2}}	   // make 3x3 covariance matrix
	gnoise1(var2,nsamples);
	wave corrdat
	Execute "Gizmo0()"
End

Window Gizmo0() : GizmoPlot
	PauseUpdate; Silent 1	// Building Gizmo 6 window...

	// Do nothing if the Gizmo XOP is not available.
	if(exists("NewGizmo")!=4)
		DoAlert 0, "Gizmo XOP must be installed"
		return
	endif

	NewGizmo/N=Gizmo0/T="Gizmo0"
	ModifyGizmo startRecMacro
	MoveWindow 17.25,53.75,317.25,289.25 
	AppendToGizmo Scatter=root:corrdat,name=scatter0
	ModifyGizmo ModifyObject=scatter0 property={ scatterColorType,3}
	ModifyGizmo ModifyObject=scatter0 property={ markerType,0}
	ModifyGizmo ModifyObject=scatter0 property={ sizeType,0}
	ModifyGizmo ModifyObject=scatter0 property={ rotationType,0}
	ModifyGizmo ModifyObject=scatter0 property={ Shape,1}
	ModifyGizmo ModifyObject=scatter0 property={ size,1}
	ModifyGizmo ModifyObject=scatter0 property={ markerCTab,Rainbow}
	ModifyGizmo ModifyObject=scatter0 property={ CTABScaling,16}
	AppendToGizmo Axes=boxAxes,name=axes0
	ModifyGizmo ModifyObject=axes0,property={0,axisRange,-1,-1,-1,1,-1,-1}
	ModifyGizmo ModifyObject=axes0,property={1,axisRange,-1,-1,-1,-1,1,-1}
	ModifyGizmo ModifyObject=axes0,property={2,axisRange,-1,-1,-1,-1,-1,1}
	ModifyGizmo ModifyObject=axes0,property={3,axisRange,-1,1,-1,-1,1,1}
	ModifyGizmo ModifyObject=axes0,property={4,axisRange,1,1,-1,1,1,1}
	ModifyGizmo ModifyObject=axes0,property={5,axisRange,1,-1,-1,1,-1,1}
	ModifyGizmo ModifyObject=axes0,property={6,axisRange,-1,-1,1,-1,1,1}
	ModifyGizmo ModifyObject=axes0,property={7,axisRange,1,-1,1,1,1,1}
	ModifyGizmo ModifyObject=axes0,property={8,axisRange,1,-1,-1,1,1,-1}
	ModifyGizmo ModifyObject=axes0,property={9,axisRange,-1,1,-1,1,1,-1}
	ModifyGizmo ModifyObject=axes0,property={10,axisRange,-1,1,1,1,1,1}
	ModifyGizmo ModifyObject=axes0,property={11,axisRange,-1,-1,1,1,-1,1}
	ModifyGizmo ModifyObject=axes0,property={-1,axisScalingMode,1}
	ModifyGizmo ModifyObject=axes0,property={-1,axisColor,0,0,0,1}
	ModifyGizmo modifyObject=axes0 property={Clipped,0}
	AppendToGizmo freeAxesCue={0,0,0,1},name=freeAxesCue0
	ModifyGizmo setDisplayList=0, opName=ortho0, operation=ortho, data={-2.5,2.5,-2,2,-2,2}
	ModifyGizmo setDisplayList=1, object=freeAxesCue0
	ModifyGizmo setDisplayList=2, object=scatter0
	ModifyGizmo setDisplayList=3, object=axes0
	ModifyGizmo SETQUATERNION={-0.197844,0.141310,0.028073,0.969590}
	ModifyGizmo autoscaling=1
	ModifyGizmo currentGroupObject=""
	ModifyGizmo compile

	ModifyGizmo userString={wmgizmo_df,"Gizmo0"}
	ModifyGizmo endRecMacro
End

The Gizmo plot should resemble the attached figure. (The Gizmo Macro and figure have been updated to include a display enhancement suggested by WaveMetrics' AG.)

Attachments Attachment Gizmo0ExportedNew.pdf (15.28 KB)

July 17, 2013 at 03:24 am - Permalink

Correlated Gaussian Random Variables

Igor Pro 10

Igor XOP Toolkit

Igor NIDAQ Tools MX