R程序改写为matlab程序

ywh19860616 发表于 2012-5-17 10:40
好的，谢谢您
那个Matlab函数可以运行，结果：
Lat=35.309444;Lon=-106.551944;

Universal Transverse Mercator coordinate system
http://www.cibolasar.org/UTMConverter.shtml

epoh 发表于 2012-5-17 10:56
Universal Transverse Mercator coordinate system
http://www.cibolasar.org/UTMConverter.shtml

谢谢epoh老师，很好的信息。
研究过再向您反馈。

ywh19860616 发表于 2012-5-17 12:05
谢谢epoh老师，很好的信息。
研究过再向您反馈。

distmat.rar (1.54 KB) 本附件包括：

• distmat.m

2012-5-17 18:31:43 上传

82楼你说,Xij_bar过滤后的值比未过滤的更大
Xij_bar=Xij*EG./GI,这个公式麻烦你再确认一下
我看EG,GI两者蛮接近的
底下程序你再试试
%%%%%%
[x,y,utmzone] = wgs2utm(yc,xc);
xy=[x y];
dmat=distmat(xy);
W=row_norm(dmat);
Xij=xlsread('matlabgdp.xls');
[n1,m1]=size(Xij);
[n,m]=size(dmat);
%k=1:1:100;
k=65:0.1:85;
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  240.8070   82.9000  180.0000
    2.0000  233.5120   77.2000  123.0000
    3.0000  253.4746   77.5000  126.0000
    4.0000  275.7307   75.0000  101.0000
    5.0000  297.0071   73.1000   82.0000
    6.0000  308.8611   72.6000   77.0000
    7.0000  304.8666   73.1000   82.0000
    8.0000  304.2774   72.8000   79.0000
    9.0000  302.5375   72.1000   72.0000
   10.0000  300.6290   71.9000   70.0000
   11.0000  292.0350   71.8000   69.0000
   12.0000  289.9479   71.5000   66.0000
   13.0000  290.0635   70.6000   57.0000
   14.0000  288.0453   70.1000   52.0000
   15.0000  286.0138   69.5000   46.0000
   16.0000  276.6872   71.9000   70.0000
   17.0000  277.0790   66.9000   20.0000
   18.0000  280.4675   66.9000   20.0000
   19.0000  271.0008   68.9000   40.0000
   20.0000  276.7987   69.5000   46.0000
   21.0000  277.7074   69.9000   50.0000

epoh 发表于 2012-5-17 18:30
82楼你说,Xij_bar过滤后的值比未过滤的更大
Xij_bar=Xij*EG./GI,这个公式麻烦你再确认一下
我看EG,GI ...

epoh老师，谢谢您一直热心帮忙。
Xij_bar=Xij*EG./GI 这句程序是我加的，原程序没有，
依据我上传的Badinger那篇文章pages 10/30的 公式（3.3），请您帮忙确认下。
恩，真心觉得这个程序已经非常完美了，依靠wgs2utm和dismat函数，非常漂亮。
对于这两个函数运算出来的各省之间距离也和实际粗略对照了，也是合理的。
然而，针对我的数据，结果还是不理想，现在更大可能是数据自身问题了。

ywh19860616 发表于 2012-5-17 19:42
epoh老师，谢谢您一直热心帮忙。
Xij_bar=Xij*EG./GI 这句程序是我加的，原程序没有，
依据我上传的Bad ...

你有没范例带数据,可以跑一次程序兼除错?
公式(3.2)有没用上?

epoh 发表于 2012-5-17 19:59
你有没范例带数据,可以跑一次程序兼除错?
公式(3.2)有没用上?

epoh老师，在Getis1中，应该没有用到公式(3.2),公式(3.1),(3.3),
(3.4),(3.5)都用上了，其中公式(3.5)就是Getis1程序的最后一句，
也是自己加上的。

这个没有范例数据，实现spatial filtering这个程序基本都是您编写的，
目前好像都没有公布出来的，只有我在短信息中和您提到的那个PPA
软件可以实现这个功能。

我验证程序是这样的：根据Matlab运行结果计算出来的WWij、Xij_bar（即过滤后变量）
以及原始数据Xij，在R中用Moran.I函数进行验证。按照理论，呵呵，或者最好的
结果应该是Xij存在spatial autocorrelation，而Xij_bar不存在。可是现在的结果
都是Xij不存在，反而Xij_bar存在了。

ywh19860616 发表于 2012-5-17 20:14
epoh老师，在Getis1中，应该没有用到公式(3.2),公式(3.1),(3.3),
(3.4),(3.5)都用上了，其中公式(3.5)就 ...

Spatial Filtering with EViews and MATLAB
All programs are available on request by email from the authors
你去函要看看哈哈

epoh 发表于 2012-5-17 22:16
请注意看
Spatial Filtering with EViews and MATLAB
5.2 The Filtering Process

epoh老师，这个应该不会有错误，因为是计算Local Gi，
就是对每个i都会计算出一个值。您让我验证时，我用
Stata也是给出这样的结果。
在Getis1中的这语句 Xij_bar=Xij*EG./GI 应该不会有
错吧？

ywh19860616 发表于 2012-5-17 22:47
epoh老师，这个应该不会有错误，因为是计算Local Gi，
就是对每个i都会计算出一个值。您让我验证时，我用 ...

在Getis1中的这语句 Xij_bar=Xij*EG./GI 应该不会有错吧？
你很care,但你说中了
请改成Xij_bar=diag(EG*Xij)./GI;   
改后以底下而言,应该接近Xij
请注意短信息
Indexmat =

    1.0000  240.8070   82.9000  180.0000

  WWij=WW(:,:,180)
  Xij_bar=ZZmat1(1,:,180)
    1.0e+004 *

    0.0551    0.0340    0.1099    0.0490    0.0346    0.1344    0.0473    0.0835   
    0.0610    0.1200    0.0778    0.0455    0.0501    0.0378    0.1875    0.0731   
    0.0658    0.0727    0.1730    0.0428    0.0093    0.1028    0.0235    0.0397   
    0.0434    0.0263    0.0071    0.0065    0.0275

  Xij(1,:)
    500.82    310.95        896.33        429.27        319.3        1062.7        425.28        715.2        
    781.66    1416.5        904.69           658        522.28        428.62        1511.19        934.65        
    824.38    744.44      1559.028     449.06        102.42      1190.77        260.14        451.67        
     404.3     242.8        69.94            64.84        261.44

epoh 发表于 2012-5-18 09:43
在Getis1中的这语句 Xij_bar=Xij*EG./GI 应该不会有错吧？
你很care,但你说中了
请改成Xij_bar=diag( ...

好的，谢谢epoh老师。
Xij_bar=EG*Xij./GI;   这句程序我是按照公式写的，
没有太注意内部结构，哎，自己添了一句都错误啊。
谢谢epoh老师指正。