Аппроксимация поверхности МНК

@Skaldy · Регистрация: 06.04.2013

Студворк — интернет-сервис помощи студентам

Доброго времени суток, уважаемые форумчане!
Передо мной стоит задача аппроксимации экспериментальных данных, представленных в виде матрицы nxn (грубо говоря, это изображение, размером nxn пикселей, значение каждого элемента матрицы - интенсивность изображения в данной точке). Необходимо аппроксимировать эти данные при помощи некоторой функции, т.е. методом наименьших квадратов определить коэффициенты, входящие в эту функцию для лучшего описания ею экспериментальных данных.
Как решать такую задачу при зависимости у(x) на форуме не раз подробно обсуждалось, да и в хелпе к матлабу понятно описано. А вот как быть, если случай трехмерный z(x,y) и моя функция (которой необходимо аппроксимировать экспериментальные данные) не банальный многочлен, а достаточно сложная отдельно вычисляемая функция в которую необходимые мне коэффициенты входят нелинейно.
Буду благодарна за любую помощь!

@Skaldy · 09.04.2013, 11:47 **[ТС]**

Скажите, пожалуйста, Curve Fitting Tool позволяет задавать для аппроксимации только те функции, которые там уже "вшиты" (типа гауссианов, экспонент, степенных многочленов) либо предлагает ввести CustomEquation вручную. Вопрос вот в чем: возможно ли вводить такую функцию не вручную, а задавать уже написанную в матлабе функцию (может, ее как-то определенно задать можно, чтобы было понятно где у нее коэффициенты для подгонки?).

@sergsh · 10.04.2013, 09:50

Сообщение от Skaldy

Доброго времени суток, уважаемые форумчане!
Передо мной стоит задача аппроксимации экспериментальных данных, представленных в виде матрицы nxn (грубо говоря, это изображение, размером nxn пикселей, значение каждого элемента матрицы - интенсивность изображения в данной точке)

А вот как быть, если случай трехмерный z(x,y)

хочу уточнить задачу

все же матрица N * M
или матрица N * M * K

матрица двумерная или трехмерная ?

@Skaldy · 10.04.2013, 09:56 **[ТС]**

Матрица двумерная NxN

@Skaldy · 10.04.2013, 10:06 **[ТС]**

Вот изображение моей поверхности.

@Skaldy · 10.04.2013, 10:16 **[ТС]**

Проблемма в том, что функцию, которой нужно аппроксимировать эти данные, невозможно просто написать ручками в CurveFitTool, её приходится рассчитывать в отдельном .m файле. В связи с этим вопрос: можно ли все-таки сделать это в CurveFitTool каким-то образом дав ссылку на мою функцию, зависящую от х,у и кучи параметров?

@Skaldy · 13.04.2013, 08:13 **[ТС]**

Я так поняла, что cftool мне использовать не удастся.
Тогда, взгляните, пожалуйста, можно ли так реализовать МНК:

Matlab M
1
2
3
4
5
6
7
function D=sqsum(theta0)
x=(1:74)*0.01; % 0.01 - разрешение картинки,74х74 - размер картинки в пикселях 
y=(1:74)*0.01;
F=preposefun(x,y,theta0); % - вычисляет функцию, которой буду аппроксимировать
data=importdata('data74x74.mat'); % - считывает матрицу данных
d=(data-F).^2; % - вычисляет квадраты отклонений в каждой точке
D=sum(sum(d,2),1); % - вычисляет сумму квадратов отклонений

После чего я просто минимизирую полученную разность D при помощи fminsearch.

@sergsh · 13.04.2013, 09:17

Все же не понятно что нужно.
theta0 - на вход функции откуда берется ?
Вы ее вычисляете ? Или это массив данных ?

@Skaldy · 14.04.2013, 19:59 **[ТС]**

Как видно на прикрепленной выше картинке, моя поверхность - это некий набор "холмов" различной высоты. В начальном приближении theta0 я задаю координаты х и у максимумов этих пиков, их высоты, полуширины и константу фона. theta0 - это вектор, количество элементов которого зависит от количества пиков на экспериментальной картинке. Я задаю его приблизительно, анализируя имеющуюся у меня картинку.

