一个简单的计量经济学程序学习分享

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

欢迎使用Markdown编辑器

经管之家：Do the best economic and management education！
Author: Daniel Tulpen Liu

2019年才认真的学习使用 Rstudio ,说实话很多年前也傻傻的使用过 R ，但是因为不是在 IDE 集成环境，特别难用，我记得以前我有一个学生叫周梦婷，国内某一个大学的经济学博士，好几年前就告诉我：老师，Rstudio 很强大，你会用吗？
我说我不会，我2019年下半年吧，开始使用 Rstudio ，才发现在 IDE 环境下使用，和傻傻的单独使用 R 是不一样的， 和大多数初学者一样。只使用 R Gui ，那么会发现数据管理，编程，都是非常困难的。 多年前我使用R 就是运行工具包里别人的代码演示，看看效果。

直到 2019年下半年，安装了 Rstudio ,也是基本不会，慢慢的摸索，下载资料学习，慢慢的才学会在 Rstudio 里会编程。 其实这样说了，论坛里很多人使用 Rstudio 的经验都比我的丰富。我算是后来学习的。

这里分享的是一个计算 TFP 的 R 工具包代码，但是从公开的工具包手册里是找不到这几个程序的，我把他贴出来。我呢，好像是从外网弄到的吧。

spfrontier: Spatial Stochastic Frontier Models

在出版物中使用程序包时引用‘spfrontier’:

Dmitry Pavlyuk (2019). spfrontier: Spatial Stochastic Frontier Models. R package version
0.2.5. https://CRAN.R-project.org/package=spfrontier

A BibTeX entry for LaTeX users is

