2015-04-30 4 views
0

Я такие векторы и х у векторов, описанные в поле ниже:Matlab, чтобы найти соответствующие коэффициенты

x = [0 5 8 15 18 25 30 38 42 45 50]; 
y = [81.94 75.94 70.06 60.94 57.00 50.83 46.83 42.83 40.94 39.00 38.06]; 

с этими значениями, как я могу найти коэффициенты y = a*(b^x) ??

Я попробовал этот код, но он находит для y = a*e^(b*x)

clear, clc, close all, format compact, format long 
% enter data 
x = [0 5 8 15 18 25 30 38]; 
y = [81.94 75.94 70.06 60.94 57.00 50.83 46.83 42.83]; 
n = length(x); 
y2 = log(y); 
j = sum(x); 
k = sum(y2); 
l = sum(x.^2); 
m = sum(y2.^2); 
r2 = sum(x .* y2); 
b = (n * r2 - k * j)/(n * l - j^2) 
a = exp((k-b*j)/n) 

y = a * exp(b * 35) 
result_68 = log(68/a)/b 

Я знаю методы интерполяции, но я не мог реализовать его к моим решениям ...

Большое спасибо заранее!

ответ

4

С у = * б^х такая же, как бревно (у) = журнал () + х лог (б) вы можете сделать

>> y = y(:); 
>> x = x(:); 
>> logy = log(y); 
>> beta = regress(logy, [ones(size(x)), x]); 
>> loga = beta(1); 
>> logb = beta(2); 
>> a = exp(loga); 
>> b = exp(logb); 

поэтому значения в и б являются

>> a, b 
a = 
     78.8627588780382 
b = 
    0.984328823937827 

откладывая подходят

>> plot(x, y, '.', x, a * b .^ x, '-') 

дает вам это -

enter image description here

NB функция regress является из панели инструментов статистики, но вы можете определить очень простая версия, которая делает то, что вам нужно

function beta = regress(y, x) 
    beta = (x' * x) \ (x' * y); 
end 
+0

Мне нравится, что в то время как это на самом деле математический вопрос, более чем вопрос программирования, она была решена очень быстро. – Henrik

+0

Я бы сказал, что я математик, больше, чем программист, поэтому этот вопрос явно хорошо подходит для меня! –

+0

Почему я не могу реализовать, как вы упомянули, например log (y) = log (a) + x log (b)? Без регресса или бета-версии, возможно, я занимаюсь этим кодом, но на самом деле я не получил решение. Не могли бы вы описать это больше? – Hayra

1

В качестве дополнения к полученному ответу Криса Тейлор, который обеспечивает наилучшую линейную подгонку в логарифмической трансформированных области, вы можете найти лучшее прилегание в исходной области пути решения непосредственно нелинейная проблема , например, Gauss-Newton algorithm

Используя, например, решение, данное Крис в качестве отправной точки:

x = [0 5 8 15 18 25 30 38 42 45 50]; 
y = [81.94 75.94 70.06 60.94 57.00 50.83 46.83 42.83 40.94 39.00 38.06]; 

regress = @(y, x) (x' * x) \ (x' * y); 

y = y(:); 
x = x(:); 
logy = log(y); 
beta = regress(logy, [ones(size(x)), x]); 
loga = beta(1); 
logb = beta(2); 
a = exp(loga) 
b = exp(logb) 
error = sum((a*b.^x - y).^2) 

Что дает:

>> a, b, error 
a = 
    78.862758878038164 
b = 
    0.984328823937827 
error = 
    42.275290442577422 

Вы можете перебирать немного дальше, чтобы найти лучшее решение

beta = [a; b]; 
iter = 20 
for k = 1:iter 
    fi = beta(1)*beta(2).^x; 
    ri = y - fi; 
    J = [ beta(2).^x, beta(1)*beta(2).^(x-1).*x ]'; 
    JJ = J * J'; 
    Jr = J * ri; 
    delta = JJ \ Jr; 
    beta = beta + delta; 
end 

a = beta(1) 
b = beta(2) 
error = sum((a*b.^x - y).^2) 

Давать:

>> a, b, error 

a = 
    80.332725222265623 
b = 
    0.983480686478288 
error = 
    35.978195088265906 
+0

Спасибо. Похоже, ваша петля делает 20 итераций градиентного спуска - это правильно? –

+0

В данном случае это Гаусс-Ньютон. Он быстрее конвергенции и не требует гессиана. Мы могли бы реализовать градиентный спуск или даже чистый Ньютон, так как hessian действительно легко вычислить для этой функции. – ibancg

0

Еще один вариант заключается в использовании MATLAB-х Curve Fitting Toolbox, который позволяет создавать подгонку в интерактивном режиме, а также может испускать код MATLAB, необходимый для запуска подгонки.

Curve fitting tool screenshot

Смежные вопросы