Регрессия на главные компоненты для анализа драгметаллов

[McConst]

Евгений Машеров
Спасибо большое. Я ещё из предыдущего сообщения понял, что выбор PCR не самый удачный. В любом случае за PLS буду браться, уже есть наработки, но в варианте нормирования это займет несколько месяцев. PCR - это первый шаг, от которого ожидал результат сразу. Реальность оказалась более сурова, но Andrey_Kireew написал, что для моего набора данных есть шансы и с PCR не всё безнадежно. Если что-то получится - информацию можно использовать для PLS. Например, интенсивности для нормирования на внутренний стандарт я хочу брать не от прибора, а вычислять из спектра. Для того, чтобы иметь матрицы для предсказания вектора нормирования, я должен их найти используя начальное приближение. Если брать начальное приближение от прибора, то для подбора коэффициентов PLS из-за согласованности X и Y потребуется очень много шагов. Если PCR даст более точный вектор для нормирования градуировочной матрицы, его можно использовать как начальное приближение для PLS, чтобы сократить объем вычислений. Может быть туманно выразился. Про ридж-регрессии, я к сожалению с ними не знаком. Вижу в вашей формуле дополнительные параметры, но чтобы вникнуть в их смысл, нужно плотнее уйти в литературу - на это требуется время. Мне сегодня один знакомый написал про sparse PLS - вводится дополнительное слагаемое, которое штрафует модель за использование большого числа каналов в одном компоненте. Я пока не вникал что там за слагаемое, но в голове кусочек вашей формулы всплыл $(X^TX+kI)$ - тоже слагаемое появляется. Может быть о разных вещах речь идет, но всё что вы говорили про PCR наверняка можно будет использовать для PLS? Вернусь к этому позже. Сейчас у меня пока перегруз по информации. Уверен, в dxdy ещё не раз придется зависнуть.

Emergency
Мне пока это видится большим объемом работы - именно правильное вычитание. Нужно много экспериментальных данных по поведению линий, которых у меня пока нет. И линии, которые нужно вычитать, они тоже перекашиваются другими линиями. Там одна линия на другую завязаны. Допустим, между двумя линиями одного элемента появляется линия нового элемента. Та базовая линия, по которой будем расчитывать перекрытый спектр, она не изменилась, так как более высокой энергии, а та, которую будем вычитать, она уже изменилась промежуточной линией. Разработчики прибора всё это вычитают согласно физическим моделям. Я пока к таким шагам не готов, так как вижу перспективу в непосредственной работе с исходным спектром. Метод более универсальный, не привязан к физике моего анализа, а вообще применим к любой области спектроскопии: ИК, АЭС, ААС, впрочем - не только в оптике - полярография, хроматография - везде.

igor_ivanov

Цитата:

В Вашем уравнении $n$ - индекс элемента в списке, $i$ - номер канала?

Нет $i$ - это номер уравнения в ММНК при расчете по линиям для хим элемента x. Это уравнение не для поканального расчета, а для расчета по линиям. Ну, вроде того, что сумма каналов в области 328,1 нм даст интенсивность линии Ag в ультрафиолете. Ссылку на это уравнение я давал, когда вы спрашивали, есть ли оно в литературе.

[Евгений Машеров]

$(X^TX+kI)$ это как раз введение штрафа за величину коэффициентов. Величина штрафа определяется значением k, строгие критерии выбора мне неизвестны. Чем больше k, тем сильнее сокращается дисперсия коэффициентов, но при этом они перестают быть несмещёнными, их матожидание смещается к нулю. Что, насколько я понимаю, для данной задачи не столь существенный недостаток, занизить оценку - слегка наколоть "сдатчика", завысить - потерять деньги. Хотя мой подход лично мне нравится больше.

[McConst]

