function [y_predict] = k_hb(x_train,y_train, x_predict, h)

n = size(x_predict, 1);
l=size(h,2);

y_predict = zeros(n,l);
c = 2 * h .^ 2;
 

for i = 1:n
    % quadratische Abweichung:
    % qd = bsxfun(@minus,x_predict(i,:),x_train)).^2;
    %
    %Summe der quadr. Abweichung, wenn x_train mehrdimensional
    % s= sum(qd,1);
    %
    %Bestimme für jedes h Wi, da quadratische Abweichung der Punkte
    %für alle h gleich bleibt, sich somit nur Wi ändert
    %Wi = bsxfun(@rdivide, -s, c);
    W = zeros(n, l);
    for I=1:length(c)
        W(:,I) = exp(-sum(bsxfun(@minus,x_predict(i,:),x_train).^2,2)/c(I));
    end
    
    K = sum(W,1);
    

    for I=1:length(c)
        if (K(I)==0)
            y_predict(i,I)=0;
        else
            y_predict(i,I)=sum(W(:,I).*y_train)/K(I);
        end
    end
%     if all(K == 0)
%         y_predict(i) = 0;       
%     else
%         index=find(K==0);
%         y_predict(i,index(:))=0;
%         %wie kann ich hier sicher stellen, dass tatsächlich die operation
%         %nur an den richtigen stellen durchgefürht wird??
%         y_predict(i,:) = bsxfun(@rdivide,sum(bsxfun(@times, W, y_train),1) , K);
%     end
end