随机过程（6）——泊松过程三大变换，更新过程引入

学弱猹

发布于 2021-08-10 11:32:23

1.7K0

发布于 2021-08-10 11:32:23

文章被收录于专栏：学弱猹的精品小屋学弱猹的精品小屋

上一节笔记：随机过程（5）——无限状态马尔科夫链的进一步探讨，泊松分布引入，复合泊松分布

————————————————————————————————————

大家好！

上一节，我们给大家介绍了泊松过程和复合泊松过程（微信公众号的版本，“过程”写成了“分布”造成了一些误会，读者注意下～）。这一节我们会介绍针对泊松过程的一些变换，掌握它们就可以洞察到更多有趣的泊松过程的性质。具体是什么意思呢？看到下面就知道了。当然，如果有空的话，还会简单介绍一点更新过程的概念。

另外，马尔可夫和马尔科夫在笔记中会混着用（取决于我编辑文档所用的电脑），因此大家看到这两个词，也不要奇怪，因为它们都是可接受的翻译。

那么我们开始吧。

泊松过程变换

稀疏

稀疏（Thinning）变换简单来说，就是考虑每一次到达中，满足特定条件的子集。把这一部分子集抽出来，就会得到一条新的链条。我们后面会提到说，这链条也是一个泊松过程。

把这一个说法用数学描述一下，就是下面这个定理

Theorem 1: 设

N(t)

是一个速率为

\lambda

的泊松过程，

Y_i

是一个取正整数的随机变量，并且假设每一次到达都会得到一个“奖赏”

Y_i

。设

N_j(t)

表示符合

Y_i = j

的

的个数（

i \le N(t)

），证明

N_j(t)

是一个速率为

\lambda P(Y_i = j)

的泊松过程，并且与原始的泊松过程相互独立。

这个过程大概可以用图表示为下面这样。

这里为了比较简单的说明问题，我们假设这里的

Y_i

只会取两个值，也就是1或2。我们后面会发现，这样并不影响一般的证明结果。也就是说，我们实际上把原来的

N(t)

拆分为了

N_1(t)

和

N_2(t)

。

注意这里要说明两点，一是两个随机变量

N_1(t), N_2(t)

满足的分布（虽然我们一般用“过程”来描述它们，但它们本质上就是两个随机变量），二是它们俩之间的独立性。所以合在一起，我们就可以通过最后计算出来的分布的可分离性，一下子说清楚这两点。

注意到

P(N_1(t+s) - N_1(s) = k, N_2(t + s) - N_2(s) = l)

=P(N(t + s) - N(s) = k + l, k ~ \text{elements} = 1, l ~ \text{elements} =2 ~ \text{in} \{ Y_1, \cdots, Y_{k + l}\})

=P(N(t + s) - N(s) = k + l) P(k ~ \text{elements} = 1, l ~ \text{elements} =2 ~ \text{in} \{ Y_1, \cdots, Y_{k + l}\})

= \frac{(\lambda t) ^{k + l}}{(k + l)!}e^{-\lambda t} \binom{k + l}{k} p^k(1-p)^l = \frac{(\lambda p t)^k}{k!}e^{-\lambda p t}\frac{(\lambda (1-p)t)^l}{l!}e^{-\lambda (1-p)t}

到最后其实就说明了我们的结论成立，因为前一半对应

Poisson(\lambda pt)

，后一半对应

Poisson(\lambda (1-p)t)

。

这一串公式推导还是挺迷惑人的，我们一行一行解释哈。

第一行到第二行，其实是相当于把两条分散的泊松过程又“拼在了一起”。这可以用图表示为如下。

可以等价过来的原因就是，这两句话说的就是同一件事。虽然现在来看还是挺显然的，但说句实话确实不容易往这个方向去想到。

第二行到第三行可以看出就是用了独立性，这是因为这两件事本质上是无关的。一条泊松过程中间到达了

k + l

次，和这

k + l

次中间到底到达几次1和几次2，当然相互之间是不会产生影响的。

第三行到最后一行其实就是一个概率公式的替代。这里就可以解释为什么我们可以只讨论取两个值的情况，因为如果我们取了多个值，也只是充其量最后的式子

\binom{k + l}{k} p^k(1-p)^l

变成多项分布的公式而已。这个细节我们留给读者补充，并不困难。

