R程序改写为matlab程序

ywh19860616 发表于 2012-5-12 16:00
epoh老师，两篇主要参考文献上传了，14楼公式我好想就是安装Badinger-Regional convergence in the Europ ...

我刚试了,修改后
k=300,都没问题
只要改两行

dista = GeogDistance(yc,xc);
d=max(max(dista))

WWij(i,j)=exp((-k(s)/d)*dista(i,j));

epoh 发表于 2012-5-12 16:41
我刚试了,修改后
k=300,都没问题
只要改两行

epoh老师，恩，我理解您说的，可是，从40楼上传的文献，作者好像并未提到W需要那样的结构。
with dij denoting the geographical distance between the centres of the regions i and  j.
很多文章直接通过threshold distance d，来确定Wij(d)，即在d范围内，Wij取1，在d范围外，Wij取0.
附件这篇文章即是，也是提出spatial filtering较原始的。里面提到了Getis1里面程序详细步骤。

spatial statistics.rar (1.12 MB) 本附件包括：

• spatial statistics.pdf

2012-5-12 20:27:59 上传

ywh19860616 发表于 2012-5-12 20:25
epoh老师，恩，我理解您说的，可是，从40楼上传的文献，作者好像并未提到W需要那样的结构。
with dij ...

两个d是不同的
第一个 d :threshold distance d
在 distance_wm.m 就用上了
dista = GeogDistance(lat,long);      % Create pairwise distance matrix
temp1 = sort(dista);                 % Sort the distance matrix (1st row will be zeros)
temp2 = temp1(2:end,:);              % Grab the second through nth row
temp3 = min(temp2);                  % Find minimum values of all columns
d = max(temp3);                      % Find largest minimum value and use that as d
% d = user specified distance for cutoff
% e.g. d = 10 would indicate a 10 mile distance cutoff
dw = dista;
[n k] = size(dw);
for i = 1:n
    for j = 1:n
        if (dw(i,j)==0)
            dw(i,j)=0;
        elseif (dw(i,j)<=d);
            dw(i,j)=1;
        else
            dw(i,j)=0;
        end
    end
end
W = sparse(normw(dw));

第二个 d : d=max(max(dista))   %82.0487
主要是用来转换成FRACTIONAL distance matrix
with dij = 1 for the maximum pairwise distance
dista
dista =
         0   26.3976   46.9654   64.7889   62.6534   36.7131   43.1493   30.9348
   26.3976         0   20.6691   38.3926   39.3141   50.1330   25.8122   29.2114
   46.9654   20.6691         0   18.0368   28.1406   67.6762   28.9354   43.4492
   64.7889   38.3926   18.0368         0   24.4042   82.0487   37.9407   56.7230
   62.6534   39.3141   28.1406   24.4042         0   68.3261   22.4538   43.2943
   36.7131   50.1330   67.6762   82.0487   68.3261         0   46.2011   25.5506
   43.1493   25.8122   28.9354   37.9407   22.4538   46.2011         0   20.8906
   30.9348   29.2114   43.4492   56.7230   43.2943   25.5506   20.8906         0
dista/d
ans =
         0    0.3217    0.5724    0.7896    0.7636    0.4475    0.5259    0.3770
    0.3217         0    0.2519    0.4679    0.4792    0.6110    0.3146    0.3560
    0.5724    0.2519         0    0.2198    0.3430    0.8248    0.3527    0.5296
    0.7896    0.4679    0.2198         0    0.2974    1.0000    0.4624    0.6913
    0.7636    0.4792    0.3430    0.2974         0    0.8328    0.2737    0.5277
    0.4475    0.6110    0.8248    1.0000    0.8328         0    0.5631    0.3114
    0.5259    0.3146    0.3527    0.4624    0.2737    0.5631         0    0.2546
    0.3770    0.3560    0.5296    0.6913    0.5277    0.3114    0.2546         0

epoh 发表于 2012-5-12 21:33
两个d是不同的
第一个 d :threshold distance d
在 distance_wm.m 就用上了

是的，epoh老师，我也觉得您现在提到的比较合理。哈哈，那我感觉40楼的文献叙述的有问题，的确提到了
dij denoting the geographical distance between the centres of the regions i and  j.
且没有说明要转换成转换成FRACTIONAL distance matrix

谢谢epoh老师，我再仔细看下原文，然后再请教您。

ywh19860616 发表于 2012-5-12 22:06
是的，epoh老师，我也觉得您现在提到的比较合理。哈哈，那我感觉40楼的文献叙述的有问题，的确提到了
di ...

底下dist_wts.m
我觉得应该对你有帮助
dist_wts.rar (2.06 KB) 本附件包括：

• dist_wts.m

2012-5-12 22:24:02 上传

exponential gravity-model form
    Wij = exp(-a d(ij))

power gravity-model form
    Wij = d(ij)^(-a)

double-power form
    Wij = [1 - (dij/b)^k]^k , dij <= b
                       = 0  ,  dij > b

      where b is a bandwidth for the function, and
      where k is a positive integer power (usually k = 2 or higher)

是的，这个资料很全，谢谢epoh老师

epoh 发表于 2012-5-12 22:23
底下dist_wts.m
我觉得应该对你有帮助

epoh老师，想请教下您，您那里运行出来的结果indexmat如何？我这里还是非常规则。
k=0.001:0.001:0.01;

Indexmat =

    1.0000   22.9703   10.0000
    2.0000   22.4016   10.0000
    3.0000   23.1407   10.0000
    4.0000   24.0234   10.0000
    5.0000   24.4507   10.0000
    6.0000   25.4244   10.0000
    7.0000   25.7352   10.0000
    8.0000   25.7857   10.0000
    9.0000   25.9833   10.0000
   10.0000   26.3326   10.0000
   11.0000   26.4050   10.0000
   12.0000   26.6748   10.0000
   13.0000   27.0254   10.0000
   14.0000   27.4590   10.0000
   15.0000   27.9593   10.0000
   16.0000   28.4113   10.0000
   17.0000   28.4975   10.0000
   18.0000   28.5247   10.0000
   19.0000   28.3209   10.0000
   20.0000   28.4200   10.0000
   21.0000   28.0158   10.0000

就你的数据而言
落点在k=3:0.01:5;
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
    5.0000  301.2764  132.0000
    6.0000  312.0266  132.0000
    7.0000  313.2366  119.0000
    8.0000  314.0910  120.0000
    9.0000  316.2186  120.0000
   10.0000  320.7879  120.0000
   11.0000  323.7873  123.0000
   12.0000  327.2818  124.0000
   13.0000  330.9873  124.0000
   14.0000  336.4309  124.0000
   15.0000  341.2427  123.0000
   16.0000  347.1262  124.0000
   17.0000  348.0859  124.0000
   18.0000  346.6988  103.0000
   19.0000  343.6780  102.0000
   20.0000  341.8913   97.0000

epoh 发表于 2012-5-13 09:43
就你的数据而言
落点在k=3:0.01:5;
Indexmat =

epoh老师，谢谢
您是如何试出来的？这里面有没有规律。
比如要求k属于0到无穷，我们运行不可能这样去设置，所以发现这个是很重要的。

ywh19860616 发表于 2012-5-13 09:47
epoh老师，谢谢
您是如何试出来的？这里面有没有规律。
比如要求k属于0到无穷，我们运行不可能这样去设 ...

这个很简单,先设
k=1:1:100
假设index都是1,表示下回试 0.1:0.01:1

假设index都是100,表示下回试 100:1:200

假设index在10附近,表示下回试 8:0.01:12

两三次之后,就可抓到最适点