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

足球门的危险区域威胁度衡量问题

一.问题的提出
二.问题分析
三.模型的建立
四.模型的求解

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

一.问题的提出

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

而已知的标准球场长为104m104m104m，宽为69m69m69m,球门高2.44m2.44m2.44m，宽为7.32m7.32m7.32m。

二.问题分析

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

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

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

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

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

三.模型的建立

我们设球门的底边中点位置为原点OOO，地面为XOYXOYXOY面，球门所在的平面为YOZYOZYOZ面。直线ABABAB在地面上的投影与球门平面的夹角为θ\thetaθ。ddd表示AAA到BBB点的直线距离，kkk表示衡量球员基本素质的指标。D(x,y)D(x,y)D(x,y)表示从(x,y)(x,y)(x,y)点射门的威胁度，p(y,z)p(y,z)p(y,z)表示从AAA点命中BBB点射门时候命中的概率。建立如下图所示的基本球门坐标系：

当素质为kkk的球员从A(x0,y0)A(x_0,y_0)A(x0,y0)向距离为ddd的目标B(y1,z1)B(y_1,z_1)B(y1,z1)射门时球在球门平面所在的YOZYOZYOZ面Ω\OmegaΩ上的落点呈现二维正态分布。其密度函数如下:
f(y,z)=12πσ2e−(y−y1)2+(z−z1)22σ2,(y,z)∈Ω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 f(y,z)=2πσ21e−2σ2(y−y1)2+(z−z1)2,(y,z)∈Ω
其中σ2\sigma^2σ2与球员素质kkk成反比，就是说，如果一个球员素质越高，那么他命中的稳定性应该更高，那么偏离标准点的平方和就越小，那么σ2\sigma^2σ2就越小。同时，如果ABABAB两点的距离ddd越长，那么命中的稳定性就会下降，就是说更容易偏离。可近似认为σ2\sigma^2σ2与ddd成正比。当偏角θ=π2\theta = \frac{\pi}{2}θ=2π时,即垂直射门时，方差与角度无关，而角度越大，那么射门的稳定性越差，说明σ2\sigma^2σ2越大。那么我们可以设比例系数为111。那么列出来的σ2\sigma^2σ2的表达式为:
σ2=dk(cotθ+1)\sigma^2 = \frac{d}{k}(cot \theta+1) σ2=kd(cotθ+1)
由几何关系可以推导出来：
cotθ=∣y1−y0∣x0d=x02+(y1−y0)2+z12cot \theta = \frac{|y_1-y_0|}{x_0}\\d =\sqrt{x_0^2+(y_1-y_0)^2+z_1^2} cotθ=x0∣y1−y0∣d=x02+(y1−y0)2+z12
这里我么需要考虑一点：这个球只能落在地面之上，那么我么考虑球门区域为DDD:D{(x,y,z)∣x=0,0≤z≤2.44m,−3.66m≤y≤3.66m}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}。
PD(x0,y0;y1,z1)=∬Df(y,z)dydz=∬D12πσ2e−(y−y1)2+(z−z1)22σ2dydzPΩ(x0,y0;y1,z1)=∬Ωf(y,z)dydz=∬Ω12πσ2e−(y−y1)2+(z−z1)22σ2dydzP_D(x_0,y_0;y_1,z_1) = \iint_Df(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 PD(x0,y0;y1,z1)=∬Df(y,z)dydz=∬D2πσ21e−2σ2(y−y1)2+(z−z1)2dydzPΩ(x0,y0;y1,z1)=∬Ωf(y,z)dydz=∬Ω2πσ21e−2σ2(y−y1)2+(z−z1)2dydz
那么我们求条件概率PD∣ΩP_{D|\Omega}PD∣Ω表示这次射门命中的概率。
PD∣Ω(x0,y0;y1,z1)=PD(x0,y0;y1,z1)PΩ(x0,y0;y1,z1)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)} PD∣Ω(x0,y0;y1,z1)=PΩ(x0,y0;y1,z1)PD(x0,y0;y1,z1)
而对球场内的任意一点(y1,z1)(y_1,z_1)(y1,z1)做求和可求得所有的点对A(x0,y0)A(x_0,y_0)A(x0,y0)的威胁度。表示如下:
D(x0,y0)=∬DPD∣Ω(x0,y0;y1,z1)dy1dz1D(x_0,y_0) = \iint_DP_{D|\Omega}(x_0,y_0;y_1,z_1)dy_1dz_1 D(x0,y0)=∬DPD∣Ω(x0,y0;y1,z1)dy1dz1

四.模型的求解

取参数k=10k=10k=10。在这里我们利用以下的关系：
PΩ(x0,y0;y1,z1)=∬Ωf(y,z)dydz=∬Ω12πσ2e−(y−y1)2+(z−z1)22σ2dydz=˙1P_{\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 PΩ(x0,y0;y1,z1)=∬Ωf(y,z)dydz=∬Ω2πσ21e−2σ2(y−y1)2+(z−z1)2dydz=˙1

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

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_0def sigma2(x_0,y_0,y_1,z_1):d = distance(x_0,y_0,y_1,z_1)k = 10cot = 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.66b = 2.44S = 2*b*asum = 0for 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*sumreturn int2Valuedef P_Omega(numbers,x_0,y_0,y_1,z_1):a = 1000S = 4*a**2sum = 0for 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*sumreturn int2Valuedef 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)/1def D(numbers,x_0,y_0):a = 3.66b = 2.44S = 2 * b * asum = 0for 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*sumreturn int2Valueprint(D(1000,3,1))

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

12.007925242429721Process finished with exit code 0