如何在Python中实现自动模型判定及双态模型拟合？

Hey there! I’ve been in your exact situation—trading manual Prism fits for a Python workflow saves so much time and lets you scale your analysis easily. Let’s break down your two questions with practical, actionable code.

1. Implementing the Two-State (Linear + Plateau) Model #

The model you described—initial linear growth followed by a flat plateau—fits perfectly with a piecewise linear function. Here’s how to implement it using scipy.optimize.curve_fit, Python’s go-to tool for curve fitting:

First, define the mathematical form of the model:

import numpy as np

def two_state_model(x, a, b, x0, plateau):
    # a: slope of the linear growth phase
    # b: intercept of the linear phase
    # x0: inflection point where growth stops and plateau starts
    # plateau: value of the flat phase
    return np.where(x < x0, a * x + b, plateau)

Now, a full example with simulated data (swap this with your actual dataset):

from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Generate simulated data matching your model behavior
x_data = np.linspace(0, 20, 50)
y_true = two_state_model(x_data, 2.5, 1.0, 10, 26.0)
y_data = y_true + np.random.normal(0, 1.2, size=len(x_data))  # Add realistic noise

# Critical: Set initial guesses for parameters (matches data trends roughly)
initial_guess = [2, 0, 8, 25]

# Fit the model to your data
params, covariance = curve_fit(two_state_model, x_data, y_data, p0=initial_guess)

# Extract fitted parameters
a_fit, b_fit, x0_fit, plateau_fit = params

# Visualize results
plt.scatter(x_data, y_data, label='Raw Data')
plt.plot(x_data, two_state_model(x_data, *params), 'r-', linewidth=2, 
         label=f'Fitted Model:\nSlope={a_fit:.2f}, Intercept={b_fit:.2f}\nInflection={x0_fit:.2f}, Plateau={plateau_fit:.2f}')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.show()

Pro tip: If fitting fails to converge, tweak your initial guesses to be closer to your data’s actual trends. You can also set parameter bounds with the bounds argument in curve_fit (e.g., ensure the plateau value is positive).

2. Automating Model Selection #

To automatically pick the best model for your data, compare multiple candidate models using metrics that balance fit quality and model complexity. The most reliable metrics are:

• R²: Measures variance explained (closer to 1 = better fit)
• AIC (Akaike Information Criterion): Penalizes overly complex models (smaller = better)
• BIC (Bayesian Information Criterion): Similar to AIC but harsher on complex models

Here’s a step-by-step workflow:

First, define candidate models (we’ll include linear, two-state, and logistic as alternatives):

def linear_model(x, m, c):
    return m * x + c

def logistic_model(x, L, x0, k, b):
    # L: plateau value, x0: growth midpoint, k: growth rate, b: intercept
    return L / (1 + np.exp(-k*(x - x0))) + b

Next, write a helper function to fit models and compute evaluation metrics:

def evaluate_model(model, x, y, initial_guess):
    try:
        params, _ = curve_fit(model, x, y, p0=initial_guess)
        y_pred = model(x, *params)
        
        # Calculate R²
        ss_res = np.sum((y - y_pred)**2)
        ss_tot = np.sum((y - np.mean(y))**2)
        r_squared = 1 - (ss_res / ss_tot)
        
        # Calculate AIC and BIC
        n = len(y)
        k = len(params)
        aic = n * np.log(ss_res / n) + 2 * k
        bic = n * np.log(ss_res / n) + k * np.log(n)
        
        return {
            'model_name': model.__name__,
            'params': params,
            'r_squared': r_squared,
            'aic': aic,
            'bic': bic
        }
    except RuntimeError:
        # Skip models that fail to fit
        return None

Finally, run evaluations and select the best model:

# Define initial guesses for each model (tailor these to your data!)
model_guesses = {
    linear_model: [2, 0],
    two_state_model: [2, 0, 8, 25],
    logistic_model: [26, 10, 0.5, 1]
}

# Evaluate all candidate models
results = []
for model, guess in model_guesses.items():
    res = evaluate_model(model, x_data, y_data, guess)
    if res:
        results.append(res)

# Sort models by AIC (smallest = best balance of fit and simplicity)
results_sorted = sorted(results, key=lambda x: x['aic'])

# Print the best model details
best_model = results_sorted[0]
print(f"Best Model: {best_model['model_name']}")
print(f"R² Score: {best_model['r_squared']:.4f}")
print(f"AIC Value: {best_model['aic']:.2f}")
print(f"Fitted Parameters: {[round(p, 2) for p in best_model['params']]}")

You can add more candidate models (like exponential growth) if you think they might apply to your dataset.

内容的提问来源于stack exchange，提问作者Isaac Wong