* ttest.equal_variances.gs

* Description:
* This script plots t-test values and shades areas that are significant at 
* a specified level. This script assumes that the variance of the two samples 
* is the same. If you are uncertain about this assumption, you should
* use a f-test to determine if you can make this assumption.
 
* For step-by-step instructions on how to run this script, 
* see http://ccr.meteor.wisc.edu/model/grads/grads_statsig.html

****************************************************************************
* This script needs to be used in concert with the web-page specified below.
* This site allows you to compute cut-off t-values, meaning that a t-value
* BELOW the cutoff at a given grid cell means that the difference between the
* two means is NOT significant and a t-test score ABOVE the value means that
* the difference IS significant. 
* I.e., if  ttest (grid cell) > ttest (cutoff), then significant
* If ttest(grid cell) < ttest (cutoff), then NOT significant
*
* TO CALCULATE CUTOFF T-VALUE, goto: 
* http://members.aol.com/johnp71/pdfs.html
* USE "Student t"
* ENTER 'n' DEGREES OF FREEDOM (= N1+N2-2) AND SIGNIFICANCE LEVEL 'p' 
* 0.01 = 99% LEVEL, 0.05 = 95% LEVEL, 0.10 = 90% LEVEL
* THEN CLICK 'Calc t' BUTTON. The cutoff 't' value will appear in the box
* in the 't' box. 
*****************************************************************************

* OPEN DATA FILES (PUT PATH AND FILE NAMES HERE) 
sdfopen /grove4/selin/foam_indyr/ha.F1.6k.djfi486605.TS1.nc 
sdfopen /grove4/selin/foam_indyr/ha.F1.0k.djfi486605.TS1.nc

* ENABLE OUTPUT FILE
enable print sig

set lon -180 180

* ENTER THE SAMPLE SIZES OF YOUR TWO SAMPLES BELOW 
*n1, n2 are number of observations for each experiment 
*degrees of freedom = n1 + n2 -2
n1 = 120 
n2 = 120 

* In the formulas for x1, x2, s1, and s2, you must set the time 
* increment for each of the formulas according to the sample size. 
* x = average, s = standard deviation
x1 = ave(ts1.1,t=1,t=120) 
s1 = sqrt(ave(pow(ts1.1-x1,2),t=1,t=120)*(n1/(n1-1)))
x2 = ave(ts1.2,t=1,t=120)
s2 = sqrt(ave(pow(ts1.2-x2,2),t=1,t=120)*(n2/(n2-1)))

* compute t test statistic
ssqr = ( (n1-1)*s1*s1 + (n2-1)*s2*s2 )/(n1 + n2 - 2)
denom = sqrt(ssqr*((1.0/n1) + (1.0/n2)))
num = x1 - x2
ttest = num/denom

* Take absolute value for a two-tailed test
ttest = abs(ttest)

* Use url above to calculate the  cutoff t.
* In this case, at the .01 (99%) significance level and df=238, 
* the cutoff t value is 2.597 

set grid off
set grads off
set gxout grfill

* set clevs equal to the cutoff t so that all grid cells with values above 
* the cutoff-t will be shaded grey (they are signif)
* and plot the significant areas
set clevs  2.597 
set ccols 0 15
d ttest

* run the color bar script
run cbar.gs

* contour difference of two runs over significance
set gxout contour
set cint 0.25
d x1-x2
draw title FOAM 6K-0K TS1 DJF Delta (contour), Sig at 99% (shaded)
print