grads异常场T检验求助

afan1996

参考了论坛里http://bbs.06climate.com/forum.php?mod=viewthread&tid=18475&extra=&page=1这个帖子。做了前秋季节（son（-1））高度异常场的合成分析。
分别选了EP年6个年份和CP年6个年份。1981-2010年均值作为气候态。
代码如下
'reinit''sdfopen f:\ice\hgt.mon.mean.nc''set lev 500''set lon -180 180''set lat -90 30'
************son-1*******************climatology************i = 21'ss=0''di=0'while(i<=368)    'set t 'i''    'ss=ss+ave(hgt,t='i',t='i+2')'    'di=di+1'    i=i+12endwhile'define cli=ss/di'
'set lev 500''set lon -180 180''set lat -90 30'*************ep***************'set t 45''aae1=ave(hgt,t=45,t=47)-cli''set t 93''aae2=ave(hgt,t=93,t=95)-cli''set t 153''aae3=ave(hgt,t=153,t=155)-cli''set t 225''aae4=ave(hgt,t=225,t=227)-cli''set t 333''aae5=ave(hgt,t=333,t=335)-cli''set t 441''aae6=ave(hgt,t=441,t=443)-cli''define aae=(aae1+aae2+aae3+aae4+aae5+aae6)/6.0'**************cp*************************'set t 105''aac1=ave(hgt,t=105,t=107)-cli''set t 189''aac2=ave(hgt,t=189,t=191)-cli''set t 285''aac3=ave(hgt,t=285,t=287)-cli''set t 309''aac4=ave(hgt,t=309,t=311)-cli''set t 369''aac5=ave(hgt,t=369,t=371)-cli''set t 429''aac6=ave(hgt,t=429,t=431)-cli''define aac=(aac1+aac2+aac3+aac4+aac5+aac6)/6.0'
********TTEST*******************************'define fce=(pow(aae1-aae,2)+pow(aae2-aae,2)+pow(aae3-aae,2)+pow(aae4-aae,2)+pow(aae5-aae,2)+pow(aae6-aae,2))/6.0''define fene=sqrt(fce)''define testep=((aae)/fene)*sqrt(6)''define tte=abs(testep)-2.571'

'define fcc=(pow(aac1-aac,2)+pow(aac2-aac,2)+pow(aac3-aac,2)+pow(aac4-aac,2)+pow(aac5-aac,2)+pow(aac6-aac,2))/6.0''define fenc=sqrt(fcc)''define testcp=((aac)/fenc)*sqrt(6)''define ttc=abs(testcp)-2.571'
**********draw************************************************************'set grid off''set grads off' 'run page q1''set gxout shaded''d tte''cbar''set gxout contour''set cthick 2''set ccolor 1''d aae''set ccolor 1''set cthick 6''d cli''draw xlab ep-son(-1)''run page q2''set gxout shaded''d ttc''cbar''set gxout contour''set ccolor 1''set cthick 2''d aac''set ccolor 1''set cthick 6''d cli''draw xlab cp-son(-1)''gxprint C:\Users\zhuofan\Desktop\hgt.eps eps white';

不知道为什么，出来的t检验值比较大。救救孩子吧。。。。。。T T。。。。。。。

afan1996

'reinit'
'sdfopen f:\ice\hgt.mon.mean.nc'
'set lev 500'
'set lon -180 180'
'set lat -90 30'
************son-1**********
*********climatology************
i = 21
'ss=0'
'di=0'
while(i<=368)
   'set t 'i''  
    'ss=ss+ave(hgt,t='i',t='i+2')'
    'di=di+1'
    i=i+12
endwhile
'define cli=ss/di'
'set lev 500'
'set lon -180 180'
'set lat -90 30'
*************ep***************
'set t 45'
'aae1=ave(hgt,t=45,t=47)-cli'
'set t 93'
'aae2=ave(hgt,t=93,t=95)-cli'
'set t 153'
'aae3=ave(hgt,t=153,t=155)-cli'
'set t 225'
'aae4=ave(hgt,t=225,t=227)-cli'
'set t 333'
'aae5=ave(hgt,t=333,t=335)-cli'
'set t 441'
'aae6=ave(hgt,t=441,t=443)-cli'
'define aae=(aae1+aae2+aae3+aae4+aae5+aae6)/6.0'
**************cp*************************
'set t 105'
'aac1=ave(hgt,t=105,t=107)-cli'
'set t 189'
'aac2=ave(hgt,t=189,t=191)-cli'
'set t 285'
'aac3=ave(hgt,t=285,t=287)-cli'
'set t 309'
'aac4=ave(hgt,t=309,t=311)-cli'
'set t 369'
'aac5=ave(hgt,t=369,t=371)-cli'
'set t 429'
'aac6=ave(hgt,t=429,t=431)-cli'
'define aac=(aac1+aac2+aac3+aac4+aac5+aac6)/6.0'
********TTEST*******************************
'define fce=(pow(aae1-aae,2)+pow(aae2-aae,2)+pow(aae3-aae,2)+pow(aae4-aae,2)+pow(aae5-aae,2)+pow(aae6-aae,2))/6.0'
'define fene=sqrt(fce)'
'define testep=((aae)/fene)*sqrt(6)'
'define tte=abs(testep)-2.571'

'define fcc=(pow(aac1-aac,2)+pow(aac2-aac,2)+pow(aac3-aac,2)+pow(aac4-aac,2)+pow(aac5-aac,2)+pow(aac6-aac,2))/6.0'
'define fenc=sqrt(fcc)'
'define testcp=((aac)/fenc)*sqrt(6)'
'define ttc=abs(testcp)-2.571'
**********draw************************************************************
'set grid off'
'set grads off'
'run page q1'
'set gxout shaded'
'd tte'
'cbar'
'set gxout contour'
'set cthick 2'
'set ccolor 1'
'd aae'
'set ccolor 1'
'set cthick 6'
'd cli'
'draw xlab ep-son(-1)'
'run page q2'
'set gxout shaded'
'd ttc'
'cbar'
'set gxout contour'
'set ccolor 1'
'set cthick 2'
'd aac'
'set ccolor 1'
'set cthick 6'
'd cli'
'draw xlab cp-son(-1)'
'gxprint C:\Users\zhuofan\Desktop\hgt.eps eps white';

第一次发帖，不晓得怎么传程序，这格式怎么不对呀，吞了回车吗，希望好心人顺便指点一下

橙那个心儿

楼主解决了吗，我也遇到了这个问题

afan1996

橙那个心儿 发表于 2020-3-21 17:07
楼主解决了吗，我也遇到了这个问题

'define ccha=sqrt((pow(aac1-aac,2)+pow(aac2-aac,2)+pow(aac3-aac,2)+pow(aac4-aac,2)+pow(aac5-aac,2)+pow(aac6-aac,2))/6)'
'define testcp=aac/(ccha/sqrt(5))'

改成这样就可以啦