@sergsh · 14.04.2013, 22:58

Как можно задать "приблизительно" ?
Не понятно, можете привести пример приблизительного задания ?
Это будет число - ?

Можете привести весь код, который рисует вашу картинку и где вы определяете theta0 ?

@Skaldy · 15.04.2013, 00:49 **[ТС]**

Видимо, я не очень понятно выразилась. Не приблизительно, конечно.
Вот моя тета. Я задаю ее, подбирая так, чтобы получаемое с помощью моей функции изображение как можно больше соответствовало экспериментальной картинке, которую я приводила выше.
teta0=[0,01 0,37 0,74 0,01 0,38 0,74 0,01 0,38 0,74 0,19 0,56 0,19 0,56 0,01 0,01 0,01 0,38 0,38 0,38 0,74 0,74 0,74 0,2 0,19 0,56 0,56 40 40 40 40 40 40 40 40 40 55 55 55 55 0,04 0,04 0,04 0,04 0,04 0,04 0,04 0,04 0,04 0,07 0,07 0,07 0,07 10]
Картинка, получаемая при помощи моей функции:

Моя задача, найти такую тета, при которой конечная картинка будет максимально соответствовать экспериментальной. Для этого я хочу использовать МНК. Т.е. я задаю начальную точку theta0 (насколько могу сама ее подобрать), а дальше запускаю минимизацию функции разности между экспериментальной картинкой и моей расчетной.

@sergsh · 15.04.2013, 06:18

А экспериментальная картинка - это набор данных ?
Или вы просто сравниваете два рисунка ?

А как тогда вы определяете их совпадение - на глазок ?

Или вы сравниваете ДВА набора данных , расчетный и экспериментальный ?

@Skaldy · 15.04.2013, 09:58 **[ТС]**

Да, экспериментальная картинка - это набор данных. И я сравниваю два набора данных.

Добавлено через 2 минуты
Мне нужно чтобы каждый элемент расчетной матрицы как можно меньше отличался от экспериментальной и знать вектор тета при котором это происходит.

@sergsh · 15.04.2013, 11:28

Приведите тут если можно набор данных и ваш листинг по расчетным данным

@Skaldy · 15.04.2013, 14:18 **[ТС]**

Эксперимент и вычисления:
calculation.zip

experiment.zip

@sergsh · 16.04.2013, 23:07

В таблицах массивы данных, понятно.
А как вы получаете расчетный массив данных ?
Программа расчета доступна, она сделана в матлаб ?

Добавлено через 42 минуты
Импортировал ваши данные в матлаб
и вычел из экспериментальных данных расчетные, получилась такая картинка

Невязка имеет плавную структуру

Ее можно подобрать сплайнами и добавлять к расчетным данным в каждой точке сетки

Получим идеальное совпадение с наблюдением ...

@Skaldy · 17.04.2013, 11:10 **[ТС]**

Функция на матлабе. Боюсь, без пояснения физики процесса она будет казаться весьма странной.

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
function F=preposefun1(x,y,theta)
P=xlsread('P21x21.xls');
%figure (1), mesh(P);
 
m=20; % кол-во точек, на которые делим электронное пятно
xp=-0.1:0.2/m:0.1;
yp=-0.1:0.2/m:0.1;
 
[K,L]=meshgrid(x,y);
Mtot=13; % общее кол-во столбцов тета, равное общему кол-ву всех пиков 
O=0;
  for m1=1:Mtot;
      Sm=theta(2*Mtot+m1)*exp((-(K-theta(m1)).^2-(L-theta(Mtot+m1)).^2)/(2*(theta(3*Mtot+m1))^2));
      O=O+Sm;
  end
  O=O+theta(4*Mtot+1); % добавляем константу фона
%figure (2), mesh(x,y,O)
 
