两个分布的组合拟合

Question

我有两个由函数

P_multiple

和

P_single

给出的独立分布。我想将两者结合成一个组合拟合函数，该函数在两者之间具有自然的“过渡”。你有什么建议吗？谢谢！

这是我得到的图与我想要得到的图（红色）。

import numpy as np
import matplotlib.pyplot as plt

def P_multiple(alpha):
    return 1 / np.sqrt(np.pi) * np.exp(-alpha**2)

def P_single(alpha, Z):
    alpha_abs = np.abs(alpha)
    return 1 / (alpha_abs**3 * 8 * np.log(204 * Z**(-1/3)))

Z_silicon = 14 
alpha_threshold = 3 
alpha_values = np.linspace(-7, 7, 300)

multiple = P_multiple(alpha_values)
single = P_single(alpha_values, Z_silicon)

multiple[np.abs(alpha_values) > alpha_threshold] = 0
single[np.abs(alpha_values) <= alpha_threshold] = 0

combined_distribution = multiple + single

plt.figure(figsize=(10, 6))
plt.plot(alpha_values, multiple, label='Multiple')
plt.plot(alpha_values, single, label='Single')
plt.plot(alpha_values, combined_distribution, label='Combined', linestyle='--')

plt.legend()
plt.yscale('log') 
plt.show()

Answer 1

有很多方法可以做到这一点，但我想我会找到一种与你在你想要的内容中展示的方法类似的方法，而我只实现了一侧，我也已经在这方面花了一点时间已经很久了，让它发挥作用应该不会太难 我决定使用 bezier 曲线 来实现这一点，它的缺点是使用 t 变量，这意味着很难找到某个 x 或某个 y 的精确值，除非你有很高的采样率，因为我真的没有知道你要实现这个的目的是什么，我只是这样做，因为这是我的第一个想法，让它看起来最接近 贝塞尔曲线本质上只是使用三点制作曲线的一种方法

import numpy as np
import matplotlib.pyplot as plt

def P_multiple(alpha):
    return 1 / np.sqrt(np.pi) * np.exp(-alpha**2)

def P_single(alpha, Z):
    alpha_abs = np.abs(alpha)
    return 1 / (alpha_abs**3 * 8 * np.log(204 * Z**(-1/3)))

def findcrosspoint(error, dotpoints, valuem, values):
    #finds the closest point of crossing by approximation
    for i in range (dotpoints-1, 0, -1):
        if abs(values[i] - valuem[i]) <= error:
            return i
    return -1

def bezierline(startpoint, endpoint, t):
    # creates a bezier line
    return startpoint + t * (endpoint-startpoint)

def beziercurve(startpoint, crosspoint, endpoint, res):
    curvevals = []
    for i in range(res):
        if i != 0:
            t = i / res
        else:
            t = i
            
        # curves between two bezier lines
        curvevals.append(bezierline(startpoint, crosspoint, t) + t * (bezierline(crosspoint, endpoint, t) - bezierline(startpoint, crosspoint, t)))
    return np.array(curvevals)

def transition(Z, dotpoints, viewsize, startpoint, thresh, transitionrange, xweight, yweight):
    # calculates all the values
    xs = np.linspace(startpoint, startpoint+viewsize, dotpoints)
    valuem = P_multiple(xs)
    values = P_single(xs, Z)
    transitionpoint = thresh

    # calculates start, cross and endpoints to calculate the bezier curve
    startpoint = transitionpoint - transitionrange
    startpoint = np.array([startpoint, np.log(P_multiple(startpoint))])
    crosspoint = findcrosspoint(0.1, dotpoints, np.log(valuem), np.log(values)) * (viewsize/dotpoints)
    
    # the x and y weight allow you to change the position of the cross point which can change the look of the curve
    crosspoint = np.array([crosspoint / xweight, np.log(P_multiple(crosspoint) / yweight)])
    plt.scatter(np.array([crosspoint[0]]), np.array([crosspoint[1]]), s=20, label='cross')
    plt.scatter(np.array([startpoint[0]]), np.array([startpoint[1]]), s=20, label='start')
    endpoint = transitionpoint + transitionrange
    endpoint = np.array([endpoint, np.log(P_single(endpoint, Z))])
    plt.scatter(np.array([endpoint[0]]), np.array([endpoint[1]]), s=20, label='end')

    # calculates the bezier curve values
    curvevals = beziercurve(startpoint, crosspoint, endpoint, 10000)
    # puts them all back into one array
    values = []
    inserted = False
    for i in range(xs.shape[0]):
        if xs[i] <= startpoint[0]:
            values.append([xs[i], np.log(P_multiple(xs[i]))])
        elif xs[i] >= endpoint[0]:
            values.append([xs[i], np.log(P_single(xs[i], Z))])
        elif inserted == False:
            for value in curvevals:
                values.append(np.array([value[0], value[1]]))
            inserted = True
    return np.array(values)


Z_silicon = 14
alpha_threshold = 3
alpha_values = np.linspace(-7, 7, 500)

multiple = P_multiple(alpha_values)
single = P_single(alpha_values, Z_silicon)

multiple[np.abs(alpha_values) > alpha_threshold] = 0
single[np.abs(alpha_values) <= alpha_threshold] = 0

combined_distribution = multiple + single

plt.figure(figsize=(10, 6))
plt.plot(alpha_values, np.log(multiple), label='Multiple')
plt.plot(alpha_values, np.log(single), label='Single')
values = transition(Z_silicon, 1000, 7, 0, 3, 1, 1, 1)
xs = np.zeros(values.shape[0])
ys = np.zeros(values.shape[0])
for i in range(values.shape[0]):
    xs[i] = values[i][0]
    ys[i] = values[i][1]
plt.plot(xs, ys, label='combined')

plt.legend()
plt.show()

虽然这段代码很粗糙，但它确实与您在图表中绘制的内容非常匹配。它是相当可调节的，所以你可以让它看起来更符合你的喜好

编辑： 忘了提及，当我将它与对数刻度一起使用时，这看起来相当奇怪，因此我记录了这些值，如果这带来不便，我很抱歉