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