关于用NCL做t检验的请教

MacTaylor

请教各位一个关于t检验的问题，我要做一个气象要素5个特征年的平均值减去另外5个特征年平均值的显著性检验，看了书，应该用两组样本平均值差异的显著性检验，然后查了ncl的ttest函数（http://www.ncl.ucar.edu/Document/Functions/Built-in/ttest.shtml）其中例子4符合我的要求，因此我依葫芦画瓢写的代码如下：

begin

    f1 = addfile("hgt_ERA40.nc","r")
      system("/bin/rm "+"test_hgt.nc")
      fo      = addfile("test_hgt.nc","c")          ; create output file
      
      
    hgt = short2flt(f1->z)
    hgt = hgt/9.80665
    printVarSummary(hgt)
    hgts = month_to_season(hgt,"JJA")
    hgts!0 = "time"
    hgts!1 = "lev"
    hgts!2 = "lat"
    hgts!3 = "lon"
    hgts@long_name  = "Geopotential Heigh"
    hgts@units      = "gpm"
    printVarSummary(hgts)
    delete(hgt)
   
    hgtt = hgts(time|0,lev|:,lat|:,lon|:)
    prob = hgtt(0,:,:,:)
    hhgt = hgtt(0,:,:,:)
    lhgt = hgtt(0,:,:,:)
   
    hhgt = (hgts(0,:,:,:)+hgts(23,:,:,:)+hgts(29,:,:,:)+hgts(38,:,:,:)+hgts(40,:,:,:))/5
    lhgt = (hgts(5,:,:,:)+hgts(18,:,:,:)+hgts(26,:,:,:)+hgts(34,:,:,:)+hgts(35,:,:,:))/5

    hhgttmp = hhgt(lev|:,lat|:,lon|:,time|:)
    lhgttmp = lhgt(lev|:,lat|:,lon|:,time|:)
    hhgtAve = dim_avg (hhgttmp)
    lhgtAve = dim_avg (lhgttmp)
    hhgtVar = dim_variance (hhgttmp)
    lhgtVar = dim_variance (lhgttmp)
    sigr = 0.05
    hhgtEqv = equiv_sample_size (hhgttmp, sigr, 0)
    lhgtEqv = equiv_sample_size (lhgttmp, sigr, 0)
    hhgtN   = wgt_areaave (hhgtEqv, wgty, 1., 0)
    lhgtN   = wgt_areaave (lhgtEqv, wgty, 1., 0)
    iflag   = False
    prob = ttest(hhgtAve,hhgtVar,hhgtN, lhgtAve,lhgtVar,lhgtN, iflag, False)
    printVarSummary(prob)
    hgtt(0,:,:,:) = prob(:,:,:)
    delete(prob)
   
;    output the correlation coeffiect and linear regression
    filedimdef(fo,"time",-1,True)

    fo->hgt   = hgtt

end

但是错误提示如下：fatal:Number of subscripts on rhs do not match number of dimensions of variable,(4) Subscripts used, (3) Subscripts expected
fatal:Execute: Error occurred at or near line “红色部分”
然后我意识到我这5年平均的相减时间等于是t=1，怀疑是这里有问题，因此将蓝色部分改为了：
    hhgt = hgts(time|0:4,lev|:,lat|:,lon|:)
    lhgt = hgts(time|0:4,lev|:,lat|:,lon|:)
    hhgt(0,:,:,:) = hgts(0,:,:,:)
    hhgt(1,:,:,:) = hgts(23,:,:,:)
    hhgt(2,:,:,:) = hgts(29,:,:,:)
    hhgt(3,:,:,:) = hgts(38,:,:,:)
    hhgt(4,:,:,:) = hgts(40,:,:,:)
    lhgt(0,:,:,:) = hgts(5,:,:,:)
    lhgt(1,:,:,:) = hgts(18,:,:,:)
    lhgt(2,:,:,:) = hgts(26,:,:,:)
    lhgt(3,:,:,:) = hgts(34,:,:,:)
    lhgt(4,:,:,:) = hgts(35,:,:,:)
然后又提示：fatal:Variable (wgty) is undefined，在hhgtN   = wgt_areaave (hhgtEqv, wgty, 1., 0)这个语句中，再次查了wgt_areaave函数，对于参数wgty和wgtx的说明没怎么看懂，貌似是一个面积权重的东东，不知道在这里应该分别设置成多少，于是根据感觉，将hhgtN和lhgtN的部分注销，将prob = ttest(hhgtAve,hhgtVar,hhgtN, lhgtAve,lhgtVar,lhgtN, iflag, False)改为prob = ttest(hhgtAve,hhgtVar,5, lhgtAve,lhgtVar,5, iflag, False)，于是运行成功，然后画出来的图很怪，我知道不理解这些出现的问题是不行的，结果没有把握，因此特来请教，谁能给我讲讲这里面出现的种种问题啊？非常非常感谢了！

酱紫哈

subtropical 发表于 2016-3-18 10:27
prob = ttest(hhgtAve,hhgtVar,hhgtN, lhgtAve,lhgtVar,lhgtN, iflag, False)得出的这个Prob怎么处理，是直 ...

我也很想问这个出来之后怎么搞

杜晓瑞

？这个我也不懂爱。。

MacTaylor

额。。。谁能帮帮俺呢

wsdeqx

好东西，指的学习

圆圆石头

这个问题我也没看懂，楼主你搞明白了么

下个，路口

不知道楼主的问题解决了木有？大家讨论下呗

深韩

rad=4.0*atan(1.0)/180.0
clat=cos(rad*lat)
hhgtN   = wgt_areaave (hhgtEqv,clat, 1., 0)
这里的wgt我想应该纬度的一个权重系数，我检验的时候是写成了clat，这样是不是对的？

subtropical

想问楼主解决了吗，我刚开始学ncl用到t检验，但是也不懂里面的过程和原理，不会用啊{:eb303:}

subtropical

prob = ttest(hhgtAve,hhgtVar,hhgtN, lhgtAve,lhgtVar,lhgtN, iflag, False)得出的这个Prob怎么处理，是直接画图吗，不懂啊、、、

梦想没有终点

楼主解决了吗，我也遇到了相同的问题