Евгений Машеров
Такая же запись используется при поиске минимума в методе Левенберга-Марквардта. В книгах пишут, что величиной константы регулируется шаг, с которым сходится решение. Я эту константу подбирал эмпирически - придумываю начальное значение как на голову ляжет, потом ввожу 8 параллельных констант, одна в два раза больше, остальные в 2, 4, 8, 16... раз меньше. Каждому ядру процессора задаю параллельный расчет одной из констант и в качестве ответа выбираю ту, при которой сходимость лучше.
В каких адекватных границах должна быть данная константа? Она положительна, или может принимать отрицательные значения? Абсолютные значения в пределах 0...1 или потолка нет?

Цитата:

занизить оценку - слегка наколоть "сдатчика", завысить - потерять деньги.

Если "занизить" до нуля - тогда уравнение принимает классический вид для расчета коэффициентов регрессии и проблем для сдатчика вроде как не должно быть. Если можно, поясните в 2-3 предложения.

[Emergency]

McConst в сообщении #1436269 писал(а):

вижу перспективу в непосредственной работе с исходным спектром

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

По аналогии с фотоаппаратом это либо формат RAW, содержащий необработанные данные об электрических сигналах с фотоматрицы, либо формат преобразованный в RGB. Из RAW можно извлечь больше, но и более трудоемко.

[Евгений Машеров]

