R程序改写为matlab程序

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

程序1-R程序
x<-c(1,2,3)
n<-length(x)
xibar <- (rep(sum(x), n) - x)/(n - 1)
si2 <- ((rep(sum(x^2), n) - x^2)/(n - 1)) - xibar^2
Wi <- sapply(listw$weights, sum)
Wi2=Wi^2

S1i <- sapply(listw$weights, function(x) sum(x^2))

lx <- lag.listw(listw, x, zero.policy = zero.policy)

res <- (lx - Wi * xibar)

res <- res/sqrt(si2 * (((n - 1) * S1i - Wi^2)/(n - 2)))

程序2-matlab程序

x=input('x=')
weights=input('weights=')
n=length(x)
xibar=(sum(x)*ones(n,1)-x)/(n-1)
si2=((sum(x.^2).*ones(n,1)-x.^2)/(n-1))-xibar.^2
Wi=weights(:)'  %这句对应Wi <- sapply(listw$weights, sum)是否正确？
Wi2=Wi.^2
S1i=sum(weights.^2)
lx=Wi*x            %这句如何从lx <- lag.listw(listw, x, zero.policy = zero.policy)改写

res=lx-Wi*xibar

res =res./sqrt(si2 .* (((n - 1) * S1i - Wi.^2)/(n - 2)))

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：MATLAB程序 MATLAB atlab matla Atl matlab

这个程序哪里不对

