Função para Ajuste Potencial para Gretl

function [a_0,a_1,rmq]=num_ajuste_potencial(x,y)

    // [a_0,a_1,rmq]=num_ajuste_potencial(x,y,w)

    // Calcula os coeficientes de ajuste potencial g(x)=a_0*(x**a_1)

    // para os pontos (x_i,y_i) e o resíduo do ajuste de mínimos quadrados

    // (C) Rudnei Dias da Cunha 2014

    if and(x>0) & and(y>0)

        m = max(size(x))

        H = [ m sum(log(x)) ; sum(log(x)) sum(log(x).^2) ]

        b = [ sum(log(y)) ; sum(log(x).*log(y)) ]

        a = H\b

        a_0 = exp(a(1))

        a_1 = a(2)

        rmq = sqrt(sum((num_avalia_potencial(x,a_0,a_1)-y).^2))

    else

        a_0 = []

        a_1 = []

        rmq = []

    end

endfunction