% граничные условия - окружаю матрицу О такими же матрицами:
O2x1=cat(1,O,O); O3x1=cat(1,O,O2x1); O3x2=cat(2,O3x1,O3x1); O3x3=cat(2,O3x2,O3x1);
 
% формируем набор отдельных матриц О по количеству пикселей нужного нам 
% сканирующего пятна mxm
Okl=[];
for l=1:74 % l - это строки
    O2=[];
    for k=1:74 % k - это столбцы
        O1=O3x3((74-m/2+l):(74+m/2+l), (74-m/2+k):(74+m/2+k));
        O2=cat(2,O2,O1);
    end
    Okl=cat(1,Okl,O2);
end
 
for l=1:74
    for k=1:74
        fkl(l,k)=trapz(xp,trapz(yp,Okl(1+(m+1)*(l-1):(m+1)*l,1+(m+1)*(k-1):(m+1)*k).*P));
    end
end
 
F=fkl;
%figure(3), mesh(x,y,F)

@Skaldy · 17.04.2013, 11:24 **[ТС]**

Не тот файл присоединила вот матрица P:
P21x21.xls
Прошу админов удалить аттач в предыдущем моем сообщении - он не нужен.

@Skaldy · 17.04.2013, 13:12 **[ТС]**

На всякий случай, попробую пояснить что эта функция делает.

Эксперимент - это сканирующая просвечивающая электронная микроскопия. Есть электронный пучок, который освещает объект исследования (образец) и есть изображение этого объекта, которое мы видим на экране (та самая матрица данных).
В программе P - это заранее вычисленное распределение интенсивности сфокусированного электронного "пятна" в плоскости образца XY.
Образец - это кристалл, т.е. периодическая структура. Эту структуру мы не знаем, но предполагаем, что ее можно задать как сумму гауссианов, для которых мы задаем вектор параметров тета, перечисляющий последовательно все координаты х максимумов, все координаты у максимумов, все высоты, все полуширины, и добавляем константу фона. При помощи тета вычисляем матрицу O.

Теперь мы имеем матрицу P размером 21x21 и матрицу O размером 74x74. Каждый элемент матрицы F(l,k) получается как интеграл от произведения матрицы P на (вот тут внимание!) часть матрицы O c центром на элементе O(l,k) и размером 21х21 (как P). Т.е. чтобы, например, вычислить F(12,12) я беру интеграл от P.*O(12-10:12+10,12-10:12+10).
Чтобы это вычислить я делаю два дополнительных шага: 1. добавляю граничные условия, т.е. окружаю матрицу O такими же матрицами со всех сторон, чтобы иметь возможность вычисления, например F(1,1); 2. заранее формирую новую матрицу Okl в которую стыкую все матрицы, которые будут перемножаться с P (этого можно и не делать, можно под интегралом задавать с какого по какой элемент надо брать O, но мне показалось, так удобнее).

Смысл всего этого действа: найти неизвестные параметры исследуемого объекта, записанные в векторе тета, зная параметры прибора (при помощи них вычисляется P) и матрицу получаемых данных.
По-этому, просто добавить невязку к каждому элементу нельзя. И начальное значение тета я тоже просто подбираю.

@sergsh · 17.04.2013, 16:48

Может я не очень понял, но идея у меня такая - задать больший диапазон изменения для ТЕТА
и для каждого значения ТЕТА вычислить величину невязки.

Тогда можно будет увидеть при каких значениях ТЕТА наступает минимум невязки ?

И в районе этого минимума невязки сделать много малых изменений ТЕТА для уточнения этого минимума.

Хотя экспериментальные данные с теорией никогда на 99% не сойдутся , ошибка в 30-50 % это хороший показатель.
Потому что теория отражает только часть факторов влияния, а эксперимент фиксирует ВСЕ , в том числе и НЕ ИЗВЕСТНЫЕ еще факторы

