clear; // *** 二次元ガウス分布 *** mu1 = 3; mu2 = 1; mu = [mu1; mu2]; sigma1 = sqrt(10); sigma2 = sqrt(10); sigma12 = 5; SIGMA = [sigma1^2, sigma12; sigma12, sigma2^2]; rho = sigma12 / (sigma1 * sigma2); function z = func(x, y) z = -1 ./ (2 * %pi * sigma1 * sigma2 * sqrt(1 - rho)) .. .* exp(-1 / (2 .* (1 - rho^2)) .. .* ((x - mu1) .^ 2 ./ (sigma1 ^ 2) - 2 .* rho .* (x - mu1) .* (y - mu2) ./ (sigma1 * sigma2) + ((y - mu2) .^ 2) ./ (sigma2 ^ 2))) endfunction // *** 数値微分 *** h = 1E-3; function dzx = dfx(x, y) dzx = (func((x+h), y) - func((x-h), y)) ./ (2 * h) endfunction function dzy = dfy(x, y) dzy = (func(x, (y+h)) - func(x, (y-h))) ./ (2 * h) endfunction // *** グラフのプロット *** x = linspace(-10,10); y = linspace(-10,10); [X,Y] = ndgrid(x,y); Z = func(X,Y); // 色の設定 //xset("colormap",jetcolormap(64)) // xset("colormap",graycolormap(64)) // Sgrayplot(x,y,Z) xset("fpf"," "); // 等高線に値を表示しない contour2d(x,y,Z,10); xmin = min(X); xmax = max(X); ymin = min(Y); ymax = max(Y); zoom_rect([xmin,ymin,xmax,ymax]); // *** 最小値の計算 *** // 停止条件 err = 1E-5; a = 200; // 初期値 x = -4; y = 4; z = func(x, y); dx = dfx(x, y); dy = dfy(x, y); plot(x, y, "xk"); // 最小値の計算 while ((abs(dx) > err) | (abs(dy) > err)) x = x - a * dx; y = y - a * dy; dx = dfx(x, y); dy = dfy(x, y); plot(x, y, ".r"); end // 計算結果 x, y plot(x, y, "xk");