Averaging Histogram

This function takes as an input irregularly spaced data (at time points in timewave) and creates evenly spaced data (with x-values at the center of each bin). The data in each bin are averaged, and the standard error of the mean is given, too. It should work with bins that have one or two points in them.

Note that the standard deviation calculated should be used with caution. In reporting the standard error of the mean, the function implicitly assumes that the variation of the data across the bin is small compared to the statistical errors in the data. In the example code below the function, this condition is not met, and the error is really a combination of a deterministic spread and statistical error.

Note that the function requires v6.1 only because it uses the new /FREE flag. If you use an earlier version, just omit the flag. You may then want to add a KillWaves statement to delete those temporary waves.

Function HistAvgSd(datawave,timewave,dt,tmin)       // average all points in interval dt; computes sd, too
Wave datawave,timewave                      // datawave = data taken at unevenly spaced times
Variable dt,tmin                    // dt = bin size, tmin = start time of binned data
Variable n = numpnts(datawave)              // number of points in datawave
Variable tmax = timewave[numpnts(timewave)-1], k=ceil((tmax-tmin)/dt)   // number of bins
String wd_string = NameOfWave(datawave)+ "_avg", ws_string = NameOfWave(datawave) + "_sd"
Make /o/n=(k) \$wd_string=NaN,\$ws_string=NaN             //  avg & sd waves
Make /o/n=(k+1) /FREE indexwave                 // left-bin positions (temp wave)
Wave wd = \$wd_string, ws = \$ws_string
SetScale/P x tmin+dt/2,dt,"", wd, ws                    // x-scaling for average, sd waves
SetScale/P x tmin,dt,"", indexwave                  // x-scaling for index wave
Duplicate /o /FREE datawave datawave1                  // sorted waves (temporary)
Duplicate /o /FREE timewave timewave1
Sort timewave1 timewave1,datawave1                  // sort according to times

indexwave = BinarySearch(timewave1,x)+1                // index vals for left edges of bins
indexwave[k] = n                            // to get all the points in last bin
wd = mean(datawave1,indexwave[p],indexwave[p+1]-1)
wd = (indexwave[p+1]== indexwave[p]) ? NaN : wd             // if no data in bin set to NaN
ws = sqrt(Variance(datawave1,indexwave[p],indexwave[p+1]-1)/(indexwave[p+1]-indexwave[p]))
End

Here is some code to test the function:
make /o/n=20 dtwave,datwave
•dtwave = enoise(1)+1
Integrate dtwave/D=twave                     // irregularly spaced time points
•datwave = p + gnoise(1)
display datwave vs twave
•HistAvgSd(datwave,twave,5,0)
append datwave_avg
ModifyGraph mode=3,marker(datwave)=8,marker(datwave_avg)=19;DelayUpdate
ModifyGraph rgb(datwave_avg)=(1,4,52428);DelayUpdate
ErrorBars datwave_avg Y,wave=(datwave_sd,datwave_sd)
ModifyGraph grid(bottom)=2

Forum

Support

Gallery