epoh老师，这对于您前面是可行的，但是这样似乎不可行。
其实运行结果最后出来的WW矩阵为29*29*41,这里应该是对应最后一个变量得到的41个矩阵。
验证好像也通不过：
Xij=xlsread('E:\youwanhai-software\New\Getis\matlabgdp.xls');
data1=Xij(1,:);
W=WW(:,:,19)
[Xij_bar,GI,ZG,ZZ]=Getis_revised(W,data1);


ZZ =

   86.2954

这里并不与
    1.0000 274.6958    3.9000   19.0000
相等。

不知道可理解正确？

不请楚你的ZZ = 86.2954是怎么算的
延用48楼结果说明:
Indexmat =

    1.0000  274.6959   92.0000
    2.0000  271.1351  146.0000
    3.0000  282.2302  119.0000
    4.0000  294.7071  121.0000

274.6959是这样得来的
也就是max(ZZmat(1,:))=ZZmat(1,92)

k=3:0.01:5;
k=3.91  % s=92
w92=WW(:,:,92)
ZZcfrm=zeros(n1,1);
for p = 1:n1
     [GGI,ZZG,ZZ]=Getis1(w92,Xij(p,:));
     ZZcfrm(p,1)=ZZ;     
   end % end p
ZZcfrm

  274.6959
  270.5487
  281.9495
  294.3711
  300.6167
  311.3500
  312.9155
  313.7616
  315.8861
  320.4374
  323.3554
  326.8118
  330.5298
  335.9608
  340.7885
  346.7261
  347.7216
  346.6318
  343.6299
  341.8789
  337.0366

epoh老师，抱歉，呵呵。
其实61楼的想法是这样的：根据程序运行，对于Xij的第一行数据，其运算得到第19个值(k=3:0.05:5;若k=3:0.01:5，应该是您算出的92)对应的ZZ值是最大的，因此，我想取出对应的权重矩阵WW，按照您说的应该是WW(:,:,19).那么按照道理，我利用这个矩阵WW(:,:,19)和Xij的第一行数据代入原来的Getis1(Wij,Xij)函数中得出的ZZ值也肯定是和indexmat的第一行相等，即为274.6959。而我验证结果却不会，呵呵，不知道我哪步想错了。

Xij=xlsread('E:\youwanhai-software\New\Getis\matlabgdp.xls');
data1=Xij(1,:);
W=WW(:,:,19)
[Xij_bar,GI,ZG,ZZ]=Getis_revised(W,data1);

%WW(:,:,s)=normw(WWij);
   WW(:,:,s)=WWij;

哈哈，epoh老师。那个是我特意在您的程序中修改的。Badinger et al 的Regional convergence in the European Union文中
pages 9/30有这样一句话（附件截图），所以一般进入Getis前，需要先对WW进行row-standardized。而且我验证时都是同步的，即WW要不都未进行标准化，要不都进行标准化。刚才把normw删除了，也是一样得不到
验证。

只要是WW(:,:,s)=WWij;
Getis1,Getis_revised结果都一样

Indexmat =

     1.0000  274.6958    3.9000   19.0000
     2.0000  271.1351    4.4500   30.0000
     3.0000  282.2281    4.2000   25.0000
     4.0000  294.7071    4.2000   25.0000

WW(:,:,s)=WWij;
k=3:0.05:5;;

k=3.9  % s=19
w19=WW(:,:,19)
ZZcfrm=zeros(n1,1);
for p = 1:n1
      [Xij_bar,GGI,ZZG,ZZ]=Getis_revised(w19,Xij(p,:));
      ZZcfrm(p,1)=ZZ;     
    end % end p
ZZcfrm

ZZcfrm =

  274.6958
  270.5437
  281.9275
  294.3465
  300.5823
  311.3228
  312.8908
  313.7366
  315.8609
  320.4114
  323.3264
  326.7815
  330.4998
  335.9301
  340.7580
  346.7105
  347.7064
  346.6194
  343.6193
  341.8733
  337.0315

