Начальные шаги
Я думаю, мы должны начать с представления в Matlab функции, которую вы пытаетесь смоделировать. Прямой перевод вашей формулы выглядит так:
function y = targetfunction(A,V,P,bc,b)
y = (A/2) * (1 - V * sin((b-bc) / P));
end
Получение данных
Следующим моим шагом будет создание данных для работы (естественно, вы будете использовать свои собственные данные). Итак, вот функция, которая генерирует зашумленные данные. Обратите внимание, что я указал некоторые значения параметров.
function [y b] = generateData(npoints,noise)
A = 2;
V = 1;
P = 0.7;
bc = 0;
b = 2 * pi * rand(npoints,1);
y = targetfunction(A,V,P,bc,b) + noise * randn(npoints,1);
end
Функция rand
генерирует случайные точки на интервале [0,1]
, и я умножил их на 2*pi
, чтобы случайным образом получить точки на интервале [0, 2*pi]
. Затем я применил целевую функцию в этих точках и добавил немного шума (функция randn
генерирует нормально распределенные случайные величины).
Подходящие параметры
Самая сложная функция — это та, которая подгоняет модель к вашим данным. Для этого я использую функцию fminunc
, выполняющую неограниченную минимизацию. Рутина выглядит так:
function [A V P bc] = bestfit(y,b)
x0(1) = 1; %# A
x0(2) = 1; %# V
x0(3) = 0.5; %# P
x0(4) = 0; %# bc
f = @(x) norm(y - targetfunction(x(1),x(2),x(3),x(4),b));
x = fminunc(f,x0);
A = x(1);
V = x(2);
P = x(3);
bc = x(4);
end
Пройдёмся по строчке. Во-первых, я определяю функцию f
, которую я хочу минимизировать. Это не слишком сложно. Чтобы минимизировать функцию в Matlab, ей нужно взять в качестве параметра один вектор. Поэтому мы должны упаковать наши четыре параметра в вектор, что я и делаю в первых четырех строках. Я использовал значения, которые близки, но не совпадают с теми, которые я использовал для генерации данных.
Затем я определяю функцию, которую хочу минимизировать. Он принимает один аргумент x
, который распаковывается и передается в targetfunction
вместе с точками b
в нашем наборе данных. Надеюсь, они близки к y
. Мы измеряем, насколько они далеки от y
, вычитая из y
и применяя функцию norm
, которая возводит в квадрат каждый компонент, складывает их и извлекает квадратный корень (т. е. вычисляет среднеквадратичную ошибку).
Затем я вызываю fminunc
с минимизацией нашей функции и начальным предположением параметров. При этом используется внутренняя процедура для поиска ближайшего совпадения для каждого из параметров и возврата их в векторе x
.
Наконец, распаковываю параметры из вектора x
.
Собираем все вместе
Теперь у нас есть все необходимые компоненты, поэтому нам нужна только одна последняя функция, которая свяжет их вместе. Вот:
function master
%# Generate some data (you should read in your own data here)
[f b] = generateData(1000,1);
%# Find the best fitting parameters
[A V P bc] = bestfit(f,b);
%# Print them to the screen
fprintf('A = %f\n',A)
fprintf('V = %f\n',V)
fprintf('P = %f\n',P)
fprintf('bc = %f\n',bc)
%# Make plots of the data and the function we have fitted
plot(b,f,'.');
hold on
plot(sort(b),targetfunction(A,V,P,bc,sort(b)),'r','LineWidth',2)
end
Если я запускаю эту функцию, я вижу, что это печатается на экране:
>> master
Local minimum found.
Optimization completed because the size of the gradient is less than
the default value of the function tolerance.
A = 1.991727
V = 0.979819
P = 0.695265
bc = 0.067431
И появляется следующий сюжет:
![введите здесь описание изображения](https://i.stack.imgur.com/cwjID.png)
Эта посадка мне кажется достаточно хорошей. Дайте мне знать, если у вас есть какие-либо вопросы о том, что я здесь делал.
person
Chris Taylor
schedule
02.11.2012