通俗易懂的MonteCarlo积分方法(六)

这次我们主要以足球门的危险区域威胁度衡量问题(无防守球员)来作为MonteCarlo积分的一个实际应用

一.问题的提出

​ 一般的足球比赛中，球员在球门前不同位置起脚对球门的威胁时不同的。具体表现在：正前方的威胁要大于在球门两侧的射门，近距离的射门的威胁要远大于远距离射门。对于一般的职业球员来说，射门时候的速度一般在$10m/s$左右。我们要考虑的问题是球员在不同位置射门时候对球门的威胁度进行建模计算，并且绘制出相关的危险区域。

​ 而已知的标准球场长为$104m$，宽为$69m$,球门高$2.44m$，宽为$7.32m$。

二.问题分析

​ 球员在无论从那个地方射门都存在进或者不进两种状态，这本身就是一个随机事件。那么无非就是进球的概率最大的地方就是最危险的区域。

​ 事实上影响球员命中率的因素有很多，最重要的两点：球员基本素质(技术水平)和射门时候的位置。对于每个球员来说，其基本素质短时间几乎是无法改变的。(具体参见国足的情况)

​ 因此假设我们已经确定了每个球员的基本素质条件都是一样的，影响射门的因素只有位置的因素。

​ 某一个球员在球门前面向球门某目标射门时候，球员的素质和球员到目标点的距离决定了命中球门的概率。以球门中心为坐标原点$(0,0)$。我们假设在$(x,y)$点处命中的概率是$p(x,y)$，如果我们先不考虑球员素质的影响。一个球员如果射门那么在$x,y$两个方向上的同时射门的概率应该是互相独立的。也就是说:
$$p(x,y)=p(x)p(y)$$
​ 而对$p(x,y)$进行极坐标变换可知：
$$p(x,y)=p(rcos\theta ,rsin\theta )=g(r,\theta )$$
​ 而很明显射门的概率与角度无关，仅仅与距离有关，则:
$$p(x,y)=g(r,\theta )=g(r)$$
​ 整理一下可以得到:
$$p(x,y)=p(x)p(y)=g(r)=g(\sqrt{x^2+y^2})$$
​ 如果我们令$y=0$,那么$p(x)p(0)=g(x)$。带入上式可以得到:
$$p(x)p(y)=p(\sqrt{x^2+y^2})p(0)$$
​ 转化一下可以得到:
$$ln\frac{p(x)}{p(0)}+ln\frac{p(y)}{p(0)}=ln\frac{p(\sqrt{x^2+y^2})}{p(0)}$$
我们设$h(x)=ln\frac{p(x)}{p(0)}$。那么就有下面的式子满足:
$$h(x)+h(y)=h(\sqrt{x^2+y^2})$$
已知$h(x)=0$。那么肉眼就能看出来这个函数方程有一个特解:
$$h(x)=A{x}^2$$
那么我们甚至可以解出来：
$$p(x)=p(0)e^{A{x}^2}$$
而由归一化条件:
$${\int }_{-\infty }^{\infty }p(x)dx={\int }_{-\infty }^{\infty }p(0)e^{A{x}^2}dx=1$$
由高斯积分${\int }_{-\infty }^{\infty }{e}^{-x^2}dx=\sqrt{\pi }$可以解得:
$$p(x)=\sqrt{\frac{A}{\pi }}{e}^{-A{x}^2}$$
很明显，$p(x)\sim N(0,\frac{1}{\sqrt{2A}})$。从而$p(x,y)$就是服从的标准二维正态分布:
$$p(x,y)=p(x)p(y)=\frac{A}{\pi }{e}^{-A(x^2+y^2)}=\frac{A}{\pi }{e}^{-A{r}^2}$$
​ 很明显，球员从球场上射门得先确定一个目标点。射门之后依据概率落入到球门所在的平面。如果我们将球门视为一个区域，在该区域内对该分布积分，就可以得到这个射门能命中球门的概率。

​ 而球员在射门时，选在射门的目标点是任意的，而命中球门的概率对目标点的选择由很强的依赖性(主要是通过起始点到目标点的距离$r$来确定的)。那么我们遍历球门内部所有的点，将该区域内的不同命中的概率密度值乘以区域面积再求和，相当于做积分，将这个指标作为威胁度为确定球门的危险区域。这就是我们分析的一个大致的思路。

三.模型的建立

​ 我们设球门的底边中点位置为原点$O$，地面为$XOY$面，球门所在的平面为$YOZ$面。直线$AB$在地面上的投影与球门平面的夹角为$\theta$。$d$表示$A$到$B$点的直线距离，$k$表示衡量球员基本素质的指标。$D(x,y)$表示从$(x,y)$点射门的威胁度，$p(y,z)$表示从$A$点命中$B$点射门时候命中的概率。建立如下图所示的基本球门坐标系：

