* ttest.unequal_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


* You will need to download the xdfopen files in order to run this script as written.
* Since the time indices of the two files (6kFix and 0k) do not match, we used
* the xdfopen option to open our files. Links to these xdfopen files can be found
* at 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 (= MIN(n1-1,n2-1)) 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)
xdfopen xdf6kFixdjf
xdfopen xdf0kdjf

* ENABLE OUTPUT FILE
enable print foam.6kfix-0k.ts1.djf.uneqvar. 

set lon -180 180

* ENTER THE SAMPLE SIZES OF YOUR TWO SAMPLES BELOW
*n1, n2 are number of observations for each experiment
n2 = 120 
n1 = 10 

* x = ave, s = std dev
* 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

x2 = ave(ts1.2,t=1,t=120)
s2 = sqrt(ave(pow(ts1.2-x2,2),t=1,t=120)*(n2/(n2-1)))
x1 = ave(ts1.1,t=1,t=10)
s1 = sqrt(ave(pow(ts1.1-x1,2),t=1,t=10)*(n1/(n1-1)))

* compute t test statistic
denom = sqrt( (pow(s1,2)/n1) + (pow(s2,2)/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 .05 (95%) significance level and df=9,
* the cutoff t value is 2.262

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.262 
set ccols 0 15
d ttest
run cbar.gs

set gxout contour
set cint 0.25
d x1-x2
draw title \FOAM Insolation Only 6kFix-0k TS1 DJF Delta (contour), Sig at 95% (shaded)\
print