Блог

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

Чтобы уточнить: я хочу позвонить logPosterior с новым значением theta , но также с выводом theta и logpdf из предыдущего вызова. И мне нужно сделать это с помощью дескриптора функции, который будет вызываться несколько раз функцией, которую я не контролирую, поэтому мне нужен либо logPosterior один, либо его дескриптор для управления хранением старых значений. В первом вызове должны быть разные значения theta и old_theta , чтобы функция могла работать.

 %% Toy implementation of hmcsampler class in Matlab
NumPredictors = 2;

trueIntercept = 2;
trueBeta = [3;0];
NumData = 100;
rng('default') %For reproducibility
X = rand(NumData,NumPredictors);
mu = X*trueBeta   trueIntercept;
y = mu;

% define the mean and variance of normal distribution of each parameter
means = [0; 0];
standevs = [1;1];


% create multivariate normal log probability distribution
[logpdf, grad_logpdf] = @(theta)logPosterior(theta, old_theta, X, y, means, standevs, old_logpdf); % <- How to write this?
% create the startpoint from which sampling starts
startpoint = randn(2, 1);

% create an HMC sampler object
smp = hmcSampler(logpdf, startpoint);

% estimate maximum of log probability density
[xhat, fitinfo] = estimateMAP(smp);

num_chains = 4;
chains = cell(num_chains, 1);
burnin = 50000;
num_samples = 2000000;


function [logpdf, grad_logpdf] = logPosterior(theta, old_theta, X, y, means, standevs, old_logpdf)
    
    % values
    intercept = theta(1);
    beta = theta(2:end);
    y_computed = X*beta   intercept; 
    log_likelihood = log(y_computed);
    del_loglikelihood = log_likelihood - old_logpdf;
    del_params = theta - old_theta;
    grad_params1 = del_loglikelihood/del_params;
    
    % compute log priors and gradients of parameters
    log_prior_params = 0;
    grad_params2 = [];
    for i = 1:3
        [lp, grad] = normalDistGrad(theta(i), means(i), standevs(i));
        log_prior_params = log_prior_params   lp;
        grad_params2 = [grad_params2; grad];
    end
    % return the log posterior and its gradient
    logpdf = log_likelihood   log_prior_params;
    grad_logpdf = grad_params1   grad_params2;
end

function [lpdf,glpdf] = normalDistGrad(X, Mu, Sigma)
    Z = (X - Mu)./Sigma;
    lpdf = sum(-log(Sigma) - .5*log(2*pi) - .5*(Z.^2));
    glpdf = -Z./Sigma;
end
 

Я бы реализовал logPosterior следующее, чтобы он отслеживал значения в последнем вызове функции. persistent Переменная является локальной для функции, но сохраняется между вызовами, что делает ее идеальным инструментом для предоставления памяти функции.

 function [logpdf, grad_logpdf] = logPosterior(theta, X, y, means, standevs)
    persistent old_theta old_logpdf
    if isempty(old_theta)
       % function hasn't been called before, initialize the old values:
       old_theta = zeros(size(theta));
       old_logpdf = 0; % adjust to be meaningful initial values
    end

    % ... (your original function code here)

    % save new values as old values for next call
    old_theta = theta;
    old_logpdf = logpdf;
end
 

Теперь вы должны взять дескриптор функции следующим образом:

 func = @(theta)logPosterior(theta, X, y, means, standevs);
 

func это дескриптор, который вы передадите в любую функцию, которая его вызовет. Вы могли бы инициализировать предыдущие переменные, запустив свою функцию один раз (я не уверен, что это подходящая инициализация!):

 func(startpoint);
smp = hmcSampler(func, startpoint);
 

справка matlab

сначала определите свою функцию в отдельном файле m, который находится в вашем текущем пути :

 function [o1,o2]=myfunc(in1,in2,in3)
    o1=in1 in2 in3;
    o2=in3-in2;
end
 

затем создайте дескриптор функции :

 f=@myfunc;
f(1,2,3)
ans= 6
 

чтобы использовать только один из входов:

 f=@(x)myfunc(x,3,5);
f(1)
ans=9
 

чтобы получить оба вывода:

 [a,b]=f(1);
 

или определите один вывод в своей функции и ссылайтесь на их индекс после их вызова :

 function o=myfunc(in1,in2,in3)
    o1=in1 in2 in3;
    o2=in3-in2;
    o=[o1,o2];
end
 

…