Новые блоги и статьи Все статьи Все блоги /
Отчёт о спецтехнике находящейся в ремонте Maks 20.04.2026 Отчёт из решения ниже размещен в конфигурации КА2. Задача: отобразить спецтехнику, которая на данный момент находится в ремонте. Есть нетиповой документ "Заявка на ремонт спецтехники" который. . .	Памятка для бота и "визитка" для читателей "Semantic Universe Layer (Слой семантической вселенной)" Hrethgir 19.04.2026 Сгенерировано для краткого описания по случаю сборки и компиляции скелета серверного приложения. И пусть после этого скажут, что статьи сгенерированные AI - туфта и не интересно. И это не реклама -. . .	Запрет удаления строк ТЧ документа при определенном условии Maks 19.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "Аккумуляторы", разработанного в конфигурации КА2. У данного документа есть ТЧ, в которой в зависимости от прав доступа. . .	Модель заражения группы наркоманов alhaos 17.04.2026 Условия задачи сформулированы тут Суть: - Группа наркоманов из 10 человек. - Только один инфицирован ВИЧ. - Колются одной иглой. - Колются раз в день. - Колются последовательно через. . .
Мысли в слух. Про "навсегда". kumehtar 16.04.2026 Подумалось тут, что наверное очень глупо использовать во всяких своих установках понятие "навсегда". Это очень сильное понятие, и я только начинаю понимать край его смысла, не смотря на то что давно. . .	My Business CRM MaGz GoLd 16.04.2026 Всем привет, недавно возникла потребность создать CRM, для личных нужд. Собственно программа предоставляет из себя базу данных клиентов, в которой можно фиксировать звонки, стадии сделки, а также. . .	Знаешь почему 90% людей редко бывают счастливыми? kumehtar 14.04.2026 Потому что они ждут. Ждут выходных, ждут отпуска, ждут удачного момента. . . а удачный момент так и не приходит.	Фиксация колонок в отчете СКД Maks 14.04.2026 Фиксация колонок в СКД отчета типа Таблица. Задача: зафиксировать три левых колонки в отчете. Процедура ПриКомпоновкеРезультата(ДокументРезультат, ДанныеРасшифровки, СтандартнаяОбработка) / / . . .

@Skaldy 0 / 0 / 0 Регистрация: 06.04.2013 Сообщений: 24

	Аппроксимация поверхности МНК 07.04.2013, 22:08. Показов 12716. Ответов 27 Метки нет (Все метки) Доброго времени суток, уважаемые форумчане! Передо мной стоит задача аппроксимации экспериментальных данных, представленных в виде матрицы nxn (грубо говоря, это изображение, размером nxn пикселей, значение каждого элемента матрицы - интенсивность изображения в данной точке). Необходимо аппроксимировать эти данные при помощи некоторой функции, т.е. методом наименьших квадратов определить коэффициенты, входящие в эту функцию для лучшего описания ею экспериментальных данных. Как решать такую задачу при зависимости у(x) на форуме не раз подробно обсуждалось, да и в хелпе к матлабу понятно описано. А вот как быть, если случай трехмерный z(x,y) и моя функция (которой необходимо аппроксимировать экспериментальные данные) не банальный многочлен, а достаточно сложная отдельно вычисляемая функция в которую необходимые мне коэффициенты входят нелинейно. Буду благодарна за любую помощь! 0

@Skaldy 0 / 0 / 0 Регистрация: 06.04.2013 Сообщений: 24
	09.04.2013, 11:47 [ТС]
	Скажите, пожалуйста, Curve Fitting Tool позволяет задавать для аппроксимации только те функции, которые там уже "вшиты" (типа гауссианов, экспонент, степенных многочленов) либо предлагает ввести CustomEquation вручную. Вопрос вот в чем: возможно ли вводить такую функцию не вручную, а задавать уже написанную в матлабе функцию (может, ее как-то определенно задать можно, чтобы было понятно где у нее коэффициенты для подгонки?). 0

@Skaldy 0 / 0 / 0 Регистрация: 06.04.2013 Сообщений: 24
	10.04.2013, 09:56 [ТС]
	Матрица двумерная NxN 0