可能有的人会疑惑在于，只证明两条泊松过程在相同时间段内的独立性，就可以把两个泊松过程的独立性证明出来吗？答案是肯定的。下面这张图可以很好的解释。

这个性质其实最有趣的地方也就在于，所有的“子过程”相互之间是互不干扰的。简单来说，改变一个泊松子过程的性质，不会影响到其他泊松子过程。这个观念还挺反直觉的，因为都是从一个泊松过程中“稀疏”出来的，结果相互之间却毫无关联。

通过这个稀疏的性质，其实还可以推出下面这个性质。

Corollary 1: Nonhomogeneous Thinning 考虑一个速率为

\lambda

的泊松过程，假设点到达数轴上的点

的概率密度函数为

p(s)

，那么这一些点服从一个速率为

\lambda p(s)

的泊松过程。

说它是theorem 1的推广，很大程度就是因为在这里，我们假设了每一次到达不同点的概率都不一样。这和之前的泊松过程的差别就在于，之前我们只讨论“到达”这一个概念，而什么时候到达仅仅取决于速率

\lambda

和指数分布

\exp(\lambda)

。但是这个情况就相当于多了一个制约因素

p(s)

，因此这里得到的泊松过程就是一个更复杂一些的结果。

用数学来说，其实就是希望证明

N(t +s ) - N(s) \sim Poisson(\int_s^{t + s} \lambda p(r) dr)

但这个的证明讲白了就是定积分求解的梯形公式那一套。我们可以把区间

[s, t + s]

拆分为无数个小块，每一块有一个小的增量。假如说某一块是

[s_i, s_i + \Delta s_i]

，那么这一块因为

\Delta s_i

很小，所以可以认为这一块基本上到达的概率都是

p(s_i)

。那么实际上这一块我们可以得到

N(s_i + \Delta s_i) - N(s_i) \sim Poisson(\lambda p(s_i) \Delta s_i)

这个如何利用上Theorem 1呢？简单来说，对于这一个区间，因为到达的概率会变化，所以我们不妨把它抽象为两个状态，一个是到达，得到奖赏为1；另外一个也是到达，但是奖赏为0。这样的话其实相当于在讨论一个同质（homogenous）的随机过程，也就是我们之前一直讨论的泊松过程。而拆出的“奖赏为1”的那一部分，就是我们这个问题所研究的泊松过程。所以为什么我们可以得到系数

\lambda p(s_i)

，其实是来源于这么一个思考。

有了这个之后，我们把各个区间加在一起，利用泊松分布的可加性不难得到

N(s + t) - N(s) \sim Poisson(\sum_i \lambda p(s_i) \Delta s_i)

接下来，我们令

\Delta s_i \to 0

，就可以得到我们的结论。

异质泊松过程（Non-homogeneous）的一个最大的特点就是到达每一个点的概率不一样，所以理解起来要比同质随机过程，要更加困难一些。

我们来举两个例子吧。

Problem 1:

M / G / \infty

考虑一个电话排队问题。有无穷多条电话线，来电话的个数是一个速率为

\lambda

的泊松过程，每一个电话所需要的时间服从一个分布

G(\cdot)

，但只知道

G(0) = 0

且期望为

\mu

。问长期来看（时间

t \to \infty

），这个系统里有多少个电话？

M/G/\infty

是对于这一个排队系统的简单描述。

表示电话到来是服从马尔可夫性的，

表示分布，而

\infty

表示电话线的个数，简单来说就是这个系统中，因为电话线无穷多，所以不会有等待的问题。

这个问题怎么建模呢？首先我们来看，在一个有限时间

内，电话到来的情况可以怎么样描述。大致可以画出下面这一张图。

那么问题就是，在

[0, t]

这个时间段内，会有多少个人还在打电话。我们不妨假设有一个人在点

处到达，那么如果它在这个时间段内有打电话，它的电话时长一定要大于

t- s

。而这个概率为

1 - G(t - s)

。总而言之，我们把问题抽象为下面的说法

有

1 - G(t -s )

的概率这个人打不完电话，所以会被计数，相当于到达，得到奖赏为1。其他情况，电话会打完，就不会被计数，相当于到达，但得到奖赏为0。

所以我们可以得到，

