Уменьшенный хи-квадрат слишком мал (близок к 0) после взвешенной подгонки — интеграл свертки — Python lmfit

Я подгоняю следующие данные, где t: время (с), G: количество импульсов в секунду, f: импульсная функция (мм/с):

    t     G     f 
    0   4.58    0
  900   11.73   (11/900)
 1800   18.23   (8.25/900)
 2700   19.33   (3/900)
 3600   19.04   (0.5/900)
 4500   17.21   0
 5400   12.98   0
 6300   11.59   0
 7200   9.26    0
 8100   7.66    0
 9000   6.59    0
 9900   5.68    0
10800   5.1     0

Используя следующий интеграл свертки:

И, более конкретно:

Где: lambda_1 = 0.000431062 и lambda_2 = 0.000580525.

Код, используемый для выполнения этой подгонки:

#Extract data into numpy arrays
t=df['t'].as_matrix()
g=df['G'].as_matrix()
f=df['f'].as_matrix()
#add parameters
params=Parameters()
params.add('a',value=1)
params.add('b',value=0.7)
params.add('c',value =1)
#define functions
def exp(x,k):
    return np.exp(-x*k)

def residuals(params,x,y):
    A=params['a'].value
    B=params['b'].value
    C=params['c'].value

    dt=x[2]-x[1]
    model = A*(np.convolve(exp(x,lambda_1), f))[:len(x)]*dt+B*np.convolve(exp(x,lambda_2), f)[:len(x)]*dt+C
    weights=1/np.sqrt(y)
    return (model - y)*weights

#perform fit using leastsq
result = minimize(residuals, params, args=(t,g))
final = g + result.residual
print(report_fit(result))

Это работает, однако я получаю очень низкий приведенный хи-квадрат (около 0), когда я умножаю остаток, который нужно минимизировать, на вес (1/np.sqrt (g) (взвешенная подгонка). Если я не приму во внимание веса (невзвешенная подгонка), я получаю приведенный хи-квадрат 0,254. Я хотел бы получить уменьшенный хи-квадрат около 1.


person Gerard    schedule 13.08.2018    source источник


Ответы (1)


Уменьшенный хи-квадрат намного ниже 1 будет означать, что ваша оценка неопределенности данных слишком велика. Если я правильно понял ваш пример, вы используете квадратный корень из G в качестве неопределенности в G. Использование квадратного корня является стандартным подходом для оценки неопределенностей в значениях, в которых преобладает статистика подсчета.

Но... ваш G - это число с плавающей запятой, которое вы описываете как число в секунду. Я мог бы предположить, что счет в секунду превышает 900 секунд.

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

g_values = [4.58 ,  11.73,  18.23] 
g_uncertainties = sqrt(g_values) = [2.1401, 3.4249, 4.2697]

но неопределенность в подсчетах будет sqrt(g_values*900), и, таким образом, неопределенность в подсчетах в секунду на sqrt(g_values*900)/900 = sqrt(g_values)/30.

Более формально неопределенности в значении, представляющем «количество отсчетов за время», добавляли бы неопределенности в счетах и ​​неопределенности во времени в квадратуре. Но опять же, неопределенность вашего времени, вероятно, очень мала (или, по крайней мере, ваши временные данные подразумевают, что оно меньше 1 секунды).

person M Newville    schedule 16.08.2018
comment
спасибо, но делая то, что вы говорите, мой уменьшенный хи-квадрат становится очень высоким... общее количество не совсем равно g_values ​​* 900, у меня есть другой массив для общего количества, и я использовал его для вычисления неопределенностей моих данных... однако я получаю приведенный хи-квадрат приблизительно 5... Может быть, это связано с тем, что ожидаемое количество импульсов в секунду в основном меньше 5??? Спасибо - person Gerard; 27.08.2018
comment
Чтобы получить уменьшенный хи-квадрат около 1,0, вам нужно а) хорошее соответствие и б) точная оценка неопределенностей в данных, которые вы фактически моделируете. Я понятия не имею, над какой проблемой вы работаете, поэтому я не могу сказать вам, каковы погрешности в ваших данных, но если вы извлекаете квадратный корень из скорости счета, вы определенно не ошибаетесь. будет иметь точную оценку неопределенностей в ваших данных - это не проходит размерный анализ.... Неопределенности в y должны иметь те же единицы измерения, что и y. - person M Newville; 27.08.2018
comment
Я моделирую данные, связанные с активностью (количество в секунду) радионуклида (радиоактивный распад), которая соответствует подсчетам сетей (общие подсчеты - фоновые подсчеты). Мои оценочные неопределенности равны квадратному корню из общего количества (чистое количество + фоновое количество), деленному на 900. Я не извлекаю квадратный корень из скорости счета... и неопределенность в отношении моих фоновых подсчетов незначительна. Мой вопрос: может быть, уменьшенный хи-квадрат не является точным статистическим тестом для проверки качества моего соответствия (поскольку ожидаемые значения очень низкие, а размер моей выборки небольшой)? Спасибо - person Gerard; 28.08.2018
comment
Что касается примера кода в вашем вопросе, похоже, что вы извлекаете квадратный корень из скорости счета. У вас есть weight = 1./np.sqrt(y) в вашей целевой функции, и вы передаете столбец с надписью «G» в вашем файле как y. Как это что-нибудь другое из извлечения квадратного корня из скорости счета? - person M Newville; 28.08.2018