0

我的主脚本包含以下代码:

%# Grid and model parameters
nModel=50;
nModel_want=1;
nI_grid1=5;
Nth=1;
nRow.Scale1=5;
nCol.Scale1=5;
nRow.Scale2=5^2;
nCol.Scale2=5^2;

theta = 90; % degrees
a_minor = 2; % range along minor direction
a_major = 5; % range along major direction
sill = var(reshape(Deff_matrix_NthModel,nCell.Scale1,1)); % variance of the coarse data matrix of size nRow.Scale1 X nCol.Scale1

%#  Covariance computation

% Scale 1
for ihRow = 1:nRow.Scale1
    for ihCol = 1:nCol.Scale1
        [cov.Scale1(ihRow,ihCol),heff.Scale1(ihRow,ihCol)] = general_CovModel(theta, ihCol, ihRow, a_minor, a_major, sill, 'Exp');
    end
end

% Scale 2
for ihRow = 1:nRow.Scale2
    for ihCol = 1:nCol.Scale2
        [cov.Scale2(ihRow,ihCol),heff.Scale2(ihRow,ihCol)] = general_CovModel(theta, ihCol/(nCol.Scale2/nCol.Scale1), ihRow/(nRow.Scale2/nRow.Scale1), a_minor, a_major, sill/(nRow.Scale2*nCol.Scale2), 'Exp');
    end
end


%# Scale-up of fine scale values by averaging
[covAvg.Scale2,var_covAvg.Scale2,varNorm_covAvg.Scale2] = general_AverageProperty(nRow.Scale2/nRow.Scale1,nCol.Scale2/nCol.Scale1,1,nRow.Scale1,nCol.Scale1,1,cov.Scale2,1);

我在我的主脚本中使用了两个函数,general_CovModel()并且general_AverageProperty(),它们给出如下:

function [cov,h_eff] = general_CovModel(theta, hx, hy, a_minor, a_major, sill, mod_type)
% mod_type should be in strings

angle_rad = theta*(pi/180); % theta in degrees, angle_rad in radians
R_theta = [sin(angle_rad) cos(angle_rad); -cos(angle_rad) sin(angle_rad)];
h = [hx; hy];
lambda = a_minor/a_major;
D_lambda = [lambda 0; 0 1];
h_2prime = D_lambda*R_theta*h;
h_eff = sqrt((h_2prime(1)^2)+(h_2prime(2)^2));

if strcmp(mod_type,'Sph')==1 || strcmp(mod_type,'sph') ==1
    if h_eff<=a
        cov = sill - sill.*(1.5*(h_eff/a_minor)-0.5*((h_eff/a_minor)^3));
    else
        cov = sill;
    end
elseif strcmp(mod_type,'Exp')==1 || strcmp(mod_type,'exp') ==1
    cov = sill-(sill.*(1-exp(-(3*h_eff)/a_minor)));
elseif strcmp(mod_type,'Gauss')==1 || strcmp(mod_type,'gauss') ==1
    cov = sill-(sill.*(1-exp(-((3*h_eff)^2/(a_minor^2)))));
end

function [PropertyAvg,variance_PropertyAvg,NormVariance_PropertyAvg]=...
    general_AverageProperty(blocksize_row,blocksize_col,blocksize_t,...
    nUpscaledRow,nUpscaledCol,nUpscaledT,PropertyArray,omega)
% This function computes average of a property and variance of that averaged
% property using power averaging

PropertyAvg=zeros(nUpscaledRow,nUpscaledCol,nUpscaledT);

%# Average of property
for k=1:nUpscaledT,
    for j=1:nUpscaledCol,
        for i=1:nUpscaledRow,
            sum=0;
            for a=1:blocksize_row,
                for b=1:blocksize_col,
                    for c=1:blocksize_t,
                        sum=sum+(PropertyArray((i-1)*blocksize_row+a,(j-1)*blocksize_col+b,(k-1)*blocksize_t+c).^omega); % add all the property values in 'blocksize_x','blocksize_y','blocksize_t' to one variable
                    end
                end
            end
            PropertyAvg(i,j,k)=(sum/(blocksize_row*blocksize_col*blocksize_t)).^(1/omega); % take average of the summed property
        end
    end
end

%# Variance of averageed property
variance_PropertyAvg=var(reshape(PropertyAvg,...
    nUpscaledRow*nUpscaledCol*nUpscaledT,1),1,1); 

%# Normalized variance of averageed property
NormVariance_PropertyAvg=variance_PropertyAvg./(var(reshape(...
    PropertyArray,numel(PropertyArray),1),1,1)); 

问题:使用 Matlab,我想通过扰动/改变以下任何(或所有)变量来优化covAvg.Scale2它以使其与cov.Scale1

1)a_minor

2)a_major

3)theta

我知道我可以使用fminsearch,但是,我无法在使用它时扰乱我想要的变量fminsearch

4

1 回答 1

2

我不会假装理解你所做的一切。但这听起来像是一个典型的最小化问题。你要做的是想出一个函数,它接受a_minora_major作为参数,并返回和theta之间的差的平方。像这样的东西:covAvg.Scale2cov.Scale1

function diff = minimize_me(a_minor, a_major, theta)
    ... your script goes here
    diff = (covAvg.Scale2 - cov.Scale1)^2;
end

然后你需要matlab来最小化这个函数。这里有不止一种选择。由于您只有三个变量要最小化,fminsearch因此是一个很好的起点。你可以这样称呼它:

opts = optimset('display', 'iter');
x = fminsearch( @(x) minimize_me(x(1), x(2), x(3)), [a_minor_start a_major_start theta_start], opts)

第一个参数fminsearch是要优化的函数。它必须采用单个参数:将被扰动以找到最小值的变量向量。在这里,我使用一个匿名函数从这个向量中提取值并将它们传递给minimize_me. to的第二个参数fminsearch是一个包含开始搜索的值的向量。第三个参数是影响搜索的选项;display当您第一次开始优化时设置为一个好主意iter,这样您就可以了解优化器正在收敛。

如果您的参数具有受限域(例如,它们都必须是正数),请查看fminsearchbnd文件交换。

如果我误解了您的问题,而这根本没有帮助,请尝试发布我们可以运行的代码来自己重现问题。

于 2012-11-15T22:11:27.920 回答