Минимизируйте разницу между индикаторными переменными в Matlab
Я новичок в Matlab и хочу написать программу, которая выбирает значение параметра (P), чтобы минимизировать разницу между двумя векторами, где каждый вектор является переменной в фрейме данных. Первый вектор (назовем его A) представляет собой заданный вектор из единиц и нулей, а второй вектор (назовем его B) имеет каждый из своих элементов, определенных как индикаторная функция, которая зависит от значения параметра P и других переменных в фрейм данных. Например, пусть C будет третьей переменной в наборе данных, поэтому
А = [1, 0, 0, 1, 0]
В = [х, у, г, и, v]
Где x = 1, если (C[1]+10)^0,5 - P > (C[1])^0,5, иначе x = 0, и аналогично, y = 1, если (C[2]+10)^0,5 - P > (C[2])^0,5, иначе y = 0 и так далее.
Я не совсем уверен, с чего начать код, за исключением того, что может быть полезно использовать команду fminsearch. Какие-либо предложения?
Редактировать: я изменил вышеизложенное, возведя в степень, которая ближе к реальному примеру, который у меня есть. Я также предоставляю полный пример в ответ на комментарий: Пусть A такое же, как указано выше, и пусть C = [10, 1, 100, 1000, 1]. Тогда моя цель с кодом Matlab состояла бы в том, чтобы выбрать значение P, чтобы минимизировать различия между координатами векторов A и B, где B[1] = 1, если (10 + 10) ^ 0,5 - P > (10) ^ 0,5 и в противном случае B[1] = 0, и аналогично B[2] = 1, если (1 + 10) ^ 0,5 - P > (1) ^ 0,5 и в противном случае B [2] = 0 и т. д. Итак, я хочу выберите P, чтобы максимизировать вероятность того, что A[1] = B[1], A[2] = B[2] и т. д.
У меня есть следующая настройка в Matlab, где ds — это имя моего набора данных:
ds.B = zeros(size(ds,1),1); % empty vector to fill
for i = 1:size(ds,1)
if ((ds.C(i) + 10)^(0.5) - P > (ds.C(i))^(0.5))
ds.B(i) = 1;
else
ds.B(i) = 0;
end
end
Теперь я хочу выбрать значение P, чтобы минимизировать разницу между A и B. Как мне это сделать?
Обновлено: мне также интересно, как это сделать, когда неравенство похоже на (C[i]+10)^0,5 - P*D[i] > (C[i])^0,5, где D — другая переменная в мой набор данных. Теперь P — это скаляр, который умножается, а не просто складывается. Это кажется более сложным, поскольку я не могу точно определить P. Как мне решить проблему в этом случае?
Я заметил, что вы разместили связанный вопрос после этого? Какое-то время это вертелось у меня в голове, и я только что нашел ответ. Но похоже вопрос снят? Если вы хотите сделать репост, я могу дать ответ через пару дней.
Спасибо. Я разместил его здесь: stackoverflow.com/questions/75194207/…
Ответы 2
РЕДАКТИРОВАТЬ 1: кажется, что fminbnd() не является оптимальным, вероятно, из-за ступенчатого характера функции индикатора. Я обновил, чтобы проверить средние точки всех областей между переворотами функции индикатора, а также конечные точки.
Обновлено еще раз: Обновлено, чтобы включить набор данных D в качестве коэффициента P.
Если вы можете упаковать расчет расстояния в одну функцию, основанную на P, вы можете найти его минимум.
arraySize = 1000;
ds.A = double(rand([arraySize,1]) > 0.5);
ds.C = rand(size(ds.A));
ds.D = rand(size(ds.A));
B = @(P)double((ds.C+10).^0.5 - P.*ds.D > ds.C.^0.5);
costFcn = @(P)sqrt(sum((ds.A-B(P)).^2));
% Solving the equation (C+10)^0.5 - P*D = C^0.5 for P, and sorting the results
BCrossingPoints = sort(((ds.C+10).^0.5-ds.C.^0.5)./ds.D);
% Taking the average of each crossing point with its neighbors
BMidpoints = (BCrossingPoints(1:end-1)+BCrossingPoints(2:end))/2;
% Appending endpoints onto the midpoints
PsToTest = [BCrossingPoints(1)-0.1; BMidpoints; BCrossingPoints(end)+0.1];
% Calculate the distance from A to B at each P to test
costResult = arrayfun(costFcn,PsToTest);
% Find the minimum cost
[~,lowestCostIndex] = min(costResult);
% Find the optimum P
optimumP = PsToTest(lowestCostIndex);
ds.B = B(optimumP);
semilogx(PsToTest,costResult)
xlabel('P')
ylabel('Distance from A to B')
Большое спасибо за ваш полезный ответ. У меня есть один уточняющий вопрос. Кажется, что даже когда я немного меняю функцию, например. изменяя мощность или добавляемую константу (10), я всегда получаю вектор ds.B, который либо содержит все 0, либо все единицы. Является ли это результатом того, как устроен код, или это характер конкретной функции, которую я выбрал? В идеале я хотел бы разрешить вариант, что элементы B могут отличаться.
Вероятно, это результат работы входных данных и функции индикатора. Мне удалось получить B = [0 0 0 1 0] с C = [50 0 100 500 3] и функцией C+10-P>C.^0.5
Что ж, похоже, что fminbnd может не подойти для вашего приложения, вероятно, из-за ступенчатого характера функции расстояния. Я изменил свой ответ, чтобы использовать другой метод нахождения минимального расстояния.
Спасибо за редактирование. Когда я применяю это к своим данным, я получаю следующую ошибку, связанную со строкой PsToTest = [BCrossingPoints(1)-0.1, BMidpoints, BCrossingPoints(end)+0.1]; Error using horzcat Dimensions of arrays being concatenated are not consistent. Не могли бы вы объяснить, как определяются BCrossingPoints и BMidpoints? Я еще не следую интуиции.
@ RK124 Глядя на ваш пример выше, кажется, что вы используете векторы-столбцы, в то время как в моем примере использовались векторы-строки. Я отредактировал свой ответ, чтобы он соответствовал вашему, что означало замену запятых в массивах точками с запятой. Теперь он должен работать с вашими данными. Я также добавил комментарии к ответу, чтобы объяснить, для чего предназначен каждый шаг.
Еще раз спасибо за вашу помощь в этом. Мне интересно, могу ли я использовать аналогичный метод для решения значения, которое умножается на вектор (т.е. на одну из других переменных). Например, если уравнение имеет вид B = @(P)double((ds.C+10).^0,5 - P*ds.D > ds.C.^0,5), где ds.D — другая переменная в наборе данных . Когда я пытаюсь разделить все это на ds.D, чтобы определить BCrossingPoints (т.е. найти P), я получаю ошибку Requested 37885x37885 (21.4GB) array exceeds maximum array size preference (8.0GB). Есть ли у вас какие-либо предложения?
Это должно работать нормально. Каковы размеры переменных непосредственно перед шагом, который вызывает ошибку? Похоже, что один из ваших наборов данных представляет собой вектор-строку, а другой — вектор-столбец, поэтому, когда вы пытаетесь выполнить поэлементные операции, они неявно расширяются в квадрат. Для получения дополнительной информации попробуйте здесь
Немного подумав об этом, я не уверен, что концептуально понимаю, как адаптировать вышеизложенное к случаю, когда B = @(P)double((ds.C+10).^0,5 - P*ds.D > ds.C.^0,5). Поскольку «деление вектора на другой вектор» на самом деле не определяется математически, я действительно не понимаю, как определить BCrossingPoints в этом случае ... Я отредактирую свой исходный вопрос, чтобы узнать об этом немного глубже.
Я так понимаю, вы хотите, чтобы каждый элемент B зависел только от соответствующих ему партнеров в C и D? Вы захотите использовать поэлементное деление, а не матричное деление. Я отредактирую свой ответ выше, чтобы приспособиться. Вы можете узнать больше здесь.
1.- x предполагается только положительным действительным, потому что с x<0 появляются комплексные значения.
Поскольку в вопросе нет комментариев, кажется разумным предположить только x real и x>0.
В соответствии с запросом P 'параметр' является скаляром, P имеет только 2 значимых состояния >0 или <0, давайте посмотрим, как это:
2.- Следующие строки генерируют случайные числа A и C.
Затем выполняется развертка p и вычисляются расстояния d1 и d2.
D1 — евклидово расстояние, а d2 — абсолютная разница между A и B при преобразовании из двоичного в десятичное:
N=10
% A=[1 0 0 1 0]
A=randi([0 1],1,N);
% C=[10 1 1e2 1e3 1]
C=randi([0 1e3],1,N)
p=[-1e4:1:1e4]; % parameter to optimize
B=zeros(1,numel(A));
d1=zeros(1,numel(p)); % euclidean distance
d2=zeros(1,numel(p)); % difference distance
for k1=1:1:numel(p)
B=(C+10).^.5-p(k1)>C.^.5;
d1(k1)=(sum((B-A).^2))^.5;
d2(k1)=abs(sum(A.*2.^[numel(A)-1:-1:0])-sum(B.*2.^[numel(A)-1:-1:0]));
end
figure;
plot(p,d1)
grid on
xlabel('p');title('d1')
figure
plot(p,d2)
grid on
xlabel('p');title('d2')
Единственная степень свободы для оптимизации, по-видимому, — это знак P независимо от |P| ценить.
3.- f(p,x) либо не имеет корня, либо имеет только один корень, в зависимости от p
Пороговая функция
if f(x)>0 then B(k)==1 else B(k)==0
Это
f(p,x)=(x+10)^.5-p-x^.5
Сейчас
(x+10).^.5-p>x.^.5 такое же, как (x+10).^.5-x.^.5>p
Существует диапазон p, который сохраняет f(p,x)=0 без какого-либо (настоящего) корня.
Для частного случая p=0 тогда (x+10).^.5 и x.^.5 не пересекаются (пока Inf не достигнуто = пересечения нет)
figure;plot(x,(x+10).^.5,x,x.^.5);grid on
[![enter image description here][3]][3]
y2=diff((x+10).^.5-x.^.5)
figure;plot(x(2:end),y2);
grid on;xlabel('x')
title('y2=diff((x+10).^.5-x.^.5)')
[![enter image description here][3]][3]
% 005
Это означает, что условие f(x)>0 всегда истинно, если все биты B=1. При B=1 d(A,B) превращается в d(A,1), константу.
Однако для определенного значения p есть один корень, а f(x)>0 всегда ложно, сохраняя все биты B=0.
В этом случае d(A,B) функция стоимости превращается в d(A,0), а это и есть сама A.
4.- P как вектор
Оптимизация выигрывает в степенях свободы, если вместо скаляра P рассматривать P как вектор.
Для данного x существует значение p, которое переключает B(k) с 0 на 1.
Любое значение p ниже этого порога сохраняет B(k)=0.
Эквивалентно, инвертирование f(x) :
g(p)=(10-p^2)^2/(4*p^2)>x
Значения x ниже этого порога приближают B к A, потому что для каждого элемента B оно переворачивается на значение элемента A.
Поэтому удобно рассматривать P как вектор, а не как скаляр, и:
Чтобы все или как можно больше (насколько это возможно) элементов C соответствовали c(k)<(10-p^2)^2/(4*p^2), чтобы получить C=A или минимизировать d(A,C)
5.- корни f(p,x)
syms t positive
p=[-1000:.1:1000];
zp=NaN*ones(1,numel(p));
sol=zeros(1,numel(p));
for k1=1:1:numel(p)
p(k1)
eq1=(t+10)^.5-p(k1)-t^.5-p(k1)==0;
s1=solve(eq1,t);
if ~isempty(s1)
zp(k1)=s1;
end
end
nzp=~isnan(zp);
zp(nzp)
Возвращается
=
620.0100 151.2900 64.5344 34.2225 20.2500 12.7211
8.2451 5.4056 3.5260 2.2500 1.3753 0.7803
0.3882 0.1488 0.0278
Ваш пункт № 2 страдает от отсутствия разрешения выборки. Если вы измените на p=[0:.01:2];, вы увидите интересное поведение.
Я не уверен, что полностью понял пункт № 3. При p==0 я вижу, что B всегда будет ложным. Однако при положительном P наличие корня означает, что при малых значениях C результат функции будет слева от этого корня и, следовательно, истинен, а при больших значениях C результат функции будет справа от этого корня и, следовательно, ложен. Затем, учитывая массив независимых значений C, вы получите массив независимых значений B.
Предлагаемое уменьшение шага для увеличения разрешения остается за автором вопроса. Иногда может потребоваться улучшение данного ответа, но решение указывает на правильное направление. Почему бы нам не оставить автору вопроса, если он заинтересован, добавить дополнительные комментарии?
Напротив, тот факт, что существует область P, где функция стоимости возвращает расстояния, отличные от экстремумов, делает ваш вывод в пункте 2 просто неверным. На самом деле не имеет значения, замечаю ли это я или ОП. Кроме того, это на самом деле ведет в неправильном направлении, предполагая, что можно просто выбрать положительное или отрицательное P и сделать это, когда фактический минимум (который искал ОП) обычно находится в этой небольшой средней области.
Итак, если P<10 -> x = 1, иначе x = 0, независимо от значения C[1], правильно ли это? Не могли бы вы привести полный пример ожидаемого результата, полностью написав значения A, C и P, а также результирующее значение B?