stata求矩阵自然对数

XIEShichen

4639

收藏 2013-11-18

两个相同行列数的矩阵 A, B其中对应的元素求(a-b)/(lna-lnb)

不知道在stata能够实现么，是不是需要自己编程，然后调用
多谢了

扫码加我拉你入群

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

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

全部回复

bbs0805

2013-11-18 08:42:51

mata
(a:-b):/(ln(a):-ln(b))
end

扫码加我拉你入群

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

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

XIEShichen

2013-11-18 09:18:24

bbs0805 发表于 2013-11-18 08:42
mata
(a:-b):/(ln(a):-ln(b))
end

非常感谢啊

如果我想对A，B里面的元素进行判断呢
比如：

cap program drop LMDI

program define LMDI

version 12

syntax namelist

local i = 1

local j = 1

local r = rowsof()

local c = colsof()

while `i'<=r{

while `j'<=c{

`a0'= Mat1[i,j]

`a1'= Mat2[i,j]

if `a0'>0 & `a1'>0{

`L'=(`a1'-`a0')/ln(`a1'/`a0')

}

else if `a0'>0 & `a1'=0{

`L'=(`a1'-`a0')/ln(1E-200/`a0')

}

else if `a0'>0 & `a1'<0{

`L'=(`a1'-`a0')/(ln(-`a1'/1E-200)+ln(1E-200/`a0'))

}

else if `a0'=0 & `a1'>0{

`L'=(`a1'-`a0')/ln(`a1'/1E-200)

}

else if `a0'=0 & `a1'=0{

`L'=0

}

else if `a0'=0 & `a1'<0{

`L'=(`a1'-`a0')/(ln(-`a1'/1E-200))

}

else if `a0'<0 & `a1'>0{

`L'=(`a1'-`a0')/(ln(`a1'/1E-200)+ln(-1E-200/`a0'))

}

else if `a0'<0 & `a1'=0{

`L'=(`a1'-`a0')/ln(-1E-200/`a0')

}

else if `a0'<0 & `a1'<0{

`L'=(`a1'-`a0')/ln(`a1'/`a0')

}

L[i,j]=`L'

local j = `j'+1

}

local i = `i'+1

}

end

这是我自己摸索编的程序，还不知道怎么同时输入两个矩阵，syntax可以么。。。
你的好简洁啊，如果能够再加上这几个条件语句就好了

扫码加我拉你入群

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

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

bbs0805

2013-11-18 16:16:34

1、程序都编好了，就直接利用吧！不用再费事了。在你的程序中稍作修改即可：
cap program drop LMDI
program define LMDI
version 12
syntax namelist(min=2 max=2) //修改
gettoken Mat0 Mat1:namelist //添加.
.
.

end

2、输入矩阵的方法如下：
mat Mat0=(1,-2,3\2,0,-1\4,5,3)
mat Mat1=(2,0,-3\-2,6,1\0,1,3)
详见matrix的帮助

扫码加我拉你入群

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

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

XIEShichen

2013-11-18 16:35:10

bbs0805 发表于 2013-11-18 16:16
1、程序都编好了，就直接利用吧！不用再费事了。在你的程序中稍作修改即可：
cap program drop LMDI
prog ...

请问如何调用这个函数啊

我把这个函数存着personal/_Myado里面了
在profile里也引用了这个地址

输入LMDI Mat0 Mat1
显示unrecognized command: = invalid command name

扫码加我拉你入群

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

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

bbs0805

2013-11-18 18:31:07

我仔细看了你的程序，发现程序中错误太多，包括：宏的定义、宏的引用、关系表达式、矩阵运算、矩阵元素提取……
还是一个重要的是：计算式也有错误，计算结果会出现缺失值，例如，当a0=2,a1=-2时，(`a1'-`a0')/(ln(-`a1'/1E-200)+ln(1E-200/`a0'))的计算结果就是缺失值。
现对计算式你自己仔细修改，以你原来的计算式为例，程序语法修改后如下：
cap program drop LMDI
program define LMDI,rclass
version 12
syntax namelist(min=2 max=2)
gettoken Mat0 Mat1:namelist
local i = 1
local r = rowsof(`Mat0')
local c = colsof(`Mat0')
mat L=J(`r',`c',.)
while `i'<=`r' {
      local j = 1
  while `j'<=`c' {
      local a0= el(Mat0,`i',`j')
      local a1= el(Mat1,`i',`j')
            if `a0'>0 & `a1'>0{
   local L=(`a1'-`a0')/ln(`a1'/`a0')
            }
            else if `a0'>0 & `a1'==0{
   local L=(`a1'-`a0')/ln(1E-200/`a0')
            }
            else if `a0'>0 & `a1'<0{
   local L=(`a1'-`a0')/(ln(-`a1'/1E-200)+ln(1E-200/`a0'))
            }
            else if `a0'==0 & `a1'>0{
   local L=(`a1'-`a0')/ln(`a1'/1E-200)
            }
            else if `a0'==0 & `a1'==0{
   local L=0
            }
            else if `a0'==0 & `a1'<0{
   local L=(`a1'-`a0')/(ln(-`a1'/1E-200))
            }
            else if `a0'<0 & `a1'>0{
   local L=(`a1'-`a0')/(ln(`a1'/1E-200)+ln(-1E-200/`a0'))
            }
            else if `a0'<0 & `a1'==0{
   local L=(`a1'-`a0')/ln(-1E-200/`a0')
            }
            else if `a0'<0 & `a1'<0{
   local L=(`a1'-`a0')/ln(`a1'/`a0')
            }
            mat L[`i',`j']=`L'
            local j = `j'+1
      }
   local i = `i'+1
}
return matrix L=L
end
mat Mat0=(1,-2,3\4,0,-5\6,7,8)
mat Mat1=(9,10,-11\-12,13,14\15,16,17)
LMDI Mat0 Mat1
mat list r(L)

运行结果如下：

r(L)[3,3]
         c1       c2       c3
r1 3.6409569 7.4560192  -10.775174
r2  -14.563828 .02807278 18.45342
r3    9.82221 10.88694 11.939953

扫码加我拉你入群

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

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

点击查看更多内容…

XIEShichen

2013-11-18 19:06:17

bbs0805 发表于 2013-11-18 18:31
我仔细看了你的程序，发现程序中错误太多，包括：宏的定义、宏的引用、关系表达式、矩阵运算、矩阵元素提取 ...

不知道怎么感谢您了。多谢多谢多谢了~

stata就自学没两个月，都是是边学边用，见笑了
这段代码学习到好多，再次感谢：）

扫码加我拉你入群

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

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

XIEShichen

2013-11-18 22:51:44

bbs0805 发表于 2013-11-18 18:31
我仔细看了你的程序，发现程序中错误太多，包括：宏的定义、宏的引用、关系表达式、矩阵运算、矩阵元素提取 ...

我这基础上添加了两个矩阵，好像不能用了。

这个syntax namelist是不是有限制啊。
我看help里面也是给出了（min=2, max=2）的情况

cap program drop LMDI

program define LMDI,rclass

version 12

syntax namelist(min=4 max=4)

gettoken Y0 Y1 X0 X1:namelist

local r = rowsof(`Y0')

local c = colsof(`Y0')

mat L=J(`r',`c',.)

local i = 1

while `i'<=`r' {

local j = 1

while `j'<=`c' {

local y0= el(Y0,`i',`j')

local y1= el(Y1,`i',`j')

local x0= el(X0,`i',`j')

local x1= el(X1,`i',`j')

*------y0*y1==0

if `y0'==0 & `y1'==0{

local l=0

}

else if `y0'*`y1'==0 & `x0'*`x1'<>0 {

local l=0

}

else if `y0'*`y1'==0 & `x0'*`x1'==0 {

local l=`y1'-`y0'

}

*------y0*y1>0

else if `y0'*`y1'>0 & `x0'*`x1'>0 {

local l=(`y1'-`y0')/(ln(`y1'/`y0'))*ln(`x1'/`x0')

}

else if `y0'*`y1'>0 & `x0'*`x1'<0 & `x0'<0 {

local l=(`y1'-`y0')/(ln(`y1'/`y0'))*(ln(`x1'/1E-200)+ln(-1E-200/`x0'))

}

else if `y0'*`y1'>0 & `x0'*`x1'<0 & `x0'>0 {

local l=(`y1'-`y0')/(ln(`y1'/`y0'))*(ln(-`x1'/1E-200)+ln(1E-200/`x0'))

}

*------y0*y1<0 & x0*x1>0

else if `y0'*`y1'<0 & `x0'*`x1'>0 & `y0'<0 {

local l=((-1E-200-`y0')/ln(-1E-200/`y0')+(`y1'-1E-200)/ln(`y1'/1E-200))*ln(`x1'/`x0')

}

else if `y0'*`y1'<0 & `x0'*`x1'>0 & `y0'>0 {

local l=((1E-200-`y0')/ln(1E-200/`y0')+(`y1'-(-1E-200))/ln(-`y1'/1E-200))*ln(`x1'/`x0')

}

*------y0*y1<0 & x0*x1<0

*------ ->+ ->+

else if `y0'*`y1'<0 & `x0'*`x1'<0 & `y0'<0 & `x0'<0 {

local l=(-1E-200-`y0')/ln(-1E-200/`y0')*ln(-1E-200/`x0') + (`y1'-1E-200)/ln(`y1'/1E-200)*ln(`x1'/1E-200)

}

*------ ->+ +>-

else if `y0'*`y1'<0 & `x0'*`x1'<0 & `y0'<0 & `x0'>0 {

local l=(-1E-200-`y0')/ln(-1E-200/`y0')*ln(1E-200/`x0') + (`y1'-1E-200)/ln(`y1'/1E-200)*ln(-`x1'/1E-200)

}

*------ +>- ->+

else if `y0'*`y1'<0 & `x0'*`x1'<0 & `y0'>0 & `x0'<0 {

local l=(1E-200-`y0')/ln(-1E-200/`y0')*ln(-1E-200/`x0') + (`y1'-(-1E-200))/ln(`y1'/1E-200)*ln(`x1'/1E-200)

}

*------ +>- +>-

else if `y0'*`y1'<0 & `x0'*`x1'<0 & `y0'>0 & `x0'>0 {

local l=(1E-200-`y0')/ln(-1E-200/`y0')*ln(1E-200/`x0') + (`y1'-(-1E-200))/ln(`y1'/1E-200)*ln(-`x1'/1E-200)

}

*------

mat L[`i',`j']=`l'

local j = `j'+1

}

local i = `i'+1

}

return matrix L=L

end

扫码加我拉你入群

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

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

XIEShichen

2013-11-19 10:21:22

gettoken 有限制，最多只能拆分成两部分
args可以实现输入多个变量，不过引用的时候需要加`'
heihei

扫码加我拉你入群

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

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

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

分享

扫码加好友，拉您进群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群