Как выполнить подгонку MCMC, используя только подмножество переменных функции? ⇐ Python

[Anonymous]

У меня есть функция, имеющая множество параметров/переменных. Для простоты можно предположить, что это f(a,b,c)=a + b * x + c * x**2, где параметры «a», «b» и «c» могут быть либо фиксированными параметрами, либо переменными. Я хочу подогнать.
Я без проблем могу подогнать f(a,b,c) под три параметра одновременно. Проблема в том, что я хочу исправить «b» и подогнать «a» и «c» или исправить «b» и «c» и просто подогнать «a». Мотивация заключается в том, что в физических системах некоторые из этих параметров могут быть известны, и поэтому я могу сосредоточиться на подгонке только других (или, может быть, у меня недостаточно точек данных для соответствия многим параметрам, и я предпочитаю исправить некоторые и оставить свободными только самые актуальные).
Это версия кода, которая работает, фиксируя все параметры:

Код: Выделить всё

import numpy as np
import emcee
import corner
import matplotlib.pyplot as plt

# True parameters of the function
true_a = 1.0
true_b = 2.0
true_c = 0.5

# Generate synthetic data with noise
np.random.seed(42)
x_data = np.linspace(-8, 10, 50)
y_true = true_a + true_b * x_data + true_c * x_data**2
y_data = y_true + np.random.normal(scale=0.05, size=len(x_data))

# Define the quadratic function
def f(x, a, b, c):
return a + b * x + c * x**2

# Define the log-likelihood function
def log_likelihood(params, x, y, y_err):
a, b, c = params
model = f(x, a, b, c)
return -0.5 * np.sum((y - model)**2 / y_err**2)

# Define the log-prior function
def log_prior(params):
# Uniform priors for a, b, c
if all(-10.0 < p < 10.0 for p in params):
return 0.0
return -np.inf

# Define the log-posterior function
def log_posterior(params, x, y, y_err):
lp = log_prior(params)
if not np.isfinite(lp):
return -np.inf
return lp + log_likelihood(params, x, y, y_err)

# Perform MCMC fitting
ndim = 3  # Number of parameters (a, b, c)
nwalkers = 50
nsteps = 2000
burnin = nsteps // 3  # Burn-in steps (1/3 of the total steps)

# Initialize walkers with random values around the true parameters
initial_params = [true_a, true_b, true_c]
per = 0.1
initial_pos = [initial_params + per * np.random.randn(ndim) for _ in range(nwalkers)]

# Set up the MCMC sampler
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=(x_data, y_data, 0.5))

# Run the sampler with burn-in
sampler.run_mcmc(initial_pos, nsteps, progress=True)

# Get the samples after burn-in
samples = sampler.get_chain(discard=burnin, flat=True)

Теперь это моя попытка определить, какие параметры подходят, используя в качестве аргумента словарь (хотя, возможно, это неправильный подход?)

Код: Выделить всё

import numpy as np
import emcee
import corner
import matplotlib.pyplot as plt

# True parameters of the function
true_a = 1.0
true_b = 2.0
true_c = 0.5

# Generate synthetic data with noise
np.random.seed(42)
x_data = np.linspace(-8, 10, 50)
y_true = true_a + true_b * x_data + true_c * x_data**2
y_data = y_true + np.random.normal(scale=0.05, size=len(x_data))

# Define the function
def f(x, params):
a = params.get('a', 0.0)
b = params.get('b', 0.0)
c = params.get('c', 0.0)
return a + b * x + c * x**2

# Define the log-likelihood function
def log_likelihood(params, x, y, y_err):
model = f(x, params)
return -0.5 * np.sum((y - model)**2 / y_err**2)

# Define the log-prior function
def log_prior(params):
# Uniform priors for a, b, c
a = params.get('a', 0.0)
b = params.get('b', 0.0)
c = params.get('c', 0.0)
if all(-10.0 < p <  10.0 for p in [a, b, c]):
return 0.0
return -np.inf

# Define the log-posterior function
def log_posterior(params, x, y, y_err):
lp = log_prior(params)
if not np.isfinite(lp):
return -np.inf
return lp + log_likelihood(params, x, y, y_err)

# Initialize walkers with random values around the true parameters
nwalkers = 50
initial_params = {'a': true_a, 'b': true_b, 'c': true_c} # here I would like to have any combination of parameters, not necessarily all three
per = 0.1
initial_pos = [{key: val + per * np.random.randn() for key, val in initial_params.items()} for _ in range(nwalkers)]

# Perform MCMC fitting
ndim = len(initial_params)  # Number of parameters
nsteps = 2000
burnin = nsteps // 3  # Burn-in steps (1/3 of the total steps)

# Set up the MCMC sampler
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=(x_data, y_data, 0.05))

# Run the sampler with burn-in
sampler.run_mcmc(initial_pos, nsteps, progress=True)

# Get the samples after burn-in
samples = sampler.get_chain(discard=burnin, flat=True)

Но при запуске он говорит:
---> 63 sampler.run_mcmc(initial_pos, nsteps, Progress=True)
ValueError: несовместимые входные размеры (1 , 50)
Я был бы очень признателен за отзывы о том, как решить эту проблему.

Подробнее здесь: https://stackoverflow.com/questions/774 ... a-function