@Manual{,
title = {spfrontier: Spatial Stochastic Frontier Models},
author = {Dmitry Pavlyuk},
year = {2019},
note = {R package version 0.2.5},
url = {https://CRAN.R-project.org/package=spfrontier},
}

注意：此引文信息是从软件包描述文件自动生成的，可能需要手工编辑，请参见‘help(“citation”)’。

routines

第一个程序

library(spdep)
library(pastecs)
library(car)
library(maptools)
library(maps)
library(ggmap)
library(corpcor)

rm(list=ls())
data(airports, envir = environment())


airports$Country<-as.factor(airports$Country)
airports$Routes <- airports$RoutesDeparture+airports$RoutesArrival
vars <- c("ICAO", "Country", "AirportName", "longitude", "latitude","Year", "PAX", "ATM", "Cargo", "Population100km", "Population200km", "Island", 
          "GDPpc","RunwayCount","CheckinCount","GateCount","ParkingSpaces", "RoutesDeparture", "RoutesArrival", "Routes")
airportsEU <- airports[vars]

airports2011 <- subset(airportsEU, Year==2011)
airports2011.complete <- with(airports2011, airports2011[complete.cases(PAX,ATM,Cargo,Population100km,Population200km,Island, GDPpc, RoutesDeparture,RoutesArrival),])
nrow(airports2011.complete)
stat.desc(airports2011.complete)

data <- airports2011.complete
nrow(data)
W <- constructW(cbind(data$longitude, data$latitude),data$ICAO)

formula <- log(PAX) ~ log(Population100km) + log(Routes) + log(GDPpc)


model000 <- spfrontier( formula , data=data)
summary(model000 )

model100 <- spfrontier(formula , data=data, W_y=W)
summary(model100 )

model010 <- spfrontier(formula , data=data, W_v=W, logging="debug")
summary(model010)

model001 <- spfrontier(formula , data=data, W_u=W, logging="debug")
summary(model001)

第二个程序, 第二个程序是这个工具包的作者公开发表的论文的代码；但是运行时间比较长。

复制粘贴这个程序，是可以运行的。数据就在这个工具包里。

rm(list=ls())

library(spdep)
library(pastecs)
library(car)
library(maptools)
library(maps)
library(ggmap)
library(corpcor)
library(spgwr)
library(raster)
library(RColorBrewer)
library(classInt)
library(frontier)
library(spfrontier)
library(gstat)



rowStdrt = function(W){
    for (j in 1:nrow(W)){
        W[j,] <- W[j,]/sum(W[j,])
    }
    return(W)
}


constructW <- function(coords, labels){
    W <- 1/as.matrix(dist(coords))
    colnames(W) <- labels
    rownames(W) <- labels
    W[which(W==Inf)] <- 0
    return(W)
}
#data <-read.csv("data/data.csv", header=T)
#save(data, file = "data/RTEpaper.rda")
data(RTEpaper, envir = environment())
data<-RTEpaper

#data <- subset(data, PAX > 4000000)





data$PAXperRunway <- data$PAX/(data$RunwayCount*1000000)
data$PAXperRunwayLength <- data$PAX/(data$RunwayLength*1000)
data$PAXperPop100 <- data$PAX/(data$Population100km)
data$TouristsPerkm<- data$Tourists/data$Area
#data <- subset(data, data$ICAO!="EGLC")
#data <- subset(data, data$ICAO!="ELLX")
stat.desc(data)

nrow(data)

data$OwnPublic <- ifelse(data$ownership=="public",1,0)
data$OwnMinorPrivate <- ifelse(data$ownership=="minor private",1,0)
data$OwnMajorPrivate <- ifelse(data$ownership=="major private",1,0)
data$OwnPrivate <- ifelse(data$ownership=="private",1,0)



W <- constructW(cbind(data$longitude, data$latitude),data$ICAO)

W <- rowStdrt(W)
listw <- mat2listw(W)

#formula <- log(PAX) ~ log(RunwayCount) +  log(Population100km)+  log(GDPpc)+  log(TouristsPerkm) + SouthIsland + hub + international
#formula <- log(PAX) ~ log(RunwayLength) +  log(Population100km)+  log(GDPpc)+  log(TouristsPerkm) + hub + international

#formula <- log(PAX) ~ log(Population100km) + log(Routes)+ log(RunwayCount) + log(CheckinCount)+ log(GateCount)+ log(GDPpc) + Island 

formula <- -log(PAX) ~ log(0.00001+cargo/PAX) + log(RunwayLength) +  log(Population100km)+  log(GDPpc)+  log(TouristsPerkm) + hub + international

eu.sfa <- spfrontier( formula , data=data, logging="info", costFrontier=T)
summary(eu.sfa)

ini <- c(coefficients(eu.sfa)$beta,0,coefficients(eu.sfa)$sigmaV, coefficients(eu.sfa)$sigmaU,0)

eu.ssf1000 <- spfrontier( formula , W_y=W,data=data, logging="info", costFrontier=T)
summary(eu.ssf1000)

eu.ols <- lm( formula , data=data)
summary(eu.ols)
vif(eu.ols)

resid <- stats::residuals(eu.ols)
dens <- density(resid)
dens <-  density(spfrontier::efficiencies(eu.sfa))
dens <-  density(spfrontier::efficiencies(eu.ssfa0010 ))


plot(dens,main="OLS residuals' kernel density")
x1 <- min(dens$x)  
x2 <- max(dens$x)

dd <- with(dens,data.fr ame(x,y))
library(ggplot2)
qplot(x,y,data=dd,geom="line",ylab="density",xlab="OLS residuals")+
  geom_ribbon(data=subset(dd,x>x1 & x<x2),aes(ymax=y),ymin=0,
              fill="red",colour=NA,alpha=0.2)

hist(spfrontier::efficiencies(eu.ssfa0010 ))
skewness(resid)



eu.sfa <- sfa( formula , data=data, ineffDecrease=FALSE)
summary(eu.sfa)


moran.test(resid,listw,randomisation=FALSE,alternative="two.sided")
moran.plot(resid,listw)
moran.plot(as.vector(spfrontier::efficiencies(eu.sfa)),listw)

moran.test(spfrontier::efficiencies(eu.sfa),listw,randomisation=FALSE,alternative="two.sided")

lm.LMtests(eu.ols, listw, test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))

summary(lagsarlm(formula ,data=data, listw))
eu.sem <- errorsarlm(formula ,data=data, listw)
summary(eu.sem)

resid <- residuals(eu.sem) - eu.sem$lambda * W %*% residuals(eu.sem)