N_1(t) \sim Poisson(\lambda \int_0 ^ t (1 - G(t -s )) ds)

，这是用到了Corollary 1的结论。

有了这个之后，我们只需要令

t \to \infty

就可以了，我们可以得到

\int_0^t (1 - G(t -s )) ds = \int_0^t (1 - G(r))dr \to \int_0^\infty (1 - G(r))dr = \mu

最后一步是用到了期望的一个经典公式

E(X) = \int_0^\infty P(X \ge x)dx

。

有了这个之后，就可以知道，长期来看，会有

\lambda \mu

这么多人正在打电话。这是因为泊松分布，参数就是它的期望。

Problem 2: 考虑一个道路上汽车行驶的问题。已知每分钟，在某地经过的汽车数目服从速率为

\frac 23

的泊松过程。其中10%是卡车，90%是轿车。如果知道刚才的一个小时内，经过了10辆卡车，问刚才的一个小时内，经过的车的数量的期望是多少？

这个题乍一看，10辆卡车，占了10%，所以张口就来会经过100辆车。如果得到了这么一个答案，就说明第一印象往往会给人带来很严重的偏差。

首先这是一个很明显的可以通过拆分（也就是稀疏）解决的问题。既然我们知道汽车到来的速率为

\frac 23

，那么根据稀疏，我们可以把汽车拆分为“卡车”和“轿车”，那么根据Theorem 1就不难得到，卡车到来的速率为

\frac 2 {30}

，轿车到来的速率为

\frac {18}{30}

，并且二者相互独立。

正经计算一下，我们可以得到

E(N(60)|N_1(60) = 10)=E(N_1(60)|N_1(60) = 10) + E(N_2(60)|N_1(60) = 10)

= 10 + E(N_2(60)) = 10 + 60 \times \frac {18}{30} = 46

所以如果得到100这个答案，就说明忽视了两条“子泊松过程”的独立性。当然如果不清楚这个独立性，也很难得到46这个答案。

叠加

叠加（superposition）就是稀疏的逆向操作，理解起来并不困难。

Theorem 2: 设

N_1(t), \ldots, N_k(t)

是相互独立的泊松过程，且速率分别为

\lambda_1, \ldots, \lambda_k

。那么

\sum_{i = 1}^k N_i

是一个速率为

\sum_{i = 1}^k \lambda _i

的泊松过程。

证明其实和我们之前说的差不多。我们考虑拆分出不同的区间

t_0 < t_1 < \cdots < t_n

。那么很明显，

N(t_{i + 1}) - N(t_i)

相互独立。另一方面，我们根据条件，又有

N(t_1) - N(t_0) = \sum_{i = 1}^k (N_i(t_1) - N_i(t_0)) \sim Poisson(\sum_{i = 1}^k \lambda_i (t_1 - t_0))

这自然是根据独立性和泊松分布可加性得到的。

接下来就是把各个区间加在一起。很容易观察到不相交区间所带来的独立性，因此再次利用泊松分布可加性就可以得到我们这里的结论。

这个证明看上去很容易，不过我们还是希望介绍另一种证明，因为这可以帮助我们更好的理解泊松过程中，到达相距时间

\tau_i

的相关分析。

可以看下面这张图。

以两条泊松过程为例，在把这两条泊松过程加在一起之后，从某一个到达点

到下一个到达点中间距离的时间，其实可以写为

\tau_1^1 \land \tau_1^2

。而这个事实上相当于，在知道

\tau_1^1 \sim \exp(\lambda_1)

，

\tau_1^2 \sim \exp(\lambda _2)

的时候，问

\tau_1^1 \land \tau_1^2 = \min(\tau_1^2, \tau_1^2)

所服从的分布。概率论告诉我们，这个分布就是

\exp(\lambda_1+ \lambda_2)

，所以这个叠加性也就不难说明了。

好的，我们再来看两个例子。

Problem 3: 已知有两条泊松过程，一条表示的是红色点的到达次数，到达速率为

\lambda

。另外一条表示的是绿色点的到达次数，到达速率为

\mu

。问有4个绿色点到达前，已经有6个红色点到达的概率。

这个题看似很复杂，但实际上可以充分利用泊松过程的变换，把这个题做一个极大的简化。