谢谢epoh老师，可能我中间有些步骤理解错了，我再仔细核对下。

epoh老师，我知道问题在哪了，还是我自己修改的WW(:,:,s)=normw(WWij)这句有问题。
for s =1:length(k)
WWij=zeros(n,m);
for i = 1:n
    for j = 1:m
     if (i~=j)
        WWij(i,j)=exp((-k(s)/d)*dista(i,j));  
     elseif (i==j)
        WWij(i,j)=0;
     end %end if
    end  %end j
      WW(:,:,s)=normw(WWij);end %end i
   for p = 1:n1
     [Xij_bar,GGI,ZZG,ZZ]=Getis1(WWij,Xij(p,:));
     ZZmat(p,s)=ZZ;
     ZZmat1(p,:)=Xij_bar;
   end % end p
end %end s
ZZmat

在上述程序中，虽然我用了WW(:,:,s)=normw(WWij); 但是Getis1中还是调用了未标准化的WWij。
而我之后在单独验证时，又用了标准化的WWij去验证，所以出现结果不一致。

因此，我还想请教您，在这个循环中，如何让Getis1真正调用标准化后的WWij呢？

Getis.m 的计算公式与公式(3.1)不同
你有没其他软件可以比较一下

SpaceStat TUTORIAL page 146/263
http://dae.unizar.es/docencia/regional/spacestat%20Tutorial.pdf

23.2.1 G Statistic
wij(d) stands for an element of the symmetric
(unstandardized) spatial weights matrix for distance d.

epoh老师，您是说GI计算公式不一致吗？即下面一段程序。
GIn=0;
Xn=0;
for m=1:N
for n=1:N
  if n==m
    GIn=GIn;
    Xn=Xn;
  else
   GIn=GIn+Wij(m,n)*Xij(n);
   Xn=Xn+Xij(n);
  end  %end if
end
GI(m)=GIn/Xn;  这一段程序应该和公式（3.1）一致吧。
我刚才用其他软件粗略对比了，GI和EG的结果基本一样。

恩，这里说wij(d)不用标准化。而我上传的几篇文章一般都提到：
W, which is conventionally row-standardized。这个方向的大部分文献都是提到一般为标准化，因为更好解释。

谢谢epoh老师了，呵呵，可能是我太小心眼了，一直盯着文献的方法。

麻烦上传你的数据
和别的软件结果

epoh老师，我利用stata验证了一下，R应该也有函数可以验证。
其中Getis.csv是对应WW，data.csv是对应Xij的第一行，RE是Stata计算结果。
如果方便，还请epoh老师帮我修改对WW 进行row-standardized的程序。

yc=[40.16304.....]
xc=[116.38513....]
Xij=xlsread('matlabgdp.xls');
[n1,m1]=size(Xij);
dista=GeogDistance(yc,xc);
[n,m]=size(dista);
W=row_norm(dista)
k=40:0.05:50;
s=length(k);
ZZmat=zeros(n1,length(k));
ZZmat1=zeros(n1,m1);
WWij=zeros(n,m);
WW=zeros(n,m,s);

epoh老师，谢谢您
不过我觉得先对WWij进行标准化再代入exp（）与在exp（）运行得到的WWij再标准化应该不一样哦。


您核对公式（3.1）有错误吗？

公式(3)是对的
主要是在12楼,我更改过Getis.m
我忘记当时有否改过GI算法
所以才想在其他软件确认一下
这跟STATA结果是一样的.
话说回来,文献上要求row-standardized
但是软件STATA,Getis.m,SpaceStat都没如此要求

yc=[40.16304....]
xc=[116.38513 ....]
result = distance_wm(yc,xc);
Xij=xlsread('matlabgdp.xls');
[n1,m1]=size(Xij);
WWij=result.dw;
     [Xij_bar,GGI,ZZG,ZZ]=Getis_revised(WWij,Xij(1,:));
GGI  %the same as stata

