Monte Carlo Uncertainty Propagation

Share icon

I try to perform a Monte Carlo Uncertainty Propagation.
I use correlated random variables using the known covariance Covar of a set of 4 parameters Pi and the cholesky decomposition of the (4,4) Covar.
The calculation being a bit complicated, I perform a loop over the NIterations, creating a set of output variables Oi[Niterations] and then use wavestats to extract the average and sdev of the Oi.

Is it faster in general to generate a set of 4 x NIterations random variables before the loop, or to draw the 4 parameters at each iteration of the loop? 

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

November 21, 2013 at 08:40 pm - Permalink