| Hogfather |
05.03.2015 09:53 |
Вложений: 1
Цитата:
Сообщение от Realwert
(Сообщение 511672)
Хог, благодарю, а скиньте сам файл, может будет время разберусь в проге.
|
(восхищенно) А скопировать в редактор и стереть лишнее религия не позволяет?
Ладно уж, пользуйтесь моей добротой.
Ниже написанное надо скопировать и вставить в окно R или RStudio.
Код:
Re<-rep(c(26575.0,22145.8,17716.6,14763.9,10334.7,7381.9),4)
T<-c(rep(1,12),rep(2,12))
Sd<-rep(c(rep(1,6),rep(2,6)),2)
Nu<-c(142.3,121.8,114.8,105.8,88.4,74.4,
139.0,125.4,110,102.1,81,73.2,
140.0,120.0,114.7,104.0,87.0,72.0,
137.0,123.3,109.0,101.0,80.1,71.0)
MyData<-data.frame(Re=Re,T=T,Sd=Sd,Nu=Nu)
MyData
summary(lm1<-lm(log(Nu)~log(Re)+log(T)+log(Sd)))
confint(lm1)
summary(lm2<-lm(log(Nu)~log(Re)+log(Sd)))
confint(lm2)
oldpar<-par(mfrow=c(2,2))
plot(lm2)
par(oldpar)
В приведенных выше скриптах все команды начинаются с ">" (символ приглашения консоли) или "+" (символ приглашения консоли, если оператор продолжается на несколько строк). Остальное все -- результаты расчетов. Поэтому можно просто копировать в любой редактор, чистить от лишнего и пользоваться.
Кстати, интересно, что если использовать информационный критерий Акаике для уменьшения числа предикторов, то оставлять надо все .
Код:
> step(lm1)
Start: AIC=-174.57
log(Nu) ~ log(Re) + log(T) + log(Sd)
Df Sum of Sq RSS AIC
<none> 0.01193 -174.568
- log(T) 1 0.00152 0.01345 -173.689
- log(Sd) 1 0.00504 0.01696 -168.112
- log(Re) 1 1.16787 1.17980 -66.305
Call:
lm(formula = log(Nu) ~ log(Re) + log(T) + log(Sd))
Coefficients:
(Intercept) log(Re) log(T) log(Sd)
-0.19233 0.50427 -0.02296 -0.04180
Нетрудно увидеть, что модель тогда почти совпадает с моделью коллеги avz. Есть расхождение в множителе А.
Код:
> lm1
Call:
lm(formula = log(Nu) ~ log(Re) + log(T) + log(Sd))
Coefficients:
(Intercept) log(Re) log(T) log(Sd)
-0.19233 0.50427 -0.02296 -0.04180
> exp(-0.19233)
[1] 0.8250346
Проверка модели в Gretl (данные в формате Gretl во вложении. В меню Gretl выбрать Модель->МНК и заполнить форму)
Код:
Модель 1: МНК, использованы наблюдения 1-24
Зависимая переменная: l_Nu
Коэффициент Ст. ошибка t-статистика P-значение
---------------------------------------------------------------
const -0,192327 0,109970 -1,749 0,0956 *
l_Re 0,504269 0,0113950 44,25 1,97e-021 ***
l_T -0,0229633 0,0143830 -1,597 0,1260
l_Sd -0,0418047 0,0143830 -2,907 0,0087 ***
Среднее зав. перемен 4,636746 Ст. откл. зав. перемен 0,227114
Сумма кв. остатков 0,011927 Ст. ошибка модели 0,024420
R-квадрат 0,989946 Испр. R-квадрат 0,988438
F(3, 20) 656,4516 Р-значение (F) 3,89e-20
Лог. правдоподобие 57,22951 Крит. Акаике -106,4590
Крит. Шварца -101,7468 Крит. Хеннана-Куинна -105,2089
Исключая константу, наибольшее р-значение получено для переменной 6 (l_T)
Тест Рамсея (РЕСЕТ) -
Нулевая гипотеза: спецификация адекватна
Тестовая статистика: Ф(2, 18) = 0,208498
р-значение = П(Ф(2, 18) > 0,208498) = 0,813736
Тест на нормальное распределение ошибок -
Нулевая гипотеза: ошибки распределены по нормальному закону
Тестовая статистика: Хи-квадрат(2) = 27,1627
р-значение = 1,26382е-006
Вторая статистическая программа дает сходные результаты.
Цитата:
Сообщение от avz
(Сообщение 511657)
Если логарифмировать, то А=0,92, n=0,50, k=-0,042, m=-0,023
|
Проверил. Excel действительно не так считает свободный член и получается (с помощью пакета анализа -0,083526622), а, соответственно,
экспонента равна 0,92
В Excel имеем следующее
Код:
Коэффициенты
Y-пересечение -0,083526622
lnL -0,022963282
ln(s/d) -0,041804663
ln(Re) 0,504268996
|