GGI =

    0.5595
    0.5642
    0.5496
    0.5925
    0.3093
    0.4133
    0.2951
    0.0839
    0.5709
    0.6024
    0.5349
    0.5767
    0.5001
    0.6090
    0.6120
    0.6768
    0.7928
    0.6395
    0.2509
    0.3381
    0.1889
    0.2577
    0.4186
    0.1527
    0.5950
    0.1186
    0.0958
    0.4154
    0.0038

epoh老师，是这样的，都是要求不同。
stata在计算时要求必须是元素为0-1且不能标准化的。SpaceStat那部分其实只是计算G指数，并没用涉及到spatial filtering。
在现在主流文献还是要求对公式WWij=exp(-k*dij)计算出来的WWij进行标准化，但是没用说一定，只是提到
“通常”进行标准化。

底下这两种结果是相同的
针对75楼copy错误修改过
请试试Spatial Filtering效果如何
%%%%%case 1
dista=GeogDistance(yc,xc);
[n,m]=size(dista);
W=row_norm(dista)
k=40:0.1:50;
s=length(k);
ZZmat=zeros(n1,length(k));
ZZmat1=zeros(n1,m1);
WWij=zeros(n,m);
WW=zeros(n,m,s);
for s =1:length(k)
WWij=zeros(n,m);
for i = 1:n
    for j = 1:m
     WWij(i,j)=exp(-k(s)*W(i,j));  
     if ( i~=j & WWij(i,j) > 0)
         WWij(i,j)= WWij(i,j);  
     else
         WWij(i,j)=0;
     end %end if
    end  %end j      
      WW(:,:,s)=WWij;
end %end i
   for p = 1:n1
     [Xij_bar,GGI,ZZG,ZZ]=Getis_revised(WWij,Xij(p,:));
     ZZmat(p,s)=ZZ;
     ZZmat1(p,:)=Xij_bar;
   end % end p
end %end s
ZZmat

Indexmat=zeros(n1,4);%3--->4
for i =1:n1
Indexmat(i,1)=i;
index=find(ZZmat(i,:)==max(ZZmat(i,:)));
Indexmat(i,2)=ZZmat(i,index);
Indexmat(i,3)=k(index);
Indexmat(i,4)=index; %Indexmat(i,3)=index;
end
Indexmat
%%%%%%%%%%
%%%%%%%%%case 2

dista=GeogDistance(yc,xc);
d=max(max(dista))   
dista=dista/d
[n,m]=size(dista);
W=row_norm(dista)
k=40:0.1:50;
s=length(k);
ZZmat=zeros(n1,length(k));
ZZmat1=zeros(n1,m1);
WWij=zeros(n,m);
WW=zeros(n,m,s);
for s =1:length(k)
WWij=zeros(n,m);
for i = 1:n
    for j = 1:m
     WWij(i,j)=exp(-k(s)*W(i,j));  
     if ( i~=j & WWij(i,j) > 0)
         WWij(i,j)= WWij(i,j);  
     else
         WWij(i,j)=0;
     end %end if
    end  %end j      
      WW(:,:,s)=WWij;
end %end i
   for p = 1:n1
     [Xij_bar,GGI,ZZG,ZZ]=Getis_revised(WWij,Xij(p,:));
     ZZmat(p,s)=ZZ;
     ZZmat1(p,:)=Xij_bar;
   end % end p
end %end s
ZZmat

Indexmat=zeros(n1,4);%3--->4
for i =1:n1
Indexmat(i,1)=i;
index=find(ZZmat(i,:)==max(ZZmat(i,:)));
Indexmat(i,2)=ZZmat(i,index);
Indexmat(i,3)=k(index);
Indexmat(i,4)=index; %Indexmat(i,3)=index;
end
Indexmat