eu.sfa <- spfrontier( formula , data=data, logging="info", costFrontier=T)
ini <- c(coefficients(eu.sfa)$beta,coefficients(eu.sfa)$sigmaV, coefficients(eu.sfa)$sigmaU)
eu.ssfa0010 <- spfrontier( formula , data=data,W_v=W, logging="debug", 
                           initialValues=c(ini,0.5)+0.01, costFrontier=T,
                           control=list(optim.control = list(tol=1e-2, iterlim=2000,repeatNM=1,printLevel=999)))
summary(eu.ssfa0010)
save(eu.ssfa0010, file="eu.ssfa0010-review2-2.rData")
load(file="eu.ssfa0010-review.rData")

ini <- c(coefficients(eu.sfa)$beta, coefficients(eu.sfa)$sigmaV, coefficients(eu.sfa)$sigmaU)
eu.ssfa0001 <- spfrontier( formula , data=data, W_u=W, logging="debug", 
                           initialValues=c(ini, 0),
                           control=list(optim.control = list(tol=1e-2, iterlim=2000,repeatNM=1)))
summary(eu.ssfa0001)
save(eu.ssfa0001, file="eu.ssfa0001.rData")


eu.ssfa1000 <- spfrontier( formula , data=data, W_y=W)
summary(eu.ssfa1000)

eu.bw <- gwr.sel(formula, data=data,coords=cbind(data$latitude, data$longitude))

eu.gwr <- gwr(formula, data=data,coords=cbind(data$latitude, data$longitude), bandwidth=eu.bw, hatmatrix=TRUE)
eu.gwr
eu.gwr$SDF$latitude <- data$latitude
eu.gwr$SDF$longitude <- data$longitude
eu.gwr$SDF$ICAO <- data$ICAO
eu.gwr$SDF$log_RunwayLength_t <- eu.gwr$SDF[["log(RunwayLength)"]]/eu.gwr$SDF[["log(RunwayLength)_se"]]
eu.gwr$SDF$log_Population100km_t <- eu.gwr$SDF[["log(Population100km)"]]/eu.gwr$SDF[["log(Population100km)_se"]]
eu.gwr$SDF$log_GDPpc_t <- eu.gwr$SDF[["log(GDPpc)"]]/eu.gwr$SDF[["log(GDPpc)_se"]]
eu.gwr$SDF$log_TouristsPerkm_t <- eu.gwr$SDF[["log(TouristsPerkm)"]]/eu.gwr$SDF[["log(TouristsPerkm)_se"]]

spplot(eu.gwr$SDF, "gwr.e", add=TRUE)
spplot(eu.gwr$SDF, "log.RunwayCount.")
spplot(eu.gwr$SDF, "log_RunwayCount_t")
spplot(eu.gwr$SDF, "log_Population100km_t")
spplot(eu.gwr$SDF, "log_GDPpc_t")
spplot(eu.gwr$SDF, "log.Population100km.")
spplot(eu.gwr$SDF, "log.GDPpc._se")


wd <- "D:/Dmitry/Dropbox/science/PhD/R/spfrontier/phd"
setwd(wd)
source("utils.R")

library(maptools)
library(classI)

world<-readShapeSpatial("TM_WORLD_BORDERS-0.3/TM_WORLD_BORDERS-0.3.shp", proj4string=CRS("+proj=longlat"))


world <- world[!is.na(world$FIPS),]

europe <- world[world$REGION == 150 | world$FIPS == "CY",]
eu <- europe

nclr <- 10
plotclr <- colorRampPalette(c("yellow", "red"))(nclr)
classI <- classIntervals(eu.gwr$SDF[["log_RunwayCount_t"]], nclr, style = "equal")
brks <- round(classI$brks, 1)
classI <- classIntervals(eu.gwr$SDF[["log_RunwayCount_t"]], nclr, style = "fixed",fixedBreaks = brks)
colcode <- findColours(classI, plotclr)

plot(eu, border="grey", xlim = c(0 , 10) , ylim = c(30 , 70))
points(eu.gwr$SDF$longitude,eu.gwr$SDF$latitude,pch=20,col=colcode)
legend("bottomleft", legend = names(attr(colcode, "table")), title="Title",fill = attr(colcode, "palette"), cex = 1, bty = "n")