Ну, собственно, в нелинейном оценивании Л-М смысл тот же. Вводится штраф за большие коэффициенты (собственно, у нас теперь минимизируется не сумма квадратов отклонений модели от наблюдений, а эта сумма плюс сумма квадратов коэффициентов, умноженных на k. Этим резко уменьшается дисперсия коэффициентов, так что случаев "перескакивания" будет меньше, но зато сами оценки вспомогательной линейной регрессии (отклонений на производные по параметрам модели) уменьшаются, а поскольку они дают поправки для нелинейной модели, приближение замедляется. То есть малый коэффициент - идём у цели быстро, но постоянно перескакиваем, большой - медленно и надёжно ползём (при $k\rightarrow\infty$ Левенберг-Макрвардт переходит в метод градиента). k неотрицательно, отрицательное k будет означать, что нам желательно увеличить абсолютные величины коэффициентов, даже ценой бессмысленных значений некоторых (бессмысленно больших или неадекватных по знаку, впрочем, можно и гармонично сочетать). Для линейной регрессии выбор оптимального k неясен. Чем оно больше, тем меньше значения оценок параметров (при $k\rightarrow\infty$ будет 0, величина с гарантированно нулевой дисперсией, но как-то малополезная), при нулевом k оценка несмещённая, но дисперсия может быть весьма велика. Предлагали строить график значений параметров в зависимости от k, и по нему интуитивно выбирать. Но как-то стрёмно...
Что до моего не вполне серьёзного замечания, то речь идёт об оценках содержания металлов, которые, если я правильно понимаю, получаются, как коэффициенты регрессии. И если взять завышенное k, то оценки содержания металлов будут занижены, что выгодно покупателю этого сырья, но невыгодно сдатчику. А при малом k (особенно нулевом) дисперсии коэффициентов велики, и за счёт этого у одних содержание будет завышено, а у других будут физически бессмысленные отрицательные коэффициенты, что убыточно для покупателя.

[Emergency]

Разве нельзя программно запретить отрицательные концентрации?

[McConst]

Цитата:

график значений параметров в зависимости от k, и по нему интуитивно выбирать. Но как-то стрёмно...

Тут нужно k подбирать так, чтобы не было вопросов у контролирующих органов. Интуитивно не пойдёт. Возможно следует собрать килограмм 100 материала, сделать итоговый анализ и исходя из того, перебор или недобор, сделать коррекцию с помощью k. В итоге методом постепенного приближения к ответу его можно подобрать. Интересная тема. Потом в литературе поищу. Без ссылок на документы такие вещи лучше не экспериментировать, особенно как у меня - с поверхностным представлением как и что работает.
Emergency

Цитата:

Разве нельзя программно запретить отрицательные концентрации?

Так и сделал. Но разброс всё портит. Способы решения мне тут уже подсказали, нужно теперь пробовать варианты.

Возможно я ещё одну ошибку при калибровке сделал. Я прописал одинаковое количество ГК для всех хим. элементов. Так получилось в примере на сайте хемометриков Там об этом не писалось, что всегда так должно быть. Я это из логики сделал, что независимые переменные с правильным значением ГК наиболее полно описывают всю систему, остальное - шум. Сейчас подумал, что в неидеальной системе так наверное нельзя. Для каждого элемента лучше использовать свои количества ГК. Особенно для PCR, в котором с ростом ГК вклад шума увеличивается. А ещё для моей реализации с нормирования получается слишком много параметров. Параметры для расчета коэффициентов нормирования и косвенно связанные с ними параметры непосредственной градуировки. Возможно для системы условно из 20 компонент 64 образца - это мало, нужно снижать размерность как минимум через PLS.

[Евгений Машеров]

Emergency в сообщении #1436370 писал(а):

Разве нельзя программно запретить отрицательные концентрации?

Можно. Но это значит, что мы вводим ограничения-неравенства, и задача превращается в задачу квадратичной оптимизации, существенно более сложную. Или в перебор всех подмножеств регрессоров. Требующий рассмотрения $2^n$ моделей, из которых выбирается наилучшая удовлетворяющая ограничениям на знаки.

[Евгений Машеров]

Почему мне не нравится идея регрессии на главные компоненты. Вот простой пример, в котором она отчётливо не работает (и априори узнать, тот ли случай, когда она полезна, или же не работает - не могу).
Зависимая переменная равна разности двух регрессоров, слегка отягощённой шумом $y=x_1-x_2+\varepsilon$, причём регрессоры сильно положительно коррелированы (для конкретизации примера - исследуется зависимость прибыли некоего предприятия от дохода и расхода, меняющихся почти синхронно в зависимости от колебаний объёма продаж, то есть правильные коэффициенты $+1$ и $-1$, но исследователю доступны только сами данные, а не экономическая их интерпретация, позволившая бы обойтись без статистического анализа).
В численном примере корреляция между регрессорами 99.6%, между регрессорами и регрессандом 2.5% и 6.8% соответственно. То есть отбрасывание одного из регрессоров даёт очень плохо объясняющую модель. Включение обоих даёт коэффициенты (для первого и второго регрессоров; свободный член опускаю) 10.6 и -11, явно бессмысленные.
Ридж-регрессия даёт, при изменении ридж-параметра k

k      a1     a2
0.001 8.65   -8.7
...
0.01   3.22   -3.2
0.02   1.89   -1.9
0.03   1.34   -1.4
0.04   1.03   -1.1
0.05   0.843 -0.86
...
0.1     0.434 -0.45

То есть, варьируя его, получаем разумное значение (при $k=0.04$), но проблема в том, что узнать, какое из них разумное, мы не можем.
Построим главные компоненты. Собственные значения для них равны 1.99565 и .004354, то есть следует отбросить вторую, оставив лишь первую.
Нагрузки для них

0.500545171   10.7159015
0.500545171   -10.7159015

То есть мы получаем бессмысленную модель регрессии на сумму (полусумму, непринципиально) регрессоров, причём первая компонента, так полученная, имеет корреляцию с регрессандом 2.1%, меньше, чем каждая переменная в отдельности, а вторая, отброшенная, компонента, имела при жизни корреляцию 99.35%
Если же мы бы взяли лишь вторую компоненту, то коэффициенты модели были бы весьма близки к "истинным" - 1.02 и -1.02. Но идеология регрессии на главные компоненты требует отбрасывания именно её, как объясняющей лишь малую часть дисперсии регрессоров.
В предлагавшейся мною модели не была бы отброшена ни одна компонента, однако они получили бы веса (собственно, веса, равные единице, это "полностью включена", а нулю - "полностью исключена", но они таких значений в реальности не принимают, хотя могут быть близки, в данном случае первая получит вес около одной тысячной, а вторая примерно 0.999). Веса определяются не собственными значениями корреляционной матрицы, а тем, насколько компоненты скоррелированы с регрессандом (зависимой переменной).

[McConst]

Вроде бы то, о чем вы говорите - взвешивание - оно и реализовано в PLS. Во всяком случае при программировании алгоритма там тоже составляется матрица взвешенных нагрузок по числу ГК и образцов. А в деталях я пока не разбирался.

[Andrey_Kireew]

Не знаю, как Вы McConst подгоняли свою модель, но у меня на Ваших данных (по Pt), лучший результат получился $R=0,81$ на тестовой выборке при $R=0,89$ на калибровочных данных. Это была линейная модель. Но по всем признакам, она получилась весьма неплохой. Добавление произведений компонент, как в Вашей модели, привело к повышению R всего до 0,9, т.е. практически ничего не даёт.
У Вас же $R^2$ доходит до 0,99.
Уточните пожалуйста, это у Вас такой прогноз? или просто качество подгонки под калибровочные данные?
Если второе, то какое значение $R^2$ у прогноза?

Что касается ошибок, то действительно, при малых концентрациях они очень большие. Это не нормально.

[Andrey_Kireew]

Попробовал другой метод, результат совсем другой:

даже лучше, чем в представленных Вами McConst результатах. И это с первой попытки, так что результат весьма надёжный. И видимо, его можно будет ещё немного улучшить. Жаль только, что c PCA ничего толком не выходит.

[McConst]

У меня это был прогноз. Но расчитанный на коэффициентах нормирования, несколько отличающихся от исходного.


Я считал в Visual Studio C++
А в какой программе вы это делали? У меня были однажды интересные наблюдения. На работе на двух машинах у меня стоял Wolfram Mathematica одной и той же версии из одного и того же установочного пакета. Но одна машина была под Windows 7, вторая под Windows 10. Один и тот же расчетный файл давал различающиеся примерно на 5% результаты. Попробовал дома запустить эту же исходную матрицу в Unscrambler и на C++.
Unscrambler дал тоже свои результаты, но по смыслу близкие к Mathematica. Больше всего отличие дал C++. Возможно, из-за другой библиотеки матричных вычислений, а ещё я там разложение делал не через SVD, а по алгоритму NIPALS. В общем, накапливающаяся в итерациях ошибка с действиями над числами с плавающей точкой себя вполне заметно проявляет. Конечно 5% ошибка не даст такого большого различия в $r^2$
Я попробую поставить себе ваше приложение и скинуть свои данные в его формате. Можно сравнить. Если это не очередной язык программирования, типа mathlab или R. Сейчас на изучение синтаксиса нового языка катастрофически нет времени - за мной тянется хвост незавершенных задач. Могу расписать по действиям алгоритм, по которому проводил вычисления сам или скинуть код на С.

Вообще странно, если бы у меня получился плохой $r^2$, я бы наверное сразу забросил PCR как неудачный выбор, а так первичные результаты были вполне приемлемы. Плохой показатель появлялся, если не делать нормирование вообще. Правда год назад у меня было меньше образцов и каждый образец снимался по 2 раза. Сейчас в выборку добавились образцы, которые на прежней калибровке давали выбросы и вели себя странно, поэтому да, не получается так хорошо.

У вас получились очень хорошие результаты на графике, который вы мне скинули. Это график проверочных образцов. На нём есть выбросы. Наверняка есть выбросы и на калибровочных образцах. Вы можете сказать, какие имена образцов себя плохо ведут на калибровке и на проверке? Возможно следует посмотреть на их линии. Наверняка стоит повторно отдать химикам на анализ. Сталкивался с подобным пару раз, где химики просто при анализе ошиблись с коэффициентами разбавления и выдали неправильный результат. Иногда химикам мешают лантаноиды, которые есть в катализаторах.

Я уже на PCR почти махнул рукой. Выборка по объему так себе, Главных компонент и связанных с их количеством параметров много, большие компоненты увеличивают вклад в ошибку. Осталось желание в качестве проверочной выборки использовать спектры реальных аналитических образцов, обсчитанные по ММНК. Или наоборот, сделать по аналитическим спектрам с данными ММНК градуировку по PCR и посмотреть устойчивость градуировки. В градуировку, конечно войдут все ошибки ММНК, но так как на большом объеме ММНК дает ошибку около 5-7% просто из факта межлабораторных сличений для накопленной партии автокатализатора, можно понаблюдать за выбросами. И параллельно нужно дописывать код для PLS с нормированием.

Поделитесь информацией, если не сложно, что за метод вы использовали. Что-то из классики на основе латентных переменных, или это лично какие-то ваши наработки?

[Andrey_Kireew]

McConst в сообщении #1436833 писал(а):

я там разложение делал не через SVD, а по алгоритму NIPALS

Лучше всего использовать SVD. Оно даёт максимальную точность, но я обычно вычисляю через ковариационную матрицу и проблему собственных значений, т.к. в вычислительном отношении это намного проще. Мною было замечено, что все отличия в способе вычисления проявляются в основном только на слабых компонентах. Здесь их быть не должно.
Поэтому уточните, о 5% различиях чего идёт речь? о значениях коэффициентов, или чего то ещё?

McConst в сообщении #1436833 писал(а):

Могу расписать по действиям алгоритм, по которому проводил вычисления сам или скинуть код на С

Так как результаты получились лучше у Вас, наверное правильнее будет мне искать ошибку или недочёты у себя. PCA я делал в Matlab, просто через матричные вычисления. Запросто мог и ошибиться, надо бы всё хорошенько перепроверить, чем в скором времени и займусь. Вообще, как я уже писал, заметных различий быть не должно. По крайней мере, язык программирования или выбор программного пакета, не должны заметно на что то влиять. Просто, скорее всего, где то есть ошибка.

McConst в сообщении #1436833 писал(а):

Вы можете сказать, какие имена образцов себя плохо ведут на калибровке и на проверке?

Могу, это определённо Met3 и Met6. Но на метод PCA они не влияют, так как при вычислении главных компонент информация о концентрациях никак не используется.

Самое интересное, при использовании PCA большие ошибки проявляются именно на малых концентрациях, как Вы и пишите. Причина этого не совсем ясна. Нужно ещё подумать - ведь зная причину, будет легко устранить проблему, и не только в этой, но и всех будущих задачах.

McConst в сообщении #1436833 писал(а):

Поделитесь информацией, если не сложно, что за метод вы использовали

Да, метод хорош и весьма. Но смотреть в его сторону не стоит. Он всё равно нигде не реализован. Я его использовал лишь для того, чтобы посмотреть, каких вообще результатов можно добиться. Оказалось, что вполне приличных. Причём, есть потенциал для их дальнейшего улучшения. Можно будет, ради интереса, поэкспериментировать. Но у Вас чисто практическая проблема, а для неё, по обозначенной выше причине, применить всё это скорее всего всё равно не получится.

-- 25.01.2020, 16:10 --

Andrey_Kireew в сообщении #1436853 писал(а):

надо бы всё хорошенько перепроверить

перепроверил - нашел небольшую неточность. Сейчас получается: $R=0.9743$ на калибровочных данных и $R=0.9662$ на проверочных. Это лучший результат, он получается на 16 главных компонентах + константа. Т.е. оптимальное число ГК такое же как у Вас. Но качество модели у меня, всё же значительно хуже.

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

[Andrey_Kireew]

Andrey_Kireew в сообщении #1436853 писал(а):

Но качество модели у меня, всё же значительно хуже.

с этим понятно - моё уравнение не включало произведения компонент, если их добавить, то получается $R^2=0,9985$, но тогда, число независимых переменных увеличивается с 17 до 122. Лично у меня, прогноз не получился вообще. Как он получился у Вас в таких условиях, да ещё такой хороший, остаётся загадкой. Есть подозрение, что Вы всё таки используете проверочные данные при поиске главных компонент.