Indexmat =

    1.0000  275.2452   43.7000   38.0000
    2.0000  272.2491   44.5000   46.0000
    3.0000  283.9725   45.1000   52.0000
    4.0000  296.3250   45.4000   55.0000
    5.0000  303.0755   46.6000   67.0000
    6.0000  314.4640   46.5000   66.0000
    7.0000  316.1042   45.1000   52.0000
    8.0000  316.9693   45.2000   53.0000
    9.0000  319.0838   45.2000   53.0000
   10.0000  323.6744   45.3000   54.0000
   11.0000  326.5613   45.6000   57.0000
   12.0000  330.1633   45.7000   58.0000
   13.0000  333.9460   45.7000   58.0000
   14.0000  339.6418   45.6000   57.0000
   15.0000  344.8175   45.5000   56.0000
   16.0000  350.8192   45.5000   56.0000
   17.0000  351.7772   45.4000   55.0000
   18.0000  349.8736   45.0000   51.0000
   19.0000  345.5357   44.7000   48.0000
   20.0000  343.3499   44.2000   43.0000
   21.0000  338.5056   41.8000   19.0000

非常感谢epoh老师
哈哈，epoh老师，对于标准化，我的看法还是和您的不同，请不要介意啊。
您看我上传的那个Badinger et al的文章，公式顺序应该是：1、先用(3.4)计算出Wij(d)，
这一步未提到需要标准化dij.2、把计算出来的Wij(d)标准化后代入公式(3.1)计算Gi统计量。
您觉得呢？

呵呵,讨论原本就应有不同想法,
假设依你的想法,标准化的话,
可以这样处理
[Xij_bar,GGI,ZZG,ZZ]=Getis_revised(row_norm(WWij),Xij(p,:));
但是抓不到k
k一直往上爬,到1000,3000都是一样的值,无法辨识
你只要作Spatial Filtering 就知道答案了

epoh老师，刚用您79楼的程序试算了一下，结果还是不理想。
其实我有点好奇，为什么计算出来的Xij_bar会和Xij差那么大呢，
而且Xij_bar都比Xij大，就是说过滤后的值比未过滤的更大。
从理论来说，没有规定Xij_bar不能比Xij大，但是从别人文献结果来看，
L=Xij-Xij_bar基本都是有正也有负的。
不明白是不是真是我的数据问题。

可能是数据转换问题
刚又看了下
这篇文献也是要cartesian coordinates
Spatial Filtering with EViews and MATLAB
2.1 Global Spatial Structures
(theta,w) are the cartesian coordinates

是的啊，我看了那么久，怎么都没有注意到这个，非常感谢epoh老师的提醒。
现在的程序调用的是distance_wm,那是根据经纬度计算的 the Great Circle distance，
而不是在cartesian coordinates下的坐标系统。

我记得这个问题答复过你
但是我忘记转换函数
你还记得吗?

epoh老师，您是指从经纬度到cartesian coordinates的转换函数吗？这个好像没有看见过，我不知道。
您还记得大致是什么时候告诉我的吗？

NDSF Utility: LatLong to XY
http://www.whoi.edu/marine/ndsf/ ... om=LatLon&to=XY

epoh老师，太感谢您了，还有这个功能
不过我试了下，怎么都转换不成功哦。比如：

35° 18′ 34″ N, 106° 33′ 7″ W

35.309444, -106.551944

转换为X和Y，得出的结果是NAN,NAN。
不知道这一栏要求填什么内容？Enter the Origin of X/Y coordinate:

你先填0,0试试
另请再看看这个function如何

好的，谢谢您
那个Matlab函数可以运行，结果：
Lat=35.309444;Lon=-106.551944;
[x,y,utmzone] = wgs2utm(Lat,Lon)
x =
3.5891e+005
y =

  3.9085e+006
utmzone =

    13从google资料来看，这里应该为6的。


另外，如果那个网页形式的第二个框都是输入0,0，结果和这个不一样。

我再学习下其他地理软件，确认下结果。

google参考：
http://www.google.com.hk/#hl=zh-CN&newwindow=1&safe=strict&site=&source=hp&q=Bernalillo&btnK=Google+%E6%90%9C%E7%B4%A2&oq=&aq=&aqi=&aql=&gs_l=&bav=on.2,or.r_gc.r_pw.,cf.osb&fp=49354fa1230e3dc5&biw=1440&bih=693