博客访问： 2305281
博文数量： 983
博客积分： 13195
博客等级： 上将
技术积分： 7431
注册时间：
C++,python,热爱算法和机器学习
IT168企业级官微
微信号：IT168qiye
系统架构师大会
微信号：SACC2013
分类： 信息化
文章来源：
PDF下载地址：
看了Stanford的Andrew Ng老师的机器学习公开课中关于Logistic Regression的讲解，然后又看了《机器学习实战》中的LogisticRegression部分，写下此篇学习笔记总结一下。
首先说一下我的感受，《机器学习实战》一书在介绍原理的同时将全部的算法用源代码实现，非常具有操作性，可以加深对算法的理解，但是美中不足的是在原理上介绍的比较粗略，很多细节没有具体介绍。所以，对于没有基础的朋友（包括我）某些地方可能看的一头雾水，需要查阅相关资料进行了解。所以说，该书还是比较适合有基础的朋友。
本文主要介绍以下三个方面的内容：
（1）Logistic Regression的基本原理，分布在第二章中；
（2）Logistic Regression的具体过程，包括：选取预测函数，求解Cost函数和J(θ)，梯度下降法求J(θ)的最小值，以及递归下降过程的向量化（vectorization），分布在第三章中；
（3）对《机器学习实战》中给出的实现代码进行了分析，对阅读该书LogisticRegression部分遇到的疑惑进行了解释。没有基础的朋友在阅读该书的Logistic Regression部分时可能会觉得一头雾水，书中给出的代码很简单，但是怎么也跟书中介绍的理论联系不起来。也会有很多的疑问，比如：一般都是用梯度下降法求损失函数的最小值，为何这里用梯度上升法呢？书中说用梯度上升发，为何代码实现时没见到求梯度的代码呢？这些问题在第三章和第四章中都会得到解答。
文中参考或引用内容的出处列在最后的“参考文献”中。文中所阐述的内容仅仅是我个人的理解，如有错误或疏漏，欢迎大家批评指正。下面进入正题。
2. 基本原理
Logistic Regression和Linear Regression的原理是相似的，按照我自己的理解，可以简单的描述为这样的过程：
（1）找一个合适的预测函数（Andrew Ng的公开课中称为hypothesis），一般表示为h函数，该函数就是我们需要找的分类函数，它用来预测输入数据的判断结果。这个过程时非常关键的，需要对数据有一定的了解或分析，知道或者猜测预测函数的“大概”形式，比如是线性函数还是非线性函数。
（2）构造一个Cost函数（损失函数），该函数表示预测的输出（h）与训练数据类别（y）之间的偏差，可以是二者之间的差（h-y）或者是其他的形式。综合考虑所有训练数据的“损失”，将Cost求和或者求平均，记为J(θ)函数，表示所有训练数据预测值与实际类别的偏差。
（3）显然，J(θ)函数的值越小表示预测函数越准确（即h函数越准确），所以这一步需要做的是找到J(θ)函数的最小值。找函数的最小值有不同的方法，Logistic Regression实现时有的是梯度下降法（Gradient Descent）。
3. 具体过程
3.1& 构造预测函数
Logistic Regression虽然名字里带“回归”，但是它实际上是一种分类方法，用于两分类问题（即输出只有两种）。根据第二章中的步骤，需要先找到一个预测函数（h），显然，该函数的输出必须是两个值（分别代表两个类别），所以利用了Logistic函数（或称为Sigmoid函数），函数形式为：
& & & & &&
对应的函数图像是一个取值在0和1之间的S型曲线（图1）。
接下来需要确定数据划分的边界类型，对于图2和图3中的两种数据分布，显然图2需要一个线性的边界，而图3需要一个非线性的边界。接下来我们只讨论线性边界的情况。
对于线性边界的情况，边界形式如下：
构造预测函数为：
hθ(x)函数的值有特殊的含义，它表示结果取1的概率，因此对于输入x分类结果为类别1和类别0的概率分别为：
3.2& 构造Cost函数
Andrew Ng在课程中直接给出了Cost函数及J(θ)函数如式（5）和（6），但是并没有给出具体的解释，只是说明了这个函数来衡量h函数预测的好坏是合理的。
实际上这里的Cost函数和J(θ)函数是基于推导得到的。下面详细说明推导的过程。（4）式综合起来可以写成：
取似然函数为：
对数似然函数为：
最大似然估计就是要求得使l(θ)取最大值时的θ，其实这里可以使用梯度上升法求解，求得的θ就是要求的最佳参数。但是，在Andrew Ng的课程中将J(θ)取为（6）式，即：
因为乘了一个负的系数-1/m，所以J(θ)取最小值时的θ为要求的最佳参数。
3.3& 梯度下降法求J(θ)的最小值
求J(θ)的最小值可以使用梯度下降法，根据梯度下降法可得θ的更新过程：
式中为α学习步长，下面来求偏导：
上式求解过程中用到如下的公式：
因此，（11）式的更新过程可以写成：
因为式中α本来为一常量，所以1/m一般将省略，所以最终的θ更新过程为：
另外，补充一下，3.2节中提到求得l(θ)取最大值时的θ也是一样的，用梯度上升法求（9）式的最大值，可得：
观察上式发现跟（14）是一样的，所以，采用梯度上升发和梯度下降法是完全一样的，这也是《机器学习实战》中采用梯度上升法的原因。
3.4& 梯度下降过程向量化
关于θ更新过程的vectorization，Andrew Ng的课程中只是一带而过，没有具体的讲解。
《机器学习实战》连Cost函数及求梯度等都没有说明，所以更不可能说明vectorization了。但是，其中给出的实现代码确是实现了vectorization的，图4所示代码的32行中weights（也就是θ）的更新只用了一行代码，直接通过矩阵或者向量计算更新，没有用for循环，说明确实实现了vectorization，具体代码下一章分析。
文献[3]中也提到了vectorization，但是也是比较粗略，很简单的给出vectorization的结果为：
且不论该更新公式正确与否，这里的Σ(...)是一个求和的过程，显然需要一个for语句循环m次，所以根本没有完全的实现vectorization，不像《机器学习实战》的代码中一条语句就可以完成θ的更新。
下面说明一下我理解《机器学习实战》中代码实现的vectorization过程。
约定训练数据的矩阵形式如下，x的每一行为一条训练样本，而每一列为不同的特称取值：
约定待求的参数θ的矩阵形式为：
& & & & &&
先求x.θ并记为A：
求hθ(x)-y并记为E：
g(A)的参数A为一列向量，所以实现g函数时要支持列向量作为参数，并返回列向量。由上式可知hθ(x)-y可以由g(A)-y一次计算求得。
再来看一下（15）式的θ更新过程，当j=0时：
同样的可以写出θj，
综合起来就是：
综上所述，vectorization后θ更新的步骤如下：
（1）求A=x.θ；
（2）求E=g(A)-y；
（3）求θ:=θ-α.x'.E,x'表示矩阵x的转置。
也可以综合起来写成：
前面已经提到过：1/m是可以省略的。
4. 代码分析
图4中是《机器学习实战》中给出的部分实现代码。
sigmoid函数就是前文中的g(z)函数，参数inX可以是向量，因为程序中使用了Python的numpy。
gradAscent函数是梯度上升的实现函数，参数dataMatin和classLabels为训练数据，23和24行对训练数据做了处理，转换成numpy的矩阵类型，同时将横向量的classlabels转换成列向量labelMat，此时的dataMatrix和labelMat就是（18）式中的x和y。alpha为学习步长，maxCycles为迭代次数。weights为n维（等于x的列数）列向量，就是（19）式中的θ。
29行的for循环将更新θ的过程迭代maxCycles次，每循环一次更新一次。对比3.4节最后总结的向量化的θ更新步骤，30行相当于求了A=x.θ和g(A)，31行相当于求了E=g(A)-y，32行相当于求θ:=θ-α.x'.E。所以这三行代码实际上与向量化的θ更新步骤是完全一致的。
总结一下，从上面代码分析可以看出，虽然只有十多行的代码，但是里面却隐含了太多的细节，如果没有相关基础确实是非常难以理解的。相信完整的阅读了本文，就应该没有问题了！^_^。
【参考文献】
[1]《机器学习实战》——【美】Peter Harington
[2] Stanford机器学习公开课（）
阅读(20373) | 评论(2) | 转发(0) |
相关热门文章
给主人留下些什么吧！~~
：好文。原书这一关键推导过程被一句话略过了。
是的，这里推导详细 |
好文。原书这一关键推导过程被一句话略过了。
请登录后评论。Logistic 分布_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
Logistic 分布
阅读已结束，下载文档到电脑
想免费下载本文？
定制HR最喜欢的简历
你可能喜欢第四章&概率分布
给出了详尽的统计表。R
还提供了相关工具来计算累计概率分布函数P(X
≤x)，概率密度函数和分位数函数(给定q，符合P(X ≤
x) & q 的最小x就是对应的分位数)，和基于概率分布的计算机模拟。
概率分布&&&&
R对应的名字&&&&
β分布 &&&&&&beta
shape1, &&&&shape2,
二项式分布& &binom size,
分布 &cauchy location, scale
卡方分布 &&&&chisq
&&&&&&&&&&df,
指数分布&&&
&exp&&&&&&&&&&&
F分布&&&&&&
&&&&&&&&&&&&&df1,
γ分布 &&&&&&gamma&&&&&&&&
几何分布 &&&&geom
&&&&&&&&&&prob
超几何分布 &&hyper
&&&&&&&&&m,
对数正态分布 lnorm &&&&&&&&&meanlog,
分布 logis &&&&&&&&location,
负二项式分布 nbinom &&&&&&&&size,
正态分布 &&&&norm&&&&&&&&&
分布 pois&&&&&&&&&&
分布 &&&&&&t
&&&&&&&&&&&&&df,
均匀分布 &&&&unif
&&&&&&&&&&min,
分布 weibull &&&&&&&shape,
分布 Wilcox&&&&&&&
不同的名字前缀表示不同的含义，d表示概率密度函数，p 表示累积分布函数（cumulative distribution function，CDF），q 表示分位函数以及r 表示随机模拟(random
deviates)或者随机数发生器。dxxx 的第一个参数是x，pxxx是q，qxxx是p，和rxxx的是n
(rhyper 和rwilcox例外，二者的参数是nn)。偏态指数（non-centrality parameter）ncp 现在仅用于累积分布函数，大多数概率密度函数和部分其他情况
：更细节的内容可以参考在线帮助文档。pxxx
函数都有逻辑参数lower.tail
和log.p, dxxx
也有一个逻辑参数log。这样就让计算函数的各种累积概率值成为可能。例如可以通过下式直接计算累计(“积分的”) 风险 （hazard）函数，H(t)
= &log(1 & F(t))，
- pxxx(t, ..., lower.tail = FALSE, log.p =
或用来计算更精确的对数似然值(dxxx(..., log =
TRUE))。此外还有函数ptukey 和qtukey 计算来自正态分布样本的标准化全距（studentizedrange）的分布。
这里是一些例子：
& ## t分布的双侧p值
& 2*pt(-2.43, df = 13)
& ## F(2, 7)分布的上1%分位数
& qf(0.99, 2, 7)
检验一个数据集的分布
我们可以用很多方法分析一个单变量数据集的分布。最简单的办法就是直接看数字。利用函数summary 和fivenum 会得到两个稍稍有点差异的汇总信息。 此外，stem (“茎叶”图)也会反映整个数据集的数字信息。
& attach(faithful)
summary(eruptions)
&&& Min. 1st Qu.
Median Mean 3rd Qu. Max.
&&& 1.600 2.163
4.000 3.488 4.454 5.100
fivenum(eruptions)
&& [1] 1.5 4.0000
&& & stem(eruptions)
&& The decimal point is 1
digit(s) to the left of the |
&& 22 | 3578
&& 24 | 00228
&&&26 | 23
&& 28 | 080
&& 32 | 2337
&& 34 | 250077
&& 38 | 5577
&& 48 | 00333
&& 50 | 0370
茎叶图和柱状图相似，R
用函数hist
绘制柱状图。
& hist(eruptions)
让箱距缩小，绘制密度图
& hist(eruptions, seq(1.6, 5.2, 0.2),
prob=TRUE)
& lines(density(eruptions, bw=0.1))
& rug(eruptions) #
显示实际的数据点
更为精致的密度图是用函数density
绘制的。在这个例子中，我们加了一条由density 产生的曲线。你可以用试错法（trial-and-error）选择带宽bw（bandwidth）因为默认的带宽值让密度曲线过于平滑（这样做常常会让你得到非常有“意思”的密度分布）。
我们可以用函数ecdf
绘制一个数据集的经验累积分布（empirical
cumulativedistribution）函数。
& plot(ecdf(eruptions), do.points=FALSE,
verticals=TRUE)
显然，这个分布和其他标准分布差异很大。那么右边的情况怎么样呢，就是火山爆发3分钟后的状况？我们可以拟合一个正态分布，并且重叠前面得到的经验累积密度分布。
& long &- eruptions[eruptions &
& plot(ecdf(long), do.points=FALSE,
verticals=TRUE)
& x &- seq(3, 5.4, 0.01)
& lines(x, pnorm(x, mean=mean(long),
sd=sqrt(var(long))), lty=3)
分位比较图（Quantile-quantile (Q-Q)
plot）便于我们更细致地研究二者的吻合程度。
par(pty="s") #
设置一个方形的图形区域
qqnorm(long); qqline(long)
上述命令得到的QQ图表明二者还是比较吻合的，但右侧尾部偏离期望的正态分布。我们可以用t 分布获得一些模拟数据以重复上面的过程
x &- rt(250, df = 5)
qqnorm(x); qqline(x)
这里得到的QQ图常常会出现偏离正态期望的长尾区域(如果是随机样本)。我们可以用下面的命令针对特定的分布绘制Q-Q图
qqplot(qt(ppoints(250), df = 5), x, xlab = "Q-Q
plot for t dsn")
最后，我们可能需要一个比较正规的正态性检验方法。R提供了Shapiro-Wilk 检验
& shapiro.test(long)
Shapiro-Wilk normality test
data: long
W = 0.9793, p-value = 0.01052
和Kolmogorov-Smirnov
& ks.test(long, "pnorm", mean = mean(long), sd
= sqrt(var(long)))
One-sample Kolmogorov-Smirnov test
data: long
D = 0.0661, p-value = 0.4284
alternative hypothesis: two.sided
(注意一般的统计分布理论（distribution
theory）在这里可能无效，因为我们是用同样的样本对正态分布的参数进行估计的。)
&单样本和双样本检验
到现在为止，我们已经学会了单样本的正态性检验。而更常见的操作是比较两个样本的特征。在R 里面，所有“传统”的检验都放在包stats 里面。这个包常常会自动载入。
下面是冰融化过程的潜热（latent
heat）(cal/gm)
数据（来自Rice (1995,
Method A: 79.98 80.04 80.02 80.04 80.03 80.03
80.04 79.97
80.05 80.03 80.02 80.00 80.02
Method B: 80.02 79.94 79.98 79.97 79.97 80.03
79.95 79.97
盒状图（boxplot）为这两组数据提供了简单的图形比较。
A &- scan()
79.98 80.04 80.02 80.04 80.03 80.03 80.04
80.05 80.03 80.02 80.00 80.02
B &- scan()
80.02 79.94 79.98 79.97 79.97 80.03 79.95
boxplot(A, B)
从图上可以直观的看出第一组数据相比第二组数据倾向给出较大的值。
为了比较两个样本的均值是否相等，我们可以使用非配对 t-检验
t.test(A, B)
&& Welch Two Sample
3.2499, df = 12.027, p-value = 0.00694
alternative hypothesis: true difference in means is not
equal to 0
percent confidence interval:
estimates:
of x mean of y
上面的结果表明在正态前提下 3
，二者有明显的统计差异。R
函数默认两个样本方差不齐，而SPLUS
类似函数t.test
则默认方差齐性。如果两个样本都是来自正态群体，我们可以用F检验来确定方差的齐性情况，
var.test(A, B)
test to compare two variances
0.5837, num df = 12, denom df = 7, p-value =
alternative hypothesis: true ratio of variances is not
equal to 1
percent confidence interval:
&0.1251097
estimates:&
of variances
&0.5837405
这表明二者方差在统计学上没有显著差异，我们可以采用传统的假设方差齐性的t-检验。
& t.test(A, B,
var.equal=TRUE)
Two Sample t-test
data: A and B
t = 3.4722, df = 19, p-value =
alternative hypothesis: true difference in
means is not equal to 0
95 percent confidence
sample estimates:
mean of x mean of y
所有这些检验都假设了数据的正态性。双样本的Wilcoxon
(或者Mann-Whitney)检验没有正态性的前提，仅仅要求在无效假设(null
hypothesis)情况下样本来自一个常规的连续分布。
wilcox.test(A, B)
Wilcoxon rank sum test with continuity
correction
89, p-value = 0.007497
alternative hypothesis: true mu is not equal to
Warning message:
compute exact p-value with ties in: wilcox.test(A,
注意警告信息：在两个样本中都有同秩现象,
这表明这些数据来自离散分布(可能由于数据的近似处理造成)。有好多种方法可以图形化的显示两个样本的差别。我们已经看过盒状图的比较。下面的命令
& plot(ecdf(A), do.points=FALSE,
verticals=TRUE, xlim=range(A, B))
& plot(ecdf(B), do.points=FALSE,
verticals=TRUE, add=TRUE)
同时显示两个样本的经验累计概率分布，而qqplot 得到的是两个样本的Q-Q图。Kolmogorov-Smirnov
检验是对两个经验累计概率分布间的最大垂直距离进行统计的。Kolmogorov-Smirnov 检验只假定数据服从一个常规的连续分布：
ks.test(A, B)
Two-sample Kolmogorov-Smirnov test
0.5962, p-value = 0.05919
alternative hypothesis: two.sided
Warning message:
compute correct p-values with ties in: ks.test(A,
已投稿到：
以上网友发言只代表其个人观点，不代表新浪网的观点或立场。Ⅱ型截尾广义Logistic分布的统计推断
广义Logistic分布(Generalized Logistic Distribution，GLD)最初由Balakrishnan和Leung在1988年首次提出，之后很多学者做了进一步研究。作为Logistic分布的推广，除了具有位置、尺度参数外，还具有能够影响偏度和峰度的形状参数，是一类偏态分布，目前在诸多领域得到广泛应用。该分布具有五种类型，本文主要研究截尾样本下Ⅰ型GLD的统计推断问题。　　本文研究了三参数GLD的参数估计问题。传统估计方法如极大似然估计法、矩法和概率加权矩法虽然已被广泛应用，但在实际使用时存在各种问题或局限，如:极大似然估计的算法可能不收敛。且大多是针对两参数GLD给...展开
广义Logistic分布(Generalized Logistic Distribution，GLD)最初由Balakrishnan和Leung在1988年首次提出，之后很多学者做了进一步研究。作为Logistic分布的推广，除了具有位置、尺度参数外，还具有能够影响偏度和峰度的形状参数，是一类偏态分布，目前在诸多领域得到广泛应用。该分布具有五种类型，本文主要研究截尾样本下Ⅰ型GLD的统计推断问题。　　本文研究了三参数GLD的参数估计问题。传统估计方法如极大似然估计法、矩法和概率加权矩法虽然已被广泛应用，但在实际使用时存在各种问题或局限，如:极大似然估计的算法可能不收敛。且大多是针对两参数GLD给出的估计方法。另外，已有的估计方法大多是在完全样本条件下给出的，而在应用领域常出现截尾样本，截尾样本的估计方法与完全样本不同。即使有关于截尾样本的参数估计方法，大多也都是对Ⅴ型GLD提出的。为解决这些问题，本文主要研究了完全样本下的矩估计、概率加权矩估计、L矩估计、LH矩估计等矩估计方法和极大似然估计，进而推广推导出了可在Ⅱ型截尾样本下做参数估计的偏概率加权矩估计、极大似然估计。并对针对各估计方法进行Monte Carlo模拟，对比研究估计效果。　　此外，本文还研究了Ⅰ型GLD的记录值。主要研究了Ⅰ型GLD低记录值的性质，以及如何通过低记录值进行参数估计。本文给出了位置参数μ和尺度参数σ的最佳线性无偏估计。另外，还给出了通过已观测数据对未观测到的低记录值的最佳线性无偏预测。之后进行数据模拟，研究估计和预测效果。收起
学位授予单位
机标分类号
本文读者也读过
加载中，请稍候
．客服电话
&&8:00-11:30,13:00-17:00(工作日)三参数Ⅰ型广义Logistic分布的参数估计--《北京工业大学》2015年硕士论文
三参数Ⅰ型广义Logistic分布的参数估计
【摘要】：本论文研究广义Logistic分布Generalized Logistic Distribution, GLD)的参数估计问题,主要分为三部分：第一部分是关于GLD的介绍及其参数估计研究现状介绍；第二部分利用多种估计方法对GLD进行参数估计,并给出部分估计的大样本性质；第三部分进行模拟研究,并对结果进行比较。在文中,主要利用L-矩估计、概率加权矩估计、广义概率加权矩估计、极大似然估计、最小二乘估计等对三参数Ⅰ型GLD进行参数估计,给出估计的一些大样本性质,进行数值模拟,并对结果进行比较与讨论。三参数Ⅰ型GLD的最小二乘估计和广义概率加权矩估计为首次提出,在模拟中发现最小二乘估计比较稳定,但没有某个参数范围内最小二乘估计绝对优于其他估计；广义概率加权矩估计对于形状参数的估计比较好,但u值的选取没有表现出一定的规律。这些均可作为继续研究的对象。在模拟研究之后,本文提出利用多种方法的特点和适用范围进行混合估计的思想。
【关键词】：
【学位授予单位】：北京工业大学【学位级别】：硕士【学位授予年份】：2015【分类号】：O212.1【目录】：
摘要4-5Abstract5-7第1章 绪论7-13 1.1 研究意义7-9 1.2 国内外研究现状9-10 1.3 准备工作10-13第2章 三参数Ⅰ型广义Logistic分布的参数估计13-24 2.1 关于广义Logistic分布13-15 2.2 概率加权矩估计15-16 2.3 L-矩估计16-17 2.4 广义概率加权矩估计17-19 2.5 极大似然估计19-21 2.6 最小二乘估计21-23 2.7 本章小结23-24第3章 模拟计算及讨论24-29 3.1 除广义概率加权矩估计外其他估计的模拟计算24-27 3.2 广义概率加权矩估计的模拟计算27-28 3.3 本章小结28-29结论29-30参考文献30-34致谢34
欢迎：、、)
支持CAJ、PDF文件格式
【相似文献】
中国期刊全文数据库
张二艳;许燕;;[J];北京印刷学院学报;2009年06期
吴慎言;;[J];西安工程科技学院学报;1992年01期
吴刘仓,黄文亮,詹金龙;[J];昆明理工大学学报(自然科学版);2001年02期
崔博文,陈剑,陈心昭,任章;[J];安徽大学学报(自然科学版);2005年03期
肖庭延;张培培;阎金华;;[J];高等学校计算数学学报;2005年S1期
车金锐;梁刚;陈敏茹;;[J];数学的实践与认识;2006年04期
刘湘蓉;;[J];应用概率统计;2006年04期
李兵;朱宁;唐文芳;段复建;;[J];统计与决策;2007年23期
吴如光;;[J];西安文理学院学报(自然科学版);2008年02期
蔺焕泉;;[J];长春大学学报;2010年02期
中国重要会议论文全文数据库
陈欢;杨德森;时胜国;;[A];2008年全国声学学术会议论文集[C];2008年
范龙;白鸽;;[A];第二届中国卫星导航学术年会电子文集[C];2011年
张崇高;;[A];模糊集理论与模糊应用专辑——中国系统工程学会模糊数学与模糊系统委员会第十届年会论文选集[C];2000年
崔博文;陈剑;陈心昭;任章;;[A];第二十三届中国控制会议论文集（上册）[C];2004年
刘娜;张端金;杨林;;[A];2011年中国智能自动化学术会议论文集（第一分册）[C];2011年
中国博士学位论文全文数据库
赵媛媛;[D];中国科学技术大学;2011年
中国硕士学位论文全文数据库
苏千叶;[D];山东大学;2015年
徐欢;[D];上海大学;2015年
刘翔;[D];北京工业大学;2015年
卢冬晖;[D];湖北师范学院;2014年
桂咏新;[D];华中师范大学;2003年
于家富;[D];山东师范大学;2006年
刘洪伟;[D];东北林业大学;2012年
王剑;[D];华中科技大学;2007年
尤太林;[D];贵州财经学院;2011年
詹小平;[D];东北师范大学;2006年
&快捷付款方式
&订购知网充值卡
400-819-9993
《中国学术期刊（光盘版）》电子杂志社有限公司
同方知网数字出版技术股份有限公司
地址：北京清华大学 84-48信箱 大众知识服务
出版物经营许可证 新出发京批字第直0595号
订购热线：400-819-82499
服务热线：010--
在线咨询：
传真：010-
京公网安备75号