nclr <- 7
dat <- residuals(eu.ols)
dat <- spfrontier::efficiencies(eu.sfa)
dat <- eu.gwr$SDF[["log(TouristsPerkm)"]]
dat <- eu.gwr$SDF[["log(Population100km)"]]
dat <- eu.gwr$SDF[["log(GDPpc)"]]
dat <- eu.gwr$SDF[["log(RunwayLength)"]]

dat <- eu.gwr$SDF$log_TouristsPerkm_t
dat <- eu.gwr$SDF$log_Population100km_t
dat <- eu.gwr$SDF$log_GDPpc_t
dat <- eu.gwr$SDF$log_RunwayLength_t
dat <- data$effSSFA-data$effSFA

par( mfrow = c(1, 1 ) )

plotclr <- colorRampPalette(c("yellow", "red"))(nclr ) 
classI <- classIntervals(dat, nclr, style = "equal")
brks <- round(classI$brks, 1)
classI <- classIntervals(dat, nclr, style = "fixed",fixedBreaks = brks)
colcode <- findColours(classI, plotclr)
plot(eu, border="darkgrey",col="gray90", xlim = c(-5 , 0) , ylim = c(29 , 70))
points(data$longitude,data$latitude,pch=21,col="darkgrey",bg=colcode,cex=1)
legend("bottomleft", legend = names(attr(colcode, "table")), title="Tourists: t-values",fill = attr(colcode, "palette"), cex = 1, bty = "n")



data$effSFA <- efficiencies(eu.sfa)
data$effSSFA <- spfrontier::efficiencies(eu.ssfa0010)
boxplot(as.vector(data$effSFA), as.vector(data$effSSFA),names=c("SFA","SSFA"),main="Efficiency levels")

d<-data[order(data$effSSFA-data$effSFA , decreasing=T),]
head(subset(d, select=c(ICAO,AirportName,Country,PAX,effSFA,effSSFA)),5)
tail(subset(d, select=c(Country,ICAO,AirportName,PAX,effSFA,effSSFA)),5)

d<-data[order(data$effSSFA , decreasing=T),]
subset(d, select=c(ICAO,effSFA,effSSFA))

d<-data[order(data$PAXperRunwayLength , decreasing=T),]
subset(d, select=c(ICAO,PAXperRunwayLength))


coordinates(data) =  ~longitude + latitude
proj4string(data) =  "+proj=longlat +datum=WGS84" 

x.range <- as.numeric(c(min(data$longitude), max(data$longitude))) 
y.range <- as.numeric(c(min(data$latitude), max(data$latitude))) 
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 0.1), y = seq(from = y.range[1], to = y.range[2], by = 0.1))  
coordinates(grd) <- ~ x + y
proj4string(grd) =  "+proj=longlat +datum=WGS84" 
gridded(grd) <- TRUE
plot(grd, cex = 1.5, col = "grey")
points(data, pch = 1, col = "red", cex = 1)
idw1 <- idw(formula = data$effSFA ~ 1, locations = data, newdata = grd,idp=2)
idw.output <- as.data.fr ame(idw1)
names(idw.output)[1:3]<- c("long", "lat", "var1.pred")
dat <- data$effSFA
nclr <- 7
plotclr <- colorRampPalette(c("cyan", "orange"))(nclr ) 
classI <- classIntervals(dat, nclr, style = "equal")
brks <- round(classI$brks, 1)
classI <- classIntervals(dat, nclr, style = "fixed",fixedBreaks = brks)
colcode <- findColours(classI, plotclr)
ggplot() + geom_tile(data = idw.output, alpha = 1, aes(x = long, y = lat, 
    fill = round(var1.pred, 1))) + scale_fill_gradient(low = "cyan", high = "orange") + 
    geom_path(data = eu, aes(long, lat, group = group), colour = "black") + xlim(x.range) +ylim(y.range)+
    geom_point(data = as.data.fr ame(data), aes(x = longitude, y = latitude), shape = 21, fill= colcode, col= colcode)+
	labs(fill = "Efficiency level", title = "Efficiency of EU airports (SF model)")+ theme_bw()