%输入一个列向量x（n行1列）和矩阵（n行n列）
x=input('x=')
weights=input('weights=')
n=length(x)
xibar=(sum(x)*ones(n,1)-x)/(n-1)
si2=((sum(x.^2).*ones(n,1)-x.^2)/(n-1))-xibar.^2
%Wi=weights(:)'
Wi=sum(weights')
Wi2=Wi.^2
S1i=sum(weights.^2')
lx=Wi*x
res=lx-Wi*xibar
res =res./sqrt(si2 .* (((n - 1).*S1i-Wi2)./(n - 2)))

library(spdep)
example(columbus)
col.listw <- nb2listw(col.gal.nb)
x=seq(1:49)
n<-length(x)
xibar <- (rep(sum(x), n) - x)/(n - 1)
si2 <- ((rep(sum(x^2), n) - x^2)/(n - 1)) - xibar^2
Wi <- sapply(col.listw$weights, sum)
#[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#[37] 1 1 1 1 1 1 1 1 1 1 1 1 1
Wi2=Wi^2
S1i <- sapply(col.listw$weights, function(x) sum(x^2))
# [1] 0.5000000 0.3333333 0.2500000 0.2500000 0.1428571 0.5000000 0.2500000
# [8] 0.1666667 0.1250000 0.2500000 0.2000000 0.1666667 0.2500000 0.1666667
#[15] 0.1666667 0.1428571 0.3333333 0.2500000 0.3333333 0.1000000 0.3333333
#[22] 0.1666667 0.3333333 0.1428571 0.1428571 0.1666667 0.2500000 0.1111111
#[29] 0.1428571 0.2500000 0.5000000 0.2500000 0.2500000 0.2500000 0.1428571
#[36] 0.2000000 0.1666667 0.1666667 0.5000000 0.2000000 0.3333333 0.5000000
#[43] 0.1666667 0.2000000 0.2500000 0.5000000 0.5000000 0.2500000 0.3333333
lx <- lag.listw(col.listw, x)
res <- (lx - Wi * xibar)
res <- res/sqrt(si2 * (((n - 1) * S1i - Wi^2)/(n - 2)))
res
[1]  -2.37332110 -2.90731303 -3.33404263 -3.10096296 -3.52712088 -1.87601018
[7]  -2.00517335 -3.07281364 -2.00272695 -1.21686333 -2.32942618 -2.51955856
[13] -2.18641859 -2.01890149 -1.80019026 -1.63006479 -0.93223336 -1.00013383
[19] -0.80158055 -0.07496224  0.52694104 -0.91849373 -0.25298002 -0.48562024
[25] -0.79300074 -0.66119098  0.11454945  1.21862840  0.95192246  0.41343356
[31]  1.01485374  0.89135794  0.38706277  1.11682531  2.06289781  2.24404180
[37]  1.96425487  2.06215004  1.64678694  1.70734475  1.88352479  1.05182914
[43]  3.05170115  3.00125143  2.90624770  1.32511796  1.63837602  3.08496360
[49]  2.70457746
############

但是function lag.listw()

要转成matlab比较麻烦

因为lag.listw()是调用C function

epoh老师，我是想根据附件中的4个公式计算Gi（d），其中已知x和w（为一个矩阵）。
第（1）个公式是第（2）个公式的标准化
呵呵，我看这个很好实现的，但是程序总是运行不出来

我写的程序：
%输入一个列向量x（n行1列）和矩阵（n行n列）
x=input('x=')
weights=input('weights=')
n=length(x)
xibar=(sum(x)*ones(n,1)-x)/(n-1)
si2=((sum(x.^2).*ones(n,1)-x.^2)/(n-1))-xibar.^2
%Wi=weights(:)'
Wi=sum(weights')
Wi2=Wi.^2
S1i=sum(weights.^2')
lx=Wi*x
res=lx-Wi*xibar
res =res./sqrt(si2 .* (((n - 1).*S1i-Wi2)./(n - 2)))

这里我暂时假定了G与d无关，比较好编写。而且假定了w的对角线全为0
后面可能还要与d相关，所以如果用R比较难实现，因此想转换为matlab

未命名3.jpg & 未命名2.jpg

有没矛盾?两个都是 j not equal to i

Wij(d)是甚么关系?要知道

如果公式对,就有方向了

epoh 发表于 2011-10-30 22:07
未命名3.jpg & 未命名2.jpg 有没矛盾?两个都是 j not equal to i Wij(d)是甚么关系?要知道如果公式对,就有方 ...

epoh老师，未命名（3）与（2）没有矛盾，（3）是（2）的标准化形式。
Wij（d）是一个矩阵，为已知的。
若x为n*1的列向量，Wij（d）就为n*n的矩阵，且x也为已知数。
附件中为参考文章
reference.rar (1.11 MB) 本附件包括：

• reference.pdf

2011-10-30 23:17:03 上传

function [GI,ZG]=Getis(Wij,Xij)

% 计算G统计值：G值是不包含自己的G统计值

% 输入：Wij—空间权值矩阵

% 输入：Xij—研究区域的空间属性数据

% 输出：GI— 空间局部自相关指标-G统计值

% 输出：ZG— 对于GIJ的检验Z值的计算结果

这里有现成程序,请先看看.

局部G统计计算的Matlab程序.pdf

   局部G统计计算的Matlab程序.pdf (256.93 KB)

2011-10-31 09:21:20 上传

epoh老师，谢谢您，我先用这个程序试下
如需修改，而自己不会，再找您

ct_data.mat,distance_wm.m,Getis.m
    Getis.rar (2.92 KB)

2011-10-31 14:46:13 上传

distance_wm.m from
  http://community.wvu.edu/~djl041/matlab.html
%%%%%%run script
load ct_data    % Load Connecticut county latitude and longitude coordinates
                       % 8 counties in all
                       % http://en.wikipedia.org/wiki/List_of_counties_in_Connecticut
result = distance_wm(yc,xc);   % Make sure that you enter latitude and longitude
                                                % in that order
Xij=result.dista                         % This martrix contains the pairwise distances
Wij=result.dw                          % This matrix contains the connectivity structure
                                               % 1 = neighbor, 0 = not a neighbor     
[GI,ZG]=Getis(Wij,Xij)
%%%%%%%%%%
Getis_revised(2011.10.31  21:03)