Code for Standard Deviation, Averages and r statistic

Share icon

Sun, 07/21/2013 - 09:20 pm

Hello,
I have completed numerous experiments for analysis of H2O (water vapour mixing ratio) and D_DH (deuterium). For these experiments (often lasting in excess of 8 hours with 5 sec data points), I have performed variable temperature increases which are sometimes 40 mins, sometimes 60 mins, sometimes longer. These temperature increases cause differences in the data which I’m testing for to understand the effect of temperature on H2O and D_DH.

Currently I’m using WaveStats to generate the average and standard deviation over each temperature increase as well as the r statistic (StatsLinearCorrelationTest/T=1/Q H2O,D_DH). However, currently I’m doing this process manually, which is not only very time consuming, but can be prone to errors during data entry.

I’m wondering if anyone has an Idea to formulate a code with the input of a data point range which captures each temperature range corresponding to a time I have recorded for each temp setting (i.e. data points 122 – 350) that returns the standard deviation, average vale and ‘r’ statistic over the designated range of points. As well as this I would need to repeat this process multiple times for each temperature range in an experiment.

Any help would be greatly appreciated!
mean(dataWave, 122, 350) and sqrt(variance(dataWave, 122, 350)) will give you the mean and s.d. of a wave for 122 < x < 350. Here x is the scaled x value. If you use default wave scaling then x is identical to the point numbers, i.e. deltax(dataWave) == 1 and leftx(dataWave) == 0.

July 21, 2013 at 11:53 pm - Permalink

Yeah, I'd say this type of problem crops up so frequently, the solution probably gets reinvented every time somebody new learns how to program an igor function.

Here's my version that I end up using all too often. You'll need to modify this to meet your own needs (add a second raw data wave, add in the code to do StatsLinearCorrelationTest on subsets of the raw data, do anything else...?). I don't know if you started down the road of learning to write a function yet, but if not this will get you started:

function processstats(raw,startpts,endpts)
    wave raw,startpts,endpts
    variable idx
    make /o /d /n = (numpnts(startpts)) avgs, sdevs
    for (idx = 0; idx < numpnts(startpts); idx += 1)
        wavestats /Q /R=[startpts[idx],endpts[idx]] raw
        avgs[idx] = V_avg
        sdevs[idx] = V_sdev
    endfor
    return 0
End

July 22, 2013 at 08:32 am - Permalink

Here's a function that will do the job of iterating.

You need to set it up by creating two waves expStart and expStop which represent the start (non-inclusive) and stop (inclusive) times of each desired averaging period. e.g., if you did an experiment from 11:00 to 12:00 today:

expStart = {date2secs(2013,07,24) + 11*3600}
expStop = {date2secs(2013,07,24) + 12*3600}


will need to be done before running averageExperiments(data, dataTimes, expStart, expStop).

I'll leave the addition of the r-values code to you...

Note that using Extract means your times don't have to be ordered, but it makes the function way slower on large waves. The second helper function isn't really necessary but I had used it in my own function so I just threw it in here.

#pragma rtGlobals=3
#pragma IgorVersion=6.2
#pragma version=20130725

function averageExperiments(data, dataTimes, expStart, expStop)
 // Given two waves 'expStart' and 'expStop' that represent expt measurement
 // periods, averages wave 'data' for each time range in expStart/expStop.
 //
 // INPUT:
        // wave expStart  = experiment start times (type "dat", in seconds)
        // wave expStop   = experiment stop times
        // wave data      = data to average
        // wave dataTimes = timestamps for wave data
 // OUTPUT:
        // wave data+"_ExpAvg" = average of data for each experiment
        //           [one experiment is the time from expStart[i] to expStop[i], for 0 < i < numpnts(expStart)]
        // wave data+"_ExpAvgSD" = standard error of each average

    wave expStart
    wave expStop
    wave data
    wave dataTimes
   
    variable nExp = numpnts(expStart)
    variable i
    wave dtemp
   
    // sanity checks
    if (!(WaveExists(expStart) && WaveExists(expStop)))
            print "ERROR in averageExperiments(): expStop and expStart do not both exist."
            return -1
    endif
   
    for (i=0; i<nExp; i+=1)
        // throw an error if expStop[i]<expStart[i]
        if (expStop[i] < expStart[i])
            printf "ERROR in averageExperiments(): expStop times cannot precede expStart times."
            printf "Problem at: %s -- %s", secs2dt(expStart[i]), secs2dt(expStop[i])
            return -1
        endif
        // throw a warning if experiments overlap
        if (i>0)
            if (expStart[i] < expStop[i-1]) // if next exp starts before previous ends
                printf "// WARNING: averageExperiments() detected overlapping experiments: "
                printf "(%s -- %s ", secs2dt(expStart[i]), secs2dt(expStop[i])
                printf "overlaps with the next expt at %s). ", secs2dt(expStop[i-1])
                printf "Continuing...\r"
            endif
        endif
    endfor
   
    // prepare waves for saving average and standard error of mean (sem)
    make /o/d/n=(numpnts(expStart)) $(nameOfWave(data)+"_ExpAvg")  /WAVE=expAvg
    make /o/d/n=(numpnts(expStart)) $(nameOfWave(data)+"_ExpAvgSD") /WAVE=expAvgSD
    expAvg  = nan
    expAvgSD = nan
   
    // Extract any points in 'dataTimes' that fall within expStart and expStop,
    // average 'data' at these points into the new waves labelled _expAvg and _expAvgSD
    for (i=0; i<nExp; i+=1)
        // select only times between expStart[i] and expStop[i]
        Extract /FREE/O data, dtemp, dataTimes>expStart[i] && dataTimes<=expStop[i]
        if (numpnts(dtemp)>0)
            WaveStats/Q dtemp                                                  
            // WaveStats sets the variables V_sdev and V_avg
        else
            V_avg = nan
            V_sdev = nan
        endif
        expAvg[i]   = V_avg
        expAvgSD[i] = V_sdev
    endfor
end

static function/s secs2dt(secs)
// Given a datetime in seconds, returns an ISO-8601--formatted datetime.
//      e.g. secs2dt(date2secs(2012,1,31)) returns "2012-01-31 00:00:00"
    variable secs
    if (secs)
        return secs2date(secs,-2) + " " + secs2time(secs,3)
    endif
    return ""
end

July 25, 2013 at 06:26 am - Permalink