v1 <- variogram(data$effSFA~1, data, cutoff=500, width=10)
lo <- loess(v1$gamma~v1$dist, span=1.5)
plot(v1$dist,v1$gamma, xlab="Distance", ylab="Semi-variance",pch=20, col="blue")
lines(v1$dist,predict(lo), col='red', lwd=2)

ggplot() + geom_point(data = as.data.fr ame(v1), aes(x = dist, y = gamma), shape = 20, fill= "blue", col="blue")+
geom_line(data = as.data.fr ame(v1), aes(x = dist, y = predict(lo)), col="red", size=2) +
labs(x = "Distance",y="Semi-variance", title = "Variogram of SF efficiencies values")


geo.dists <- dist(cbind(data$longitude, data$latitude))
value.dists <- dist(data$effSFA)
mantel <- mantel.rtest(geo.dists, value.dists, nrepet = 999)

Enjoy

Yours :

Daniel Tulpen Liu.

Date： Today .

扫码加我 拉你入群

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

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

分享0 收藏3 回帖

关键词：计量经济学 计量经济 程序学习 经济学 environment

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R是自由软件，不带任何担保。
在某些条件下你可以将其自由散布。
用'license()'或'licence()'来看散布的详细条件。

R是个合作计划，有许多人为之做出了贡献.
用'contributors()'来看合作者的详细情况
用'citation()'会告诉你如何在出版物中正确地引用R或R程序包。

用'demo()'来看一些示范程序，用'help()'来阅读在线帮助文件，或
用'help.start()'通过HTML浏览器来看帮助文件。
用'q()'退出R.


R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R是自由软件，不带任何担保。
在某些条件下你可以将其自由散布。
用'license()'或'licence()'来看散布的详细条件。

R是个合作计划，有许多人为之做出了贡献.
用'contributors()'来看合作者的详细情况
用'citation()'会告诉你如何在出版物中正确地引用R或R程序包。

用'demo()'来看一些示范程序，用'help()'来阅读在线帮助文件，或
用'help.start()'通过HTML浏览器来看帮助文件。
用'q()'退出R.

> citation("spfrontier")

在出版物中使用程序包时引用‘spfrontier’:

  Dmitry Pavlyuk (2019). spfrontier: Spatial Stochastic Frontier Models. R package version
  0.2.5. https://CRAN.R-project.org/package=spfrontier

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {spfrontier: Spatial Stochastic Frontier Models},
    author = {Dmitry Pavlyuk},
    year = {2019},
    note = {R package version 0.2.5},
    url = {https://CRAN.R-project.org/package=spfrontier},
  }

注意：此引文信息是从软件包描述文件自动生成的，可能需要手工编辑，请参见‘help("citation")’。

$$
\tag{spec test}
\nabla (nq) \rightarrow \chi^2(n-k) \quad n \in  \quad   \Re
$$

给学生的教学，暂时发论坛数学公式。截图给学生。

$ \chi $

今天测试新的 LATEX 公式和使用方法。

$$
(1 -\phi_{1}\operatorname{B} )\ (1 -\Phi_{1}\operatorname{B}^{\operatorname{4}} )\ (1 - \operatorname{B}) (y_{t} -\delta\operatorname{t}) = (1 +\theta_{1}\operatorname{B} )\ (1 +\Theta_{1}\operatorname{B}^{\operatorname{4}} )\ \varepsilon_{t}
$$

在R 语言环境里做了新的测试， 再继续测试数学公式的生成。 哈哈

\begin{alignat}{2}
&\text{let}\quad &&y_{t} = \operatorname{y}_{\operatorname{0}} +\delta\operatorname{t} +\beta_{1}\operatorname{x1}_{\operatorname{t}} +\beta_{2}\operatorname{x2}_{\operatorname{t}} +\eta_{t} \\
&\text{where}\quad  &&(1 -\phi_{1}\operatorname{B} )\ (1 -\Phi_{1}\operatorname{B}^{\operatorname{4}} )\ (1 - \operatorname{B}) \eta_{t}  \\
& &&= (1 +\theta_{1}\operatorname{B} )\ (1 +\Theta_{1}\operatorname{B}^{\operatorname{4}} )\ \varepsilon_{t} \\
&\text{where}\quad &&\varepsilon_{t} \sim{WN(0, \sigma^{2})}
\end{alignat}

