Векторизация функции содержит численное вычисление определенного интеграла в октаве

Я пытаюсь подогнать свои данные, используя формулу, содержащую численный расчет определенных интегралов с бесконечными пределами интегрирования. Для подгонки я использую октавную функцию leasqr, для которой требуется функция модели векторизации. Следующий код генерирует ошибку, возникающую при вызове численного интегрирования.

несовместимые аргументы (op1 – 1 x 387, op2 – 10 x 2)

function  [fGsAb] = GsAbs (x, p) 
Hw = 3108.0 ;  
fGsAb = Hw ./ (2.4 .* p(1) .*p(2)) .^2 .* (exp ( - (Hw - p(1) .* x) .^2 ... 
./ (2.4 .* p(1) .*p(2)) .^2 )  - exp ( - (Hw + p(1) .* x) .^2  ...
./ (2.4 .* p(1) .*p(2)) .^2 )) ;
endfunction 
function [GsDisp] = gauss_disp(x,  p) 
[GsDisp1, err] = quadgk ( @(z)  GsAbs(z, p) ./(z-x), - inf, x - 0.000001 );
[GsDisp2, err] = quadgk ( @(z)  GsAbs(z, p) ./(z-x), x + 0.000001 , inf  );
GsDisp = GsDisp1 +GsDisp2;
endfunction     
h = [200:15:6000];
pin =[1 250];
dd   =  gauss_disp (h,  pin);

Если я использую цикл:

for i = 1 : length (h)
dd (i) = gauss_disp (h (i),  pin) 
endfor

У меня нет ошибок, но я не могу использовать эту конструкцию в leasqr. Как я могу обойти это ограничение?

Заранее спасибо!


data-fitting octave numerical-integration
person Igor    schedule 25.04.2017    source источник
comment
Вопрос не ясен. Вы определяете эти функции в файле? Или в октавном терминале напрямую? И какая строка выдает вам эту ошибку?   -  person Tasos Papastylianou    schedule 25.04.2017
comment
Я дал ответ, который работает вокруг вашей проблемы. Тем не менее, в будущем, пожалуйста, предоставьте всю необходимую информацию, например, тот факт, что leasqr является функцией из пакета optim, и что это ожидает функцию определенной формы, поэтому это не удается. Мне пришлось прилично поработать детективом, чтобы догадаться об этом.   -  person Tasos Papastylianou    schedule 25.04.2017

Ответы (1)


arrow_upward
1
arrow_downward

Измените свой файл gauss_disp.m следующим образом:

function [GsDisp] = gauss_disp(x,  p) 
  [GsDisp1, err] = quadgk ( @(z)  GsAbs(z, p) ./(z-x), - inf, x - 0.000001 );
  [GsDisp2, err] = quadgk ( @(z)  GsAbs(z, p) ./(z-x), x + 0.000001 , inf  );
  GsDisp = GsDisp1 +GsDisp2;
endfunction

к этому:

function [GsDisp] = gauss_disp(x,  p)
  GsDisp  = arrayfun(@(z) gauss_disp_elementwise(z, p), x) 
endfunction 

function GsDisp = gauss_disp_elementwise(x, p)
  [GsDisp1, err] = quadgk ( @(z)  GsAbs(z, p) ./(z-x), - inf, x - 0.000001 );
  [GsDisp2, err] = quadgk ( @(z)  GsAbs(z, p) ./(z-x), x + 0.000001 , inf  );
  GsDisp = GsDisp1 +GsDisp2;
endfunction

т.е. превратите исходную функцию во вспомогательную подфункцию, которая имеет дело только со скалярными значениями, а затем используйте ее с arrayfun, чтобы получить выходные данные для всего диапазона x. Таким образом, вы можете использовать это в своей функции leasqr, например:

y = leasqr(h, dd, pin, @gauss_disp);

PS: синтаксис arrayfun предназначен только для удобства. вместо этого вы могли бы легко использовать цикл for, аналогичный тому, который вы использовали в своем вопросе.

person Tasos Papastylianou    schedule 25.04.2017
comment
прокрастинация во всей красе :-) +1 - person Andy; 26.04.2017
comment
@Энди, я только что представил диссертацию, так что мне можно, ахахах. - person Tasos Papastylianou; 26.04.2017
comment
Спасибо, такую ​​упаковку фитинговой функции я искал. Я думаю, что цикл for бесполезен для leasqr. - person Igor; 27.04.2017
comment
@Igor, я имел в виду, что arrayfun можно преобразовать в цикл for внутри функции gauss_disp. Вы правы, что это должно быть частью функции, поскольку leasqr ожидает, что эта функция будет иметь определенный формат. - person Tasos Papastylianou; 28.04.2017
comment
@Igor также, если этот ответ решит вашу проблему, пожалуйста, примите / проголосуйте соответствующим образом :) - person Tasos Papastylianou; 28.04.2017