R程序改写为matlab程序

epoh 发表于 2012-5-13 11:24
这个很简单,先设
k=1:1:100
假设index都是1,表示下回试 0.1:0.01:1

谢谢，这个方法很好。
epoh老师，请您抽空再帮我看下一段类似功能的R 程序。这里的主要区别是矩阵生成机制不同，这里我想根据这样一个原则。给定一个门限值d，若d(i,j)>d，d(i,j)=1，若小于1，d(i,j)=0。而我给的不同的d，可以得到不同的权重矩阵。这里为什么选择用R，是因为R里面有一个moranI函数可以计算moran指数。我根据您编写的matlab修改，但是总不能出结果。

yc=as.numeric(t(c(40.16304,39.28829,39.53356,37.56977,44.08654,41.29838,43.6616,47.84179,31.16541,32.98073,29.16594,32.98073,26.07515,27.60815,36.32174,33.87427,30.96772,27.60539,23.35067,23.83548,19.19263,30.61157,26.8115,24.97481,35.20325,37.8208,35.74278,37.26806,41.11203)));
xc=as.numeric(t(c(116.38513,117.30503,116.08708,112.26255,113.89869,122.56549,126.14836,127.72512,121.39731,119.4195,120.02307,119.4195,117.96558,115.69939,118.10671,113.58064,112.23893,111.69403,113.40108,108.78058,109.74489,102.72802,106.87227,101.4963,108.85362,100.93231,96.02154,106.15805,85.20092)));
data<-read.csv('data1.csv')
p<-nrow(data)
LL<-cbind(xc,yc)
dista<-rdist.earth(LL)
d = seq(10,200,10)
s<-length(d)
dw<-dista
n<-ncol(dw)
WW<-matrix(0,n,n,s)
for (p in 1:s){
WW<-matrix(0,n,n)
for (i in 1:n){
    for (j in 1:n) {
        if (dw(i,j)==0)
            dw(i,j)<-0
        else if (dw(i,j)<=d(s))
            dw(i,j)<-1
        else
            dw(i,j)<-0}
}
WW(:,:,s)=dw
moranI(WW,data)   ##data为一列数据，每次计算出来的结果也只有一个数值
}
data1.rar (226 Bytes) 本附件包括：

• data1.csv

2012-5-13 12:12:27 上传

epoh 发表于 2012-5-13 11:24
这个很简单,先设
k=1:1:100
假设index都是1,表示下回试 0.1:0.01:1

epoh老师，若在调用Getis1之前，先对Wij进行row-standardized，是不是修改为：
......
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
      WWij=normw(WWij);
      WW(:,:,s)=WWij;
end %end i
   for p = 1:n1
     [GGI,ZZG,ZZ]=Getis_revised(WWij,Xij(p,:));
     ZZmat(p,s)=ZZ;
   end % end p
end %end s
ZZmat
......
哈哈，还有，如果是这样，这里调试不出最适合点，为啥，epoh老师？

ywh19860616 发表于 2012-5-13 12:10
谢谢，这个方法很好。
epoh老师，请您抽空再帮我看下一段类似功能的R 程序。这里的主要区别是矩阵生成机 ...