首先我们考虑做一个叠加，然后再把这个叠加出来的新的泊松过程再做稀疏，这样就相当于得到了一个拆分，这个拆分满足

每一次到达，有

\frac {\lambda }{\lambda + \mu}

的概率是红色的点，有

\frac \mu{\lambda + \mu}

的概率是绿色的点。

那么如果要求出现4个绿色点之前，要先出现6个红色点，事实上有很多情况。但不管怎么说，其实只需要前9个点，至少有6个点是红色点就可以了。但是在这么拆分之后，这个问题就完全变成了一个二项分布的问题，因为每一个点是什么颜色的概率都已经知道了，且每一个颜色出现的概率，又具有相互独立的性质。

所以剩下的任务，就是把它拆分成很多情况，然后利用二项分布公式求概率就可以了。所以容易得到答案是

\sum_{k = 6}^9 \binom{9}{k} p^ k(1-p)^{9 - k}

，这里的

p = \frac{\lambda}{\lambda + \mu}

。

取条件

取条件相对前两个来说不是特别好理解。大致可以解释为，在规定了泊松过程的一些条件之后，泊松过程内部的到达时间服从的分布近似于均匀分布。具体来说，设

T_1, T_2, \ldots

为泊松过程的到达时间，设

U_1,\ldots, U_n

服从均匀分布

U[0, t]

，且两两独立。

V_1, \ldots, V_n

是将

U_1, \ldots, U_n

作升序排序之后的结果。

那么取条件，其实就是这么一个定理

Theorem 3: 在

N(t) = n

的时候，

(T_1, \ldots, T_n)

与

(V_1, \cdots, V_n)

同分布。

这个定理其实要说明的就是，在

0 < t_1 < \cdots < t_n < t

的时候，

P(T_1 = t_1, \ldots, T_n = t_n |N(t) = n) = P(V_1 = t_1, \ldots, V_n = t_n)

但是

P(V_1 = t_1, \ldots, V_n = t_n) = \frac {n!}{t^n}

（想想为什么？），所以我们只需要推出之前那个概率

P(T_1 = t_1, \ldots, T_n = t_n |N(t) = n)

就可以了。注意到我们有

P(T_1 = t_1, \ldots, T_n = t_n |N(t) = n) = \frac{P(T_1 = t_1, \ldots, T_n = t_n ,N(t) = n)}{P(N(t) = n)}

=\frac{P(\tau_1 = t_1, \tau_2 = t_2 -t_1, \cdots, \tau_{n + 1} = t - t_n )}{P(N(t) = n)}

注意到

\tau_i

相互之间独立，代入指数分布和泊松分布的公式，就可以得到结论。这个计算我们交给读者完成。

事实上，在取条件上，还有另外一个定理。

Corollary 2: 若

s < t, 0 \le m \le n

，那么有

P(N(s) = m | N(t) = n) = \binom{n}{m} (\frac s t)^m (1 - \frac s t)^{n - m}

简单来说，这个定理就是在告诉我们，在

N(t) = n

的时候，有

N(s) \sim Binomial(n, \frac s t)

，也就是服从二项分布的意思。

这个定理的证明就是利用了Theorem 3。在

N(t)= n

的时候，我们研究

N(s)

，其实就是希望观察在

之前到达的个数，而因为“到达时间”这个量，在取了条件的时候服从的是一个均匀分布。所以相当于把概率变成了

P(m ~ \text{elements} \le s ~ \text {in} ~ \{U_1, \ldots, U_n\})

这其实也是一个二项分布问题（

\le s

的概率为

\frac t s

）。所以这个概率的计算很容易，而这就是我们的结论。

好的，我们来看一个取条件的例子吧。

Problem 4: 设

N(t)

服从一个速率为2的泊松过程。计算

P(N(3) = 4|N(1) = 1)

P(N(1) = 1|N(3) = 4)

这两个题非常好的对比出了泊松过程的两种取条件的情况，对应的做法也有很大不同。

对于第一个题，我们注意到，相当于知道了之前时间的情况，问之后的泊松过程的期望。因此可以考虑做一个变换，得到

P(N(3) = 4|N(1) = 1) = P(N(3) - N(1) = 3) = \frac{4^3}{3!}e^{-4}

这是因为两个泊松过程的差依然服从一个泊松分布。