自动生成的，真先进。 嘿嘿，今天发财了。

$$
\operatorname{bill\_length\_mm} = \alpha + \beta_{1}(\operatorname{bill\_depth\_mm}) + \beta_{2}(\operatorname{flipper\_length\_mm}) + \epsilon
$$

这次多试一下几个；

$$
\operatorname{bill\_length\_mm} = \alpha + \beta_{1}(\operatorname{bill\_depth\_mm}) + \beta_{2}(\operatorname{island}_{\operatorname{Dream}}) + \beta_{3}(\operatorname{island}_{\operatorname{Torgersen}}) + \beta_{4}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Dream}}) + \beta_{5}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Torgersen}}) + \epsilon
$$


\begin{aligned}
bill\_length\_mm &= \alpha + \beta_{1}(bill\_depth\_mm) + \beta_{2}(island_{Dream}) + \beta_{3}(island_{Torgersen})\ + \\
&\quad \beta_{4}(bill\_depth\_mm \times island_{Dream}) + \beta_{5}(bill\_depth\_mm \times island_{Torgersen}) + \epsilon
\end{aligned}

$$
\begin{aligned}
\operatorname{bill\_length\_mm} &= \hat{\phi} + \hat{\gamma}_{1}(\operatorname{bill\_depth\_mm}) + \hat{\gamma}_{2}(\operatorname{island}_{\operatorname{Dream}}) + \hat{\gamma}_{3}(\operatorname{island}_{\operatorname{Torgersen}})\ + \\
&\quad \hat{\gamma}_{4}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Dream}}) + \hat{\gamma}_{5}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Torgersen}}) + \epsilon
\end{aligned}
$$

$$
\begin{aligned}
\log\left[ \frac { P( \operatorname{sex} = \operatorname{male} ) }{ 1 - P( \operatorname{sex} = \operatorname{male} ) } \right] &= \alpha + \beta_{1}(\operatorname{species}_{\operatorname{Chinstrap}}) + \beta_{2}(\operatorname{species}_{\operatorname{Gentoo}}) + \beta_{3}(\operatorname{bill\_length\_mm})\ + \\
&\quad \beta_{4}(\operatorname{species}_{\operatorname{Chinstrap}} \times \operatorname{bill\_length\_mm}) + \beta_{5}(\operatorname{species}_{\operatorname{Gentoo}} \times \operatorname{bill\_length\_mm})
\end{aligned}
$$

$$
\begin{aligned}
\operatorname{sex} &\sim Bernoulli\left(\operatorname{prob}_{\operatorname{sex} = \operatorname{male}}= \hat{P}\right) \\
\log\left[ \frac { \hat{P} }{ 1 - \hat{P} } \right]
&= \alpha + \beta_{1}(\operatorname{species}_{\operatorname{Chinstrap}}) + \beta_{2}(\operatorname{species}_{\operatorname{Gentoo}}) + \beta_{3}(\operatorname{bill\_length\_mm})\ + \\
&\quad \beta_{4}(\operatorname{species}_{\operatorname{Chinstrap}} \times \operatorname{bill\_length\_mm}) + \beta_{5}(\operatorname{species}_{\operatorname{Gentoo}} \times \operatorname{bill\_length\_mm})
\end{aligned}
$$

挺多的。 哈哈

新的测试：

$$
\widehat{Y}=\arg\max_{k=1,\cdots,K}{(\mathbf{X}-\frac{\boldsymbol{μ}_k}{2})^T\boldsymbol{β}_k+\logπ_k}.
$$

这个比较厉害

太极无极 发表于 2021-6-17 20:44
这个比较厉害

谢谢， 其实还有一些其他代码，也很厉害，还没分享论坛。

也是 R 的程序代码。