【matplotlib】ガウス分布とローレンツ分布を合わせたフォークト関数(voigt)の作成方法と左右非対称化の方法[Python]

2024年9月11日

URLをコピーしました！

matplotlib

前回、Pythonのmatplotlibでピークを境に左右の形状が非対称な分布の作成方法を紹介しました。

フォークト関数の作成方法

フォークト関数はガウス分布とローレンツ分布がそれぞれある割合で足し合わされた分布になります。

つまりガウス分布がX％含まれているとしたら、ローレンツ分布は100-X％含まれていて、それらが足し合わされた分布というわけです。

ということでフォークト関数を自作関数とするとこんな感じになります。

中でガウス分布とローレンツ分布の関数も使っているので、忘れない様に追加してあります。

def gauss(x, a, m, s):
    curve = np.exp(-(x-m)**2/(2*s**2))
    return a * curve/max(curve)

def lorentz(x, a, m, g):
    curve = (1/np.pi)*(g/((x-m)**2 + g**2))
    return a * curve/max(curve)

def voigt(x, a, m, s, g, r):
    gauss_curve = gauss(x, a, m, s)
    lorentz_curve = lorentz(x, a, m, g)
    curve = r*gauss_curve + (1-r)*lorentz_curve
    return a * curve/max(curve)

表示まで含めたプログラム全体としてはこんな感じです。

import numpy as np
import matplotlib.pyplot as plt

x_list = np.arange(0, 10, 0.001)

a = 1
m = 5
s = 1
g = 1
r = 0.5

def gauss(x, a, m, s):
    curve = np.exp(-(x-m)**2/(2*s**2))
    return a * curve/max(curve)

def lorentz(x, a, m, g):
    curve = (1/np.pi)*(g/((x-m)**2 + g**2))
    return a * curve/max(curve)

def voigt(x, a, m, s, g, r):
    gauss_curve = gauss(x, a, m, s)
    lorentz_curve = lorentz(x, a, m, g)
    curve = r*gauss_curve + (1-r)*lorentz_curve
    return a * curve/max(curve)
    
y_gauss = gauss(x_list, a, m, s)
y_lorentz = lorentz(x_list, a, m, g)
y_voigt = voigt(x_list, a, m, s, g, r)

fig = plt.figure()
plt.clf()

plt.plot(x_list, y_gauss, label="Gauss")
plt.plot(x_list, y_lorentz, label="Lorentz")
plt.plot(x_list, y_voigt, label="Voigt")

plt.legend()
plt.show()

実行結果

フォークト関数を左右非対称化する方法

左右非対称な分布を作成する方法は前回紹介しましたので、詳しくはそちらの記事もご覧いただければと思います。

あわせて読みたい

簡単にいうと２つの異なる形状の分布を作成し、ピーク位置で分割し、左右異なる組み合わせで繋ぎ合わせるということをやります。

ということで関数部分だけ見てみるとこんな感じです。

def gauss(x, a, m, s):
    curve = np.exp(-(x-m)**2/(2*s**2))
    return a * curve/max(curve)

def lorentz(x, a, m, g):
    curve = (1/np.pi)*(g/((x-m)**2 + g**2))
    return a * curve/max(curve)

def voigt(x, a, m, s, g, r):
    gauss_curve = gauss(x, a, m, s)
    lorentz_curve = lorentz(x, a, m, g)
    curve = r*gauss_curve + (1-r)*lorentz_curve
    return a * curve/max(curve)

def asymmetric_voigt(x, a, m, s1, s2, g1, g2, r1, r2):
    peak_index = int(m / x_param[2])
    voigt_R = voigt(x, a, m, s1, g1, r1)
    voigt_L= voigt(x, a, m, s2, g2, r2)
    curve = np.concatenate([voigt_R[:peak_index], voigt_L[peak_index:]])
    return a * curve/max(curve)

表示まで含めたプログラム全体としてはこんな感じです。

import numpy as np
import matplotlib.pyplot as plt

x_param = [0, 10, 0.02]

x_list = np.arange(x_param[0], x_param[1], x_param[2])

a = 1
m = 5
s = 1
g = 1
r = 0.5

s1 = 1
s2 = 2
g1 = 1
g2 = 2
r1 = 0.3
r2 = 0.7
scale = 0

peak_index = m /x_param[2]

print(peak_index)


def gauss(x, a, m, s):
    curve = np.exp(-(x-m)**2/(2*s**2))
    return a * curve/max(curve)

def lorentz(x, a, m, g):
    curve = (1/np.pi)*(g/((x-m)**2 + g**2))
    return a * curve/max(curve)

def voigt(x, a, m, s, g, r):
    gauss_curve = gauss(x, a, m, s)
    lorentz_curve = lorentz(x, a, m, g)
    curve = r*gauss_curve + (1-r)*lorentz_curve
    return a * curve/max(curve)

def asymmetric_voigt(x, a, m, s1, s2, g1, g2, r1, r2):
    peak_index = int(m / x_param[2])
    voigt_R = voigt(x, a, m, s1, g1, r1)
    voigt_L= voigt(x, a, m, s2, g2, r2)
    curve = np.concatenate([voigt_R[:peak_index], voigt_L[peak_index:]])
    return a * curve/max(curve)
    
y_gauss = gauss(x_list, a, m, s)
y_lorentz = lorentz(x_list, a, m, g)
y_voigt = voigt(x_list, a, m, s, g, r)
y_asymmetric_voigt = asymmetric_voigt(x_list, a, m, s1, s2, g1, g2, r1, r2)

fig = plt.figure()
plt.clf()

plt.plot(x_list, y_gauss, label="Gauss")
plt.plot(x_list, y_lorentz, label="Lorentz")
plt.plot(x_list, y_voigt, label="Voigt")
plt.plot(x_list, y_asymmetric_voigt, label="Asymmetric Voigt")

plt.legend()
plt.show()

実行結果