第二个题就是反过来，知道了之后时间的情况，问之前的泊松过程的期望。这就是我们这一部分所介绍的内容，根据Corollary 2我们直接可以得到

P(N(1) = 1|N(3) = 4) = \binom{4}{1}\frac 13 (1- \frac 13)^3 = \frac {32}{81}

好的，到此为止，关于泊松过程的变换我们就介绍完了。

更新过程

更新过程（Renewal Process）是泊松过程的一个推广。在泊松过程中，我们假设了相邻两个到达之间相距的时间服从指数分布

\exp(\lambda)

，那么去掉这个假设，得到的就是一个更新过程。

具体用数学来说，就是下面这一段话。

Definition 1: Renewal Process 设

t_1, t_2, \ldots

为相邻两个点到达的时间，相互之间独立同分布，服从分布

F(\cdot)

，设

T_n = t_1 + \cdots + t_n

，

N(t) = \max \{n : T_n \le t\}

，那么

N(t)

就是一个更新过程。

这里的

t_i

就和泊松过程里的

\tau_i

是一个意思。另外容易得到的是

T_n \sim F^{*n}(t)

P(N(t) = n) = F^{*n}(t) - F^{*(n + 1)}(t)

这里的

F^{*n}(t)

表示

F(t)

的

次卷积。学过概率论的一定明白我在说什么，不过即使你对卷积毫无概念，也不影响后面的理解。

一个很有趣的例子就是之前经常讨论的马尔科夫链，假如说我们研究相邻两个状态

之间所满足的分布。那么不管从哪个

出发，到达下一个

中间经过的时间（当然在离散马尔科夫链中就是步数）所服从的分布一定都是相同的，否则就没有马尔科夫性了。下面这一张图表达了这个意思。

这一张图我们也在之前介绍马尔科夫链的时候给大家用过。

讲到这里其实最容易产生的疑问就是，一个一般的分布

F(\cdot)

所引出的更新过程，会有什么非常特殊或者可用的性质吗？研究起来，会不会难度太大？事实上因为这个推广的存在，对于更新过程的研究总是需要依赖一些更加一般，更加抽象的工具。有的时候相对泊松过程来说，也不是特别好理解。因此可能确实需要大家多一些耐心……

首先我们想先介绍一下在马尔科夫链中已经提到过的大数定律（Law of Large Numbers）。这个大数定律和概率论中的略有区别。

Theorem 4: Law of Large Numbers 设

\mu = E(t_1)

，若

P(t_1 > 0) > 0

，那么有

\frac{N(t)}{t} \to \frac 1 \mu ~ a.s. ~ t \to \infty

。

这里的

a.s.

是almost surely的意思，之前也已经提过两次了。

我们在之前证明大数定律（随机过程（3）——无限状态的平稳测度，返回时间，访问频率：几个定理的证明）的时候，就已经说了，一般这种极限会通过夹逼定理来完成。而一个很显然的结果是

T_{N(t)} \le t < T_{N(t) + 1}

所以简单来说，我们希望研究的就是

\frac{N(t)}{T_{N(t)}}

的极限情况，或者更一般的说，研究

\frac{n}{T_n}

的极限情况。

注意到

\frac{T_n}{n} = \frac{\sum_{i = 1}^n t_i}{n} \to \mu

这里用到的就是概率论中，同分布情况下的强大数定律。所以我们有

\frac {n}{T_n} \to \frac 1 \mu

。这个说明清楚之后，注意到

\frac{T_{N(t)}}{N(t)} \le \frac t {N(t)} \le \frac{T_{N(t) + 1}}{N(t) + 1} \cdot \frac{N(t) + 1}{N(t)}

所以令

t \to \infty

就可以了。

但是这个证明是不是万无一失的？有没有发现有一个条件

P(t_1 > 0) > 0

没用上？这个条件其实是用来保证

\frac{N(t) + 1}{N(t)} \to 1

的。这个大家可能已经习以为常的结论，还需要额外保证？当然，因为你并不知道

N(t) \to \infty

是不是对的。

现在我们假如说

N(t) \not \to \infty

，那么这样的话，相当于说，存在一个

满足

N(t) < M, \forall t

。这样的话就有

t_1 + \cdots + t_M > t, t \to \infty

这个必须要保证

t_{i} \to 0

