Monte Carlo Uncertainty Propagation
millot1
Thu, 11/21/2013 - 05:11 pm
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?
//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