使用 scipy 拟合 Rice 分布

Question

我正在尝试为我拥有的一些大米分布式数据编写一个拟合程序，但由于某些（可能是愚蠢的）原因，它不起作用。

分布创建得很好，并且拟合路由似乎按照我习惯的高斯分布进行工作。然而，当我拟合曲线时，我只是得到了废话。似乎看不出我哪里出错了。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import rice
from scipy.optimize import curve_fit

# Custom Rice PDF function
def rice_pdf(x, nu, amplitude, b, scale):
    return (x / b) * np.exp(-(x**2 + scale**2) / (2 * b**2)) * amplitude

# Function to fit a Rice distribution to a histogram using curve_fit
def fit_rice_distribution_to_histogram(hist_data, bins):
    # Calculate bin centers
    bin_centers = (bins[:-1] + bins[1:]) / 2

    # Initial guesses for the parameters (nu, amplitude, b, scale)
    initial_guess = [8.4, 1.0, 1.0, np.mean(bin_centers)]

    # Fit the Rice distribution to the histogram data using curve_fit
    params, covariance = curve_fit(rice_pdf, bin_centers, hist_data, p0=initial_guess)
    nu, amplitude, b, scale = params

    # Create the fitted Rice distribution
    fitted_distribution = rice(nu, loc=scale, scale=np.sqrt(b**2 + scale**2))

    return fitted_distribution, nu, amplitude, b, scale

# Example usage:
if __name__ == "__main__":
    # Parameters for the Rice distribution
    nu = 8.5
    sigma = 10.5
    sample_size = 100

    # Calculate b from nu and sigma
    b = nu / sigma

    # Generate random data points from the Rice distribution
    data = rice.rvs(b=b, scale=sigma, size=sample_size)

    # Create a histogram of the generated data
    hist_data, bins, _ = plt.hist(data, bins=20, density=True, alpha=0.5, label="Generated Data")
    plt.xlabel("Value")
    plt.ylabel("Probability Density")

    # Fit a Rice distribution to the histogram using curve_fit
    fitted_distribution, fitted_nu, amplitude, fitted_b, fitted_scale = fit_rice_distribution_to_histogram(hist_data, bins)

    # Plot the original histogram and the fitted distribution
    x = np.linspace(min(bins), max(bins), 1000)
    pdf_values = fitted_distribution.pdf(x)
    plt.plot(x, pdf_values, 'r', label="Fitted Rice Distribution")
    plt.legend()
    plt.show()

    # Print fitted parameters
    print("Fitted Nu:", fitted_nu)
    print("Fitted Amplitude:", amplitude)
    print("Fitted b:", fitted_b)
    print("Fitted Scale:", fitted_scale)

Answer 1

根据您提供的试验数据集：

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats, optimize

nu = 8.5
sigma = 10.5
n = 100
b = nu / sigma

np.random.seed(12345)

# Reference Distribution:
X = stats.rice(b=b, scale=sigma)


# Binned dataset:
data = X.rvs(size=n)
density, bins = np.histogram(data, density=1.)
centers = (bins[:-1] + bins[1:]) / 2

从那里我们可以从直方图数据估计参数：

def model(x, b, loc, scale):
    return stats.rice.pdf(x, b=b, loc=loc, scale=scale)

popt, pcov = optimize.curve_fit(model, centers, density)

模型函数与您的模型函数显着不同，因为它缺少文档中指定的贝塞尔函数项。

xlin = np.linspace(0, 50, 200)
fig, axe = plt.subplots()
axe.hist(data, alpha=0.5, density=1.0, label="Data")
axe.plot(xlin, X.pdf(xlin), label="Model")
axe.plot(xlin, model(xlin, *popt), label="Fit")
axe.legend()
axe.grid()