，但是很明显这个是做不到的，因为这个情况至少要要求

P(t_1>0) = 0

，这和条件是矛盾的。

这个定理揭示了更新过程最基本的极限性质。

更新函数与更新方程

提出更新函数和更新方程的原因是，对于更新过程而言，其核心的性质自然就依赖于

F(\cdot)

这个函数。而更新函数和更新方程就是为了更好的去描述这么一个概念。

Definition 2: Renewal Function 定义

m(t) = E[N(t)]

为更新函数。 Proposition 1: Renewal Equation 更新方程为

m(t) = F(t) + \int_0^t m(t-s)f(s)ds

更新方程为什么可以写成那样，自然需要一些证明。这里的一个想法是，既然我们在研究

E[N(t)]

，那么不妨框定一个范围，只讨论一个

T_1

的情况。这是因为每一次到达之后，下一次就和上一次无关（每个

t_i

相互独立）了。

注意到

m(t) = E[N(t)] = E[N(t)]I(T_1 \le t) = E[I(T_1 \le t) E(N(t) |T_1)]

=\int_0^t E[N(t)|T_1 = s]f(s)ds

倒数第二个式子其实就是用到了重期望公式，我们希望考虑

T_1

的情况，所以自然要对

T_1

取个条件。

有了这个之后，注意到，可以得到

E[N(t)|T_1 = s] = E[N(t - s) + 1]

我们用到的，还是“马尔科夫性”（注意这里严格来说并不是马尔科夫性，我们用这个词，只是想简单的描述“前面与后面无关”这个含义）。虽然已经反反复复，但我们还是决定提醒大家仔细理解这一点，因为真的太重要了。

把这个式子代回，就可以得到最终的更新方程。这个代回的过程我们留给读者。

好的，来看一个计算的例子吧。

Problem 5: 若

F \sim U(0, 1)

，证明

m(t) = e^t - 1, 0 \le t \le 1

。

首先根据均匀分布的公式，我们有

F(t) = \begin{cases}0 & t \le 0 \\ t & 0 < t < 1 \\ 1 & t \ge 1\end{cases}

所以实际上的更新方程为

m(t) = t + \int_0^t m(t-s) ds, 0 < t < 1

简单化简一下我们有

m(t) = t + \int_0^t m(u)du

那么求导就可以得到

m'(t) = 1 + m(t)

，结合

m(0) = 0

（这是根据

N(0) = 0

），根据基本常微分方程的求解，就可以得到

m(t) = e^t - 1

。

事实上读者可以自己发现，对于泊松过程，我们有

m(t) = \lambda t

，而这是符合我们之前的认知的。

最后我们介绍一个后面也会用到，类似于大数定律的一个定理。这个定理我们不证明，原因在于它所依赖的工具比较高阶，对于《随机过程》这个系列来说没有太多了解的必要了。

Theorem 5: Elementary Renewal Theorem

\frac {m(t)} t \to \frac 1 \mu, t \to \infty

说它“类似于”大数定律，是因为确实长得很像。但是注意不能使用控制收敛定理的证明方式，这是因为

\lim_{t \to \infty} E(\frac{N(t)}t) \not \to E(\lim_{n \to \infty}\frac{N(t)}{ t})

如果使用控制收敛定理的话。原因也很简单，因为找不到一个函数

G(t)

，使得

|\frac{N(t)} t| \le G(t) ,\forall t

。

好的，那么关于更新过程的初步介绍，我们先说到这里。

小结

本节主要介绍了泊松过程的三大变换和更新过程的一些基本概念与性质。泊松过程的三大变换带来的一些应用很多时候也是有些反直觉的，其通过拆分得到的过程相互之间的独立性，以及取条件中得到的二项分布，都是泊松过程的亮点。更新过程作为泊松过程的拓展，在研究上具备更大的挑战，但我们到后面会发现，这样的一个拓展不仅仅是有必要的，也可以给我们带来研究问题的很多新的视角。

在下一节我们会开始研究更新奖赏过程，如果有空的话，会继续介绍排队模型，并引出排队论的一些定理与应用。

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2020-12-28，如有侵权请联系 cloudcommunity@tencent.com 删除

网站

本文分享自学弱猹的精品小屋微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

网站

登录后参与评论

0 条评论

热度