1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# %load "E:桌面space临时数据python个人自定义模块ImportFile.py"
# Standard Scientific Import
import numpy as np
import scipy as sp
import pandas as pd
import matplotlib as mpl
from matplotlib import pyplot as plt
from matplotlib.pyplot import plot as plot
import sklearn
import seaborn as sns
import sys
import patsy

# 个人代码测试路径
sys.path.append(r"C:UsersAdministratorPycharmProjectsQY_TS_Quant")
from QY_plot import *

plt.rcParams['font.sans-serif'] = ['SimHei']  # 中文字体设置-黑体
plt.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
sns.set(font='SimHei',font_scale=1.25,style="ticks",rc={"xtick.major.size": 3, "ytick.major.size": 3})# 解决Seaborn中文显示问题

1
2
3
4
5
# Jupyter 默认设置
%matplotlib inline
%config InlineBackend.figure_format="retina"
%config InlineBackend.rc = {"figure.figsize": (7.5,4.5)}

1
2
3
4
# 多列输出
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

1
2
3
from functools import reduce
import operator

1
2
3
4
# 假设在alpha,beta的射手，随机等概率可以往各个方向设计
def StrafingSpeed(alpha, beta, x):
    return (np.square(beta)+np.square(alpha-x))/beta

1
2
3
def normlize(x):
    return x/np.sum(x)

1
2
3
4
def MakePmf(alpha, beta, w=30):
    pmf = list(map(lambda v:1./ StrafingSpeed(alpha, beta, v), np.arange(w+1)))
    return normlize(pmf)

1
2
3
4
5
def likehood(AX, AY, x, w=30):
    AZS = reduce(operator.add, map(lambda v:1./StrafingSpeed(AX, AY, v),np.arange(w+1)))
    AZ =  1./StrafingSpeed(AX, AY, x)
    return  AZ/AZS

1
2
3
4
# 围场宽30， 长50
w = 30
h = 50

1
2
3
# 宽边 中弹位置
x = [15, 16, 18, 21]

1
2
3
4
alpha = 10
beta = 10
plt.plot(MakePmf(alpha, beta))

1
2
3
4
5
6
7
8
9
10
def update(AP,x):
    if not np.isscalar(x):
        for xi in x:
            AP = update(AP, xi)
    else:
        AZ = likehood(AX, AY, x)
        AP = AP*AZ
        AP /= AP.sum()
    return AP

1
2
3
4
X = np.arange(w+1)
Y = np.arange(1, h+1)
AX, AY = np.meshgrid(X, Y)

1
2
3
AP = np.ones_like(AX)
AP = update(AP, x)

1
2
3
plt.contour(AX,AY,AP,80)
plt.colorbar()

1
2
3
4
5
6
7
fig, axes = plt.subplots(1, 2, figsize=(18, 6))
axes[0].plot(AX[0], AP.sum(axis=0), label="alpha")
axes[1].plot(AY[:,0], AP.sum(axis=1), label="beta")
axes[0].legend();
axes[1].legend();
plt.suptitle("概率密度");

1
2
3
4
5
6
7
fig, axes = plt.subplots(1, 2, figsize=(18, 6))
axes[0].plot(AX[0], AP.sum(axis=0).cumsum(), label="alpha")
axes[1].plot(AY[:,0],AP.sum(axis=1).cumsum(), label="beta")
axes[0].legend();
axes[1].legend();
plt.suptitle("分布函数");

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def confintv(cdf, p):
    t = []
    flag = 0
    for i,v in enumerate(cdf):
        if flag==0:
            if v>=0.5-p/2:
                t.append(i)
                flag=1
        elif flag==1:
            if v>=0.5+p/2:
                t.append(i)
                break
    return t

1
2
3
# alpha的 50%置信区间
AX[0][confintv(AP.sum(axis=0).cumsum(), 0.5)]

1
2
array([14, 21])

1
2
3
# beta的50%置信区间
AY[:,1][confintv(AP.sum(axis=1).cumsum(), 0.5)]

1
2
array([ 5, 31])

1
2
3
4
5
Acolor = np.zeros_like(AP)
s = np.argsort(AP.ravel())
for v in np.searchsorted(AP.flat[s].cumsum(),[0.25, 0.5, 0.75]):
    Acolor.flat[s[:v]] += 0.25

1
2
plt.imshow(Acolor, origin=True, aspect=0.5)

1
2
3
4
clf = plt.contour(AX, AY, Acolor,3);
plt.clabel(clf, fontsize=10);
plt.contourf(AX, AY, Acolor, 3);