0

q実数のベクトルを入力として受け取り、出力として複素行列を返すアルゴリズムを実装するためのコードをいくつか作成しましたR。q以下のMatlabコードは、入力ベクトルと出力行列を示すプロットを生成しますR。

複素行列出力のみが与えられた場合R、入力ベクトルを取得したいと思いますq。最小二乗最適化を使用してこれを行うことはできますか？rs_rコード（および）には再帰的な実行合計があるrs_iため、出力行列の列の計算は、前の列の計算に依存します。

おそらく、非線形最適化を設定してq、出力行列から入力ベクトルを再構成することができますRか？

これを別の方法で見ると、アルゴリズムを使用して行列を計算しましたR。qアルゴリズムを「逆に」実行して、出力行列から入力ベクトルを取得したいと思いますR。

出力から開始値を再構成して問題を「ブラックボックス」として扱う方法がない場合は、モデル自体の数学を最適化に使用できる可能性がありますか？プログラムは次の方程式を評価します。

方程式

Utilde（tau、omega）は出力マトリックスRです。タウ（時間）変数は応答行列の列を構成しR、オメガ（頻度）変数は応答行列の行を構成しますR。積分は、tau=0から現在のtauタイムステップまでの再帰的な実行合計として実行されます。

以下に掲載されているプログラムによって作成されたプロットは次のとおりです。

q値入力 行列出力

完全なプログラムコードは次のとおりです。

N = 1001;
q = zeros(N, 1); % here is the input
q(1:200) = 55;
q(201:300) = 120;
q(301:400) = 70;
q(401:600) = 40;
q(601:800) = 100;
q(801:1001) = 70;
dt = 0.0042;
fs = 1 / dt;
wSize = 101;
Glim = 20;
ginv = 0;
R = get_response(N, q, dt, wSize, Glim, ginv); % R is output matrix
rows = wSize; 
cols = N;

figure; plot(q); title('q value input as vector'); 
ylim([0 200]); xlim([0 1001])

figure; imagesc(abs(R)); title('Matrix output of algorithm')
colorbar

計算を実行する関数は次のとおりです。

function response = get_response(N, Q, dt, wSize, Glim, ginv)

fs = 1 / dt; 
Npad = wSize - 1; 
N1 = wSize + Npad;
N2 = floor(N1 / 2 + 1);
f = (fs/2)*linspace(0,1,N2);
omega = 2 * pi .* f';
omegah = 2 * pi * f(end);
sigma2 = exp(-(0.23*Glim + 1.63));

sign = 1;
if(ginv == 1)
    sign = -1;
end

ratio = omega ./ omegah;
rs_r = zeros(N2, 1);  
rs_i = zeros(N2, 1);   
termr = zeros(N2, 1);
termi = zeros(N2, 1);
termr_sub1 = zeros(N2, 1);
termi_sub1 = zeros(N2, 1);
response = zeros(N2, N);

 % cycle over cols of matrix
for ti = 1:N               

    term0 = omega ./ (2 .* Q(ti));
    gamma = 1 / (pi * Q(ti));

    % calculate for the real part
    if(ti == 1)
        Lambda = ones(N2, 1);
        termr_sub1(1) = 0;  
        termr_sub1(2:end) = term0(2:end) .* (ratio(2:end).^-gamma);  
    else
        termr(1) = 0; 
        termr(2:end) = term0(2:end) .* (ratio(2:end).^-gamma); 
        rs_r = rs_r - dt.*(termr + termr_sub1);
        termr_sub1 = termr;
        Beta = exp( -1 .* -0.5 .* rs_r );

        Lambda = (Beta + sigma2) ./ (Beta.^2 + sigma2);  % vector
    end 

    % calculate for the complex part  
    if(ginv == 1)  
        termi(1) = 0;
        termi(2:end) = (ratio(2:end).^(sign .* gamma) - 1) .* omega(2:end);
    else
        termi = (ratio.^(sign .* gamma) - 1) .* omega;
    end
    rs_i = rs_i - dt.*(termi + termi_sub1);
    termi_sub1 = termi;
    integrand = exp( 1i .* -0.5 .* rs_i );

    if(ginv == 1) 
        response(:,ti) = Lambda .* integrand;
    else        
        response(:,ti) = (1 ./ Lambda) .* integrand;
    end  
end % ti loop
4

2 に答える 2

3

いいえ、このプロセスの「モデル」自体を知らない限り、そうすることはできません。プロセスを完全なブラックボックスとして扱う場合、一般的には不可能ですが、特定の場合には何かが起こる可能性があります。

基礎となるプロセスを知っていても、最小二乗推定量は開始値に依存しているため、機能しない可能性があります。そのため、そこに適切な推測がない場合は、間違ったパラメーターのセットに収束する可能性があります。

于 2012-08-22T23:43:47.200 に答える
0

モデルの数学を使用することにより、入力を推定できることがわかります。これは一般的には当てはまりませんが、私の問題ではうまくいくようです。累積積分は偏導関数によって除去されます。

N = 1001;
q = zeros(N, 1);
q(1:200) = 55;
q(201:300) = 120;
q(301:400) = 70;
q(401:600) = 40;
q(601:800) = 100;
q(801:1001) = 70;
dt = 0.0042;
fs = 1 / dt;
wSize = 101;
Glim = 20;
ginv = 0;
R = get_response(N, q, dt, wSize, Glim, ginv);
rows = wSize; 
cols = N;
cut_val = 200;

imagLogR = imag(log(R));

Mderiv = zeros(rows, cols-2);
for k = 1:rows
   val = deriv_3pt(imagLogR(k,:), dt);
   val(val > cut_val) = 0;
   Mderiv(k,:) = val(1:end-1);
end

disp('Running iteration');
q0 = 10;
q1 = 500;
NN = cols - 2;
qout = zeros(NN, 1);
for k = 1:NN
    data = Mderiv(:,k); 
    qout(k) = fminbnd(@(q) curve_fit_to_get_q(q, dt, rows, data),q0,q1);
end

figure; plot(q); title('q value input as vector'); 
ylim([0 200]); xlim([0 1001])

figure;
plot(qout); title('Reconstructed q')
ylim([0 200]); xlim([0 1001])

サポート機能は次のとおりです。

function output = deriv_3pt(x, dt)

% Function to compute dx/dt using the 3pt symmetrical rule
% dt is the timestep

N = length(x);
N0 = N - 1;
output = zeros(N0, 1);
denom = 2 * dt;

for k = 2:N0 
   output(k - 1) = (x(k+1) - x(k-1)) / denom;  
end


function sse = curve_fit_to_get_q(q, dt, rows, data)

fs = 1 / dt;
N2 = rows;
f = (fs/2)*linspace(0,1,N2);  % vector for frequency along cols
omega = 2 * pi .* f';
omegah = 2 * pi * f(end);
ratio = omega ./ omegah;

gamma = 1 / (pi * q);

termi = ((ratio.^(gamma)) - 1) .* omega;

Error_Vector =  termi - data;
sse = sum(Error_Vector.^2);

オリジナル 推定

于 2012-08-24T02:39:39.890 に答える