@Skaldy 0 / 0 / 0 Регистрация: 06.04.2013 Сообщений: 24
	10.04.2013, 10:06 [ТС]
	Вот изображение моей поверхности. Миниатюры 0

@Skaldy 0 / 0 / 0 Регистрация: 06.04.2013 Сообщений: 24
	10.04.2013, 10:16 [ТС]
	Проблемма в том, что функцию, которой нужно аппроксимировать эти данные, невозможно просто написать ручками в CurveFitTool, её приходится рассчитывать в отдельном .m файле. В связи с этим вопрос: можно ли все-таки сделать это в CurveFitTool каким-то образом дав ссылку на мою функцию, зависящую от х,у и кучи параметров? 0

@sergsh 319 / 258 / 30 Регистрация: 30.03.2013 Сообщений: 755
	13.04.2013, 09:17
	Все же не понятно что нужно. theta0 - на вход функции откуда берется ? Вы ее вычисляете ? Или это массив данных ? 0

@Skaldy 0 / 0 / 0 Регистрация: 06.04.2013 Сообщений: 24
	14.04.2013, 19:59 [ТС]
	Как видно на прикрепленной выше картинке, моя поверхность - это некий набор "холмов" различной высоты. В начальном приближении theta0 я задаю координаты х и у максимумов этих пиков, их высоты, полуширины и константу фона. theta0 - это вектор, количество элементов которого зависит от количества пиков на экспериментальной картинке. Я задаю его приблизительно, анализируя имеющуюся у меня картинку. 0

@sergsh 319 / 258 / 30 Регистрация: 30.03.2013 Сообщений: 755
	14.04.2013, 22:58
	Как можно задать "приблизительно" ? Не понятно, можете привести пример приблизительного задания ? Это будет число - ? Можете привести весь код, который рисует вашу картинку и где вы определяете theta0 ? 0

@Skaldy 0 / 0 / 0 Регистрация: 06.04.2013 Сообщений: 24
	15.04.2013, 00:49 [ТС]
	Видимо, я не очень понятно выразилась. Не приблизительно, конечно. Вот моя тета. Я задаю ее, подбирая так, чтобы получаемое с помощью моей функции изображение как можно больше соответствовало экспериментальной картинке, которую я приводила выше. teta0=[0,01 0,37 0,74 0,01 0,38 0,74 0,01 0,38 0,74 0,19 0,56 0,19 0,56 0,01 0,01 0,01 0,38 0,38 0,38 0,74 0,74 0,74 0,2 0,19 0,56 0,56 40 40 40 40 40 40 40 40 40 55 55 55 55 0,04 0,04 0,04 0,04 0,04 0,04 0,04 0,04 0,04 0,07 0,07 0,07 0,07 10] Картинка, получаемая при помощи моей функции: Моя задача, найти такую тета, при которой конечная картинка будет максимально соответствовать экспериментальной. Для этого я хочу использовать МНК. Т.е. я задаю начальную точку theta0 (насколько могу сама ее подобрать), а дальше запускаю минимизацию функции разности между экспериментальной картинкой и моей расчетной. 0

@sergsh 319 / 258 / 30 Регистрация: 30.03.2013 Сообщений: 755
	15.04.2013, 06:18
	А экспериментальная картинка - это набор данных ? Или вы просто сравниваете два рисунка ? А как тогда вы определяете их совпадение - на глазок ? Или вы сравниваете ДВА набора данных , расчетный и экспериментальный ? 0

@Skaldy 0 / 0 / 0 Регистрация: 06.04.2013 Сообщений: 24
	15.04.2013, 09:58 [ТС]
	Да, экспериментальная картинка - это набор данных. И я сравниваю два набора данных. Добавлено через 2 минуты Мне нужно чтобы каждый элемент расчетной матрицы как можно меньше отличался от экспериментальной и знать вектор тета при котором это происходит. 0

@sergsh 319 / 258 / 30 Регистрация: 30.03.2013 Сообщений: 755
	15.04.2013, 11:28
	Приведите тут если можно набор данных и ваш листинг по расчетным данным 0

