R程序改写为matlab程序

epoh 发表于 2012-5-15 20:39
Getis.m 的计算公式与公式(3.1)不同
你有没其他软件可以比较一下

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的结果基本一样。

epoh 发表于 2012-5-15 21:05
SpaceStat TUTORIAL page 146/263
http://dae.unizar.es/docencia/regional/spacestat%20Tutorial.pdf ...

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

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

ywh19860616 发表于 2012-5-15 21:43
恩，这里说wij(d)不用标准化。而我上传的几篇文章一般都提到：
W, which is conventionally row-standar ...

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

epoh 发表于 2012-5-15 21:46
麻烦上传你的数据
和别的软件结果

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


RE.doc (29 KB)

2012-5-15 22:15:06 上传

Getis.rar (274 Bytes) 本附件包括：

• Getis.csv

2012-5-15 22:12:49 上传

data.rar (214 Bytes) 本附件包括：

• data.csv

2012-5-15 22:12:49 上传

ywh19860616 发表于 2012-5-15 22:16
epoh老师，我利用stata验证了一下，R应该也有函数可以验证。
其中Getis.csv是对应WW，data.csv是对应Xij ...

row_norm.rar (364 Bytes) 本附件包括：

• row_norm.m

2012-5-16 10:34:44 上传

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 发表于 2012-5-16 10:33
yc=[40.16304.....]
xc=[116.38513....]
Xij=xlsread('matlabgdp.xls');

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


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

ywh19860616 发表于 2012-5-16 10:48
epoh老师，谢谢您
不过我觉得先对WWij进行标准化再代入exp（）与在exp（）运行得到的WWij再标准化应该 ...

公式(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 发表于 2012-5-16 14:43
公式(3)是对的
主要是在12楼,我更改过Getis.m
我忘记当时有否改过GI算法

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

ywh19860616 发表于 2012-5-16 14:52
epoh老师，是这样的，都是要求不同。
stata在计算时要求必须是元素为0-1且不能标准化的。SpaceStat那部分 ...

底下这两种结果是相同的
针对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 发表于 2012-5-16 20:01
底下这两种结果是相同的
针对75楼copy错误修改过
仔细推敲文献,应该是标准化再代入exp()没错

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