​ 当素质为$k$的球员从$A(x_0,y_0)$向距离为$d$的目标$B(y_1,z_1)$射门时球在球门平面所在的$YOZ$面$\Omega$上的落点呈现二维正态分布。其密度函数如下:
$$f(y,z)=\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{(y-y_1{)}^2+(z-z_1{)}^2}{2{\sigma }^2}},(y,z)\in \Omega$$
​ 其中${\sigma }^2$与球员素质$k$成反比，就是说，如果一个球员素质越高，那么他命中的稳定性应该更高，那么偏离标准点的平方和就越小，那么${\sigma }^2$就越小。同时，如果$AB$两点的距离$d$越长，那么命中的稳定性就会下降，就是说更容易偏离。可近似认为${\sigma }^2$与$d$成正比。当偏角$\theta =\frac{\pi }{2}$时,即垂直射门时，方差与角度无关，而角度越大，那么射门的稳定性越差，说明${\sigma }^2$越大。那么我们可以设比例系数为$1$。那么列出来的${\sigma }^2$的表达式为:
$${\sigma }^2=\frac{d}{k}(cot\theta +1)$$
​ 由几何关系可以推导出来：
$$\begin{gathered} cot\theta =\frac{|y_1-y_0|}{x_0} \\ d=\sqrt{x_0^2+(y_1-y_0{)}^2+z_1^2} \end{gathered}$$
​ 这里我么需要考虑一点：这个球只能落在地面之上，那么我么考虑球门区域为$D$: D { ( x , y , z ) ∣ x = 0 , 0 ≤ z ≤ 2.44 m , − 3.66 m ≤ y ≤ 3.66 m } D\{(x,y,z)|x=0,0\leq z \leq2.44m,-3.66m\leq y\leq 3.66m\} D{(x,y,z)∣x=0,0≤z≤2.44m,−3.66m≤y≤3.66m}。
$$\begin{gathered} P_D(x_0,y_0;y_1,z_1)={\iint }_{D}f(y,z)dydz={\iint }_D\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{(y-y_1{)}^2+(z-z_1{)}^2}{2{\sigma }^2}}dydz \\ {P}_{\Omega }(x_0,y_0;y_1,z_1)={\iint }_{\Omega }f(y,z)dydz={\iint }_{\Omega }\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{(y-y_1{)}^2+(z-z_1{)}^2}{2{\sigma }^2}}dydz \end{gathered}$$
​ 那么我们求条件概率$P_{D|\Omega }$表示这次射门命中的概率。
$$P_{D|\Omega }(x_0,y_0;y_1,z_1)=\frac{P_D(x_0,y_0;y_1,z_1)}{P_{\Omega }(x_0,y_0;y_1,z_1)}$$
​ 而对球场内的任意一点$(y_1,z_1)$做求和可求得所有的点对$A(x_0,y_0)$的威胁度。表示如下:
$$D(x_0,y_0)={\iint }_D{P}_{D|\Omega }(x_0,y_0;y_1,z_1)d{y}_{1}d{z}_1$$

四.模型的求解

​ 取参数$k=10$。在这里我们利用以下的关系：
$$P_{\Omega }(x_0,y_0;y_1,z_1)={\iint }_{\Omega }f(y,z)dydz={\iint }_{\Omega }\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{(y-y_1{)}^2+(z-z_1{)}^2}{2{\sigma }^2}}dydz\dot{=}1$$

利用$MentoCarlo$积分的方法求解，代码如下：

import math
import random
def distance(x_0,y_0,y_1,z_1):
    return math.sqrt(x_0**2+(y_1-y_0)**2+z_1**2)

def cotA(x_0,y_0,y_1,z_1):
    return math.fabs(y_1-y_0)/x_0

def sigma2(x_0,y_0,y_1,z_1):
    d = distance(x_0,y_0,y_1,z_1)
    k = 10
    cot = cotA(x_0,y_0,y_1,z_1)
    return d/k*(1+cot)

def f(y,z,x_0,y_0,y_1,z_1):
    return 1/(2*math.pi*sigma2(x_0,y_0,y_1,z_1))*math.pow(math.e,-((y-y_1)**2+(z-z_1)**2)/(2*sigma2(x_0,y_0,y_1,z_1)))

def P_D(numbers,x_0,y_0,y_1,z_1):
    a = 3.66
    b = 2.44
    S = 2*b*a
    sum = 0
    for i in range(numbers):
        y = -a+2*a*random.random()
        z = b*random.random()
        if (z<=b)&(z>=0)&(y>=-a)&(y<=a):
            sum += f(y,z,x_0,y_0,y_1,z_1)
    int2Value = S*1/numbers*sum
    return int2Value

def P_Omega(numbers,x_0,y_0,y_1,z_1):
    a = 1000
    S = 4*a**2
    sum = 0
    for i in  range(numbers):
        y = -a+2*a*random.random()
        z = -a+a*2*random.random()
        if (z<=a)&(z>=-a)&(y>=-a)&(y<=a):
            sum += f(y,z,x_0,y_0,y_1,z_1)
    int2Value = S*1/numbers*sum
    return int2Value

def P_D0(numbers,x_0,y_0,y_1,z_1):
    # return P_D(numbers,x_0,y_0,y_1,z_1)/P_Omega(numbers,x_0,y_0,y_1,z_1)
    return P_D(numbers,x_0,y_0,y_1,z_1)/1

def D(numbers,x_0,y_0):
    a = 3.66
    b = 2.44
    S = 2 * b * a
    sum = 0
    for i in range(numbers):
        y_1 = -a+2*a*random.random()
        z_1 = b*random.random()
        if (z_1<=b)&(z_1>=0)&(y_1>=-a)&(y_1<=a):
            sum += P_D0(numbers,x_0,y_0,y_1,z_1)
    int2Value = S*1/numbers*sum
    return int2Value

print(D(1000,3,1))

在$(3,1)$处射处球的威胁度可以计算为:

12.007925242429721

Process finished with exit code 0