library(ape)
library(fields)
yc=as.numeric(t(c(40.16304......)
xc=as.numeric(t(c(116.38513,....)
data<-read.csv('data1.csv')
data=data[,1]
p<-nrow(data)
LL<-cbind(xc,yc)
dista<-rdist.earth(LL)
d = seq(50,1000,50)
s<-length(d)
n<-ncol(dista)
WW= array(NA, dim = c(n,n,s))
moranImat=matrix(NA,s,1);
for (p in 1:s){
dw<-dista
for (i in 1:n){
     for (j in 1:n) {

         if (dw[i,j] <= d[p]) dw[i,j] =0  else
             dw[i,j] = 1
         }  #j
  } #i
WW[,,p]=dw
moranImat[p,1]=Moran.I(data,dw)$observed   ##data为一列数据，每次计算出来的结果也只有一个数值
} #p

WW
moranImat
             [,1]
[1,] -0.03562998
[2,] -0.03505226
[3,] -0.03505226
[4,] -0.03815054
[5,] -0.05173708
[6,] -0.05570482
[7,] -0.06194603
[8,] -0.05337594
[9,] -0.06442231
[10,] -0.05243863
[11,] -0.05588271
[12,] -0.06776555
[13,] -0.07897514
[14,] -0.07869158
[15,] -0.10539806
[16,] -0.12439059
[17,] -0.12607884
[18,] -0.13276167
[19,] -0.15463911
[20,] -0.19021641

epoh 发表于 2012-5-13 16:18
library(ape)
library(fields)
yc=as.numeric(t(c(40.16304......)

library(ape)
library(fields)
yc=as.numeric(t(c(40.16304......)
xc=as.numeric(t(c(116.38513,....)
data<-read.csv('data1.csv')
data=data[,1]
p<-nrow(data)
LL<-cbind(xc,yc)
dista<-rdist.earth(LL)
d = seq(50,1000,50)
s<-length(d)
n<-ncol(dista)
WW= array(NA, dim = c(n,n,s))
moranImat=matrix(NA,s,4);
for (p in 1:s){
dw<-dista
for (i in 1:n){
     for (j in 1:n) {

         if (dw[i,j] <= d[p]) dw[i,j] =0  else
             dw[i,j] = 1
         }  #j
  } #i
WW[,,p]=dw
moranImat[p,4]=Moran.I(data,dw)   ##data为一列数据，每次计算出来的结果也只有一个数值
} #p

epoh老师，我还想请教一下，这句命令data=data[,1]的作用是什么？我试了下，好像还是和原来data一样，只是排列形式变了。

还有，我想把Moran.I函数得到的4个值都输出，我对应改了几句，您看哪里错误？

ywh19860616 发表于 2012-5-13 16:42
library(ape)
library(fields)
yc=as.numeric(t(c(40.16304......)

data<-read.csv('moran.csv')   # "data.frame"
data=data[,1]  #"numeric"
换成"numeric"是为了配合function Moran.I(x, weight)
x : a numeric vector.的要求

data<-read.csv('data1.csv')   # "data.frame"
data=data[,1]  #"numeric"
p<-nrow(data)
LL<-cbind(xc,yc)
dista<-rdist.earth(LL)
d = seq(50,1000,50)
s<-length(d)
n<-ncol(dista)
WW= array(NA, dim = c(n,n,s))
moranImat=matrix(NA,s,4);
for (p in 1:s){
dw<-dista
for (i in 1:n){
     for (j in 1:n) {

         if (dw[i,j] <= d[p]) dw[i,j] =0  else
             dw[i,j] = 1
         }  #j
  } #i
WW[,,p]=dw
moranImat[p,1]=Moran.I(data,dw)$observed  
moranImat[p,2]=Moran.I(data,dw)$expected   
moranImat[p,3]=Moran.I(data,dw)$sd
moranImat[p,4]=Moran.I(data,dw)$p.value   
} #p
colnames(moranImat) = c("observed","expected","sd", "p.value")
print(moranImat)

         observed    expected          sd    p.value
[1,] -0.03562998 -0.03571429 0.003497052 0.98076591
[2,] -0.03505226 -0.03571429 0.004949200 0.89358897
[3,] -0.03505226 -0.03571429 0.004949200 0.89358897
[4,] -0.03815054 -0.03571429 0.007036459 0.72916722
[5,] -0.05173708 -0.03571429 0.011300924 0.15624031
[6,] -0.05570482 -0.03571429 0.014758350 0.17556940
[7,] -0.06194603 -0.03571429 0.018383120 0.15359511
[8,] -0.05337594 -0.03571429 0.019858576 0.37380342
[9,] -0.06442231 -0.03571429 0.022407499 0.20013071
[10,] -0.05243863 -0.03571429 0.026124726 0.52206003
[11,] -0.05588271 -0.03571429 0.029924213 0.50032175
[12,] -0.06776555 -0.03571429 0.033569987 0.33969935
[13,] -0.07897514 -0.03571429 0.036233464 0.23249847
[14,] -0.07869158 -0.03571429 0.039075959 0.27140198
[15,] -0.10539806 -0.03571429 0.042226614 0.09889518
[16,] -0.12439059 -0.03571429 0.045660253 0.05212624
[17,] -0.12607884 -0.03571429 0.050839595 0.07549515
[18,] -0.13276167 -0.03571429 0.054816061 0.07665742
[19,] -0.15463911 -0.03571429 0.059795454 0.04671659
[20,] -0.19021641 -0.03571429 0.064245291 0.01617784

非常感谢您，现在明白了。

epoh 发表于 2012-5-11 09:54
index   1              2              3             4             5             6              7 ...

epoh老师，谢谢您的解答
我还有一个问题，现在已经运行得到：
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
    5.0000  301.2758    4.3000   27.0000
    6.0000  312.0261    4.3000   27.0000
    7.0000  313.2352    4.2000   25.0000
    8.0000  314.0902    4.2000   25.0000
    9.0000  316.2178    4.2000   25.0000
   10.0000  320.7876    4.2000   25.0000
   11.0000  323.7853    4.2000   25.0000
   12.0000  327.2809    4.2500   26.0000
   13.0000  330.9854    4.2500   26.0000
   14.0000  336.4292    4.2500   26.0000
   15.0000  341.2405    4.2000   25.0000
   16.0000  347.1244    4.2500   26.0000
   17.0000  348.0834    4.2500   26.0000
   18.0000  346.6960    4.0000   21.0000
   19.0000  343.6778    4.0000   21.0000
   20.0000  341.8908    3.9500   20.0000
   21.0000  337.0468    3.9500   20.0000
如何根据Indexmat的第一列返回求WW？现在运行得到的WW为29*29*41，应该是最后一个变量的，那么根据WW(:,:19)(比如现在想取第一个19对应的)就会有问题。

ywh19860616 发表于 2012-5-15 13:49
epoh老师，谢谢您的解答
我还有一个问题，现在已经运行得到：
Indexmat =

现在想取第一个19对应的 ??
不太清楚你的意思

epoh 发表于 2012-5-15 14:02
现在想取第一个19对应的 ??
不太清楚你的意思

嗯，epoh老师，上面indexmat的21行分别是对应数据文件的21个变量所得出的最优值，我现在想根据这里的最优值返回确定对应的WW矩阵。

ywh19860616 发表于 2012-5-15 14:10
嗯，epoh老师，上面indexmat的21行分别是对应数据文件的21个变量所得出的最优值，我现在想根据这里的最优 ...

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
    5.0000  301.2758    4.3000   27.0000
   .........

WW(:,:,19)