@Skaldy 0 / 0 / 0 Регистрация: 06.04.2013 Сообщений: 24
	15.04.2013, 14:18 [ТС]
	Эксперимент и вычисления: calculation.zip experiment.zip 0

@sergsh 319 / 258 / 30 Регистрация: 30.03.2013 Сообщений: 755
	16.04.2013, 23:07
	В таблицах массивы данных, понятно. А как вы получаете расчетный массив данных ? Программа расчета доступна, она сделана в матлаб ? Добавлено через 42 минуты Импортировал ваши данные в матлаб и вычел из экспериментальных данных расчетные, получилась такая картинка Невязка имеет плавную структуру Ее можно подобрать сплайнами и добавлять к расчетным данным в каждой точке сетки Получим идеальное совпадение с наблюдением ... 0

@Skaldy 0 / 0 / 0 Регистрация: 06.04.2013 Сообщений: 24
	17.04.2013, 11:24 [ТС]
	Не тот файл присоединила вот матрица P: P21x21.xls Прошу админов удалить аттач в предыдущем моем сообщении - он не нужен. 0

@Skaldy 0 / 0 / 0 Регистрация: 06.04.2013 Сообщений: 24
	17.04.2013, 13:12 [ТС]
	На всякий случай, попробую пояснить что эта функция делает. Эксперимент - это сканирующая просвечивающая электронная микроскопия. Есть электронный пучок, который освещает объект исследования (образец) и есть изображение этого объекта, которое мы видим на экране (та самая матрица данных). В программе P - это заранее вычисленное распределение интенсивности сфокусированного электронного "пятна" в плоскости образца XY. Образец - это кристалл, т.е. периодическая структура. Эту структуру мы не знаем, но предполагаем, что ее можно задать как сумму гауссианов, для которых мы задаем вектор параметров тета, перечисляющий последовательно все координаты х максимумов, все координаты у максимумов, все высоты, все полуширины, и добавляем константу фона. При помощи тета вычисляем матрицу O. Теперь мы имеем матрицу P размером 21x21 и матрицу O размером 74x74. Каждый элемент матрицы F(l,k) получается как интеграл от произведения матрицы P на (вот тут внимание!) часть матрицы O c центром на элементе O(l,k) и размером 21х21 (как P). Т.е. чтобы, например, вычислить F(12,12) я беру интеграл от P.*O(12-10:12+10,12-10:12+10). Чтобы это вычислить я делаю два дополнительных шага: 1. добавляю граничные условия, т.е. окружаю матрицу O такими же матрицами со всех сторон, чтобы иметь возможность вычисления, например F(1,1); 2. заранее формирую новую матрицу Okl в которую стыкую все матрицы, которые будут перемножаться с P (этого можно и не делать, можно под интегралом задавать с какого по какой элемент надо брать O, но мне показалось, так удобнее). Смысл всего этого действа: найти неизвестные параметры исследуемого объекта, записанные в векторе тета, зная параметры прибора (при помощи них вычисляется P) и матрицу получаемых данных. По-этому, просто добавить невязку к каждому элементу нельзя. И начальное значение тета я тоже просто подбираю. 0

@sergsh 319 / 258 / 30 Регистрация: 30.03.2013 Сообщений: 755
	17.04.2013, 16:48
	Может я не очень понял, но идея у меня такая - задать больший диапазон изменения для ТЕТА и для каждого значения ТЕТА вычислить величину невязки. Тогда можно будет увидеть при каких значениях ТЕТА наступает минимум невязки ? И в районе этого минимума невязки сделать много малых изменений ТЕТА для уточнения этого минимума. Хотя экспериментальные данные с теорией никогда на 99% не сойдутся , ошибка в 30-50 % это хороший показатель. Потому что теория отражает только часть факторов влияния, а эксперимент фиксирует ВСЕ , в том числе и НЕ ИЗВЕСТНЫЕ еще факторы 1