Портал аспирантов - Уравнение множественной регрессии

Портал аспирантов (http://www.aspirantura.spb.ru/forum/index.php)

- Физико-математические науки (http://www.aspirantura.spb.ru/forum/forumdisplay.php?f=128)

- - Уравнение множественной регрессии (http://www.aspirantura.spb.ru/forum/showthread.php?t=13474)

Linka, да ладно. вряд ли кто-то серьезно о чем-то или о ком-то судит. просто многие сначала любят спросить, а потом искать ответы самостоятельно. я тоже такая. а мистер Хог таких осуждает слегка.

Цитата:

Сообщение от avz (Сообщение 511585)

ТС, хотите результат - прицепите табличку с XYZ, расскажите, о чем данные. Пришлю обратно с готовыми коэффициентами регресси, из файла будет понятно, как они получены. Там делов на минуту.

Дай человеку рыбу -- он будет сыт один день, научи его ловить рыбу -- и он сопьётся.

avz, буду очень обязан если поможете.

Нужно вывести зависимость критерия Нуссельта при течении воды в трубе от Re и геометрических параметров турбулизатора потока типа "витая лента", таких как шаг закрутки ленты например S/d и ее толщины 1 и 2 мм
Уравнение должно быть вида Nu = A•Re^n •(S/d)^k•(толщина)^m

Вот данные: зависимости Nu=f(Re)

Re 26575,0 22145,8 17716,6 14763,9 10334,7 7381,9
Nu при 1 и s/d=1 142,3 121,8 114,8 105,8 88,4 74,4
Nu при 1 и s/d=2 139,0 125,4 110,0 102,1 81,0 73,2
Nu при 2 и s/d=1 140,0 120,0 114,7 104,0 87,0 72,0
Nu при 2 и s/d=2 137,0 123,3 109,0 101,0 80,1 71,0

Код:

> # Нужно вывести зависимость критерия Нуссельта при течении воды в трубе от Re и геометрических параметров турбулизатора потока типа "витая лента", таких как шаг закрутки ленты например S/d и ее толщины 1 и 2 мм

> # Уравнение должно быть вида Nu = A•Re^n •(S/d)^k•(толщина)^m

> # Вот данные: зависимости Nu=f(Re)

> # 

> # Re  26575,0        22145,8        17716,6        14763,9        10334,7        7381,9

> # Nu при 1 и s/d=1        142,3        121,8        114,8        105,8        88,4        74,4

> # Nu при 1 и s/d=2        139,0        125,4        110,0        102,1        81,0        73,2

> # Nu при 2 и s/d=1        140,0        120,0        114,7        104,0        87,0        72,0

> # Nu при 2 и s/d=2        137,0        123,3        109,0        101,0        80,1        71,0

> 

> 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

        Re T Sd    Nu

1  26575.0 1  1 142.3

2  22145.8 1  1 121.8

3  17716.6 1  1 114.8

4  14763.9 1  1 105.8

5  10334.7 1  1  88.4

6   7381.9 1  1  74.4

7  26575.0 1  2 139.0

8  22145.8 1  2 125.4

9  17716.6 1  2 110.0

10 14763.9 1  2 102.1

11 10334.7 1  2  81.0

12  7381.9 1  2  73.2

13 26575.0 2  1 140.0

14 22145.8 2  1 120.0

15 17716.6 2  1 114.7

16 14763.9 2  1 104.0

17 10334.7 2  1  87.0

18  7381.9 2  1  72.0

19 26575.0 2  2 137.0

20 22145.8 2  2 123.3

21 17716.6 2  2 109.0

22 14763.9 2  2 101.0

23 10334.7 2  2  80.1

24  7381.9 2  2  71.0

> summary(lm1<-lm(log(Nu)~log(Re)+log(T)+log(Sd)))



Call:

lm(formula = log(Nu) ~ log(Re) + log(T) + log(Sd))



Residuals:

      Min        1Q    Median        3Q       Max 

-0.050707 -0.004867  0.009426  0.012959  0.023083 



Coefficients:

            Estimate Std. Error t value Pr(>|t|)    

(Intercept) -0.19233    0.10997  -1.749  0.09564 .  

log(Re)      0.50427    0.01140  44.253  < 2e-16 ***

log(T)      -0.02296    0.01438  -1.597  0.12605    

log(Sd)     -0.04180    0.01438  -2.907  0.00873 ** 

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



Residual standard error: 0.02442 on 20 degrees of freedom

Multiple R-squared:  0.9899,        Adjusted R-squared:  0.9884 

F-statistic: 656.5 on 3 and 20 DF,  p-value: < 2.2e-16



> confint(lm1)

                  2.5 %       97.5 %

(Intercept) -0.42172064  0.037066329

log(Re)      0.48049934  0.528038653

log(T)      -0.05296572  0.007039156

log(Sd)     -0.07180710 -0.011802225

> summary(lm2<-lm(log(Nu)~log(Re)+log(Sd)))



Call:

lm(formula = log(Nu) ~ log(Re) + log(Sd))



Residuals:

      Min        1Q    Median        3Q       Max 

-0.057637 -0.005429  0.004838  0.016099  0.031042 



Coefficients:

            Estimate Std. Error t value Pr(>|t|)    

(Intercept) -0.20029    0.11384  -1.759   0.0931 .  

log(Re)      0.50427    0.01181  42.706   <2e-16 ***

log(Sd)     -0.04180    0.01490  -2.805   0.0106 *  

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



Residual standard error: 0.0253 on 21 degrees of freedom

Multiple R-squared:  0.9887,        Adjusted R-squared:  0.9876 

F-statistic: 915.8 on 2 and 21 DF,  p-value: < 2.2e-16



> confint(lm2)

                  2.5 %      97.5 %

(Intercept) -0.43702165  0.03645041

log(Re)      0.47971327  0.52882473

log(Sd)     -0.07279929 -0.01081003

> oldpar<-par(mfrow=c(2,2))

> plot(lm2)

> par(oldpar)

> exp(-0.20029)

[1] 0.8184934

И диагностические графики.

http://www.aspirantura.spb.ru/forum/...pictureid=1672

Сравниваем расчетные и фактические значения

Код:

> Nu2<-0.8184934*Re^0.50427*Sd^-0.04180

> Nu2

 [1] 139.36198 127.12040 113.59153 103.61395  86.55759  73.04935 135.38211 123.49012

 [9] 110.34760 100.65496  84.08569  70.96322 139.36198 127.12040 113.59153 103.61395

[17]  86.55759  73.04935 135.38211 123.49012 110.34760 100.65496  84.08569  70.96322

> Nu2-Nu

 [1] -2.93801638  5.32040349 -1.20846824 -2.18605295 -1.84240704 -1.35065290 -3.61789408

 [8] -1.90988108  0.34760260 -1.44504382  3.08569478 -2.23678471 -0.63801638  7.12040349

[15] -1.10846824 -0.38605295 -0.44240704  1.04934710 -1.61789408  0.19011892  1.34760260

[22] -0.34504382  3.98569478 -0.03678471

Примерно так...

Если же делать с помощью нелинейного МНК

Код:

> summary(nlm1<-nls(Nu~A*Re^n*Sd^k*T^m,start=list(A=1,n=1,k=1,m=0)))



Formula: Nu ~ A * Re^n * Sd^k * T^m



Parameters:

  Estimate Std. Error t value Pr(>|t|)    

A  0.79250    0.10084   7.859 1.53e-07 ***

n  0.50806    0.01296  39.198  < 2e-16 ***

k -0.03333    0.01434  -2.325   0.0307 *  

m -0.02086    0.01434  -1.455   0.1611    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



Residual standard error: 2.631 on 20 degrees of freedom



Number of iterations to convergence: 16 

Achieved convergence tolerance: 3.872e-07



> summary(nlm2<-nls(Nu~A*Re^n*Sd^k,start=list(A=1,n=1,k=1)))



Formula: Nu ~ A * Re^n * Sd^k



Parameters:

  Estimate Std. Error t value Pr(>|t|)    

A  0.78664    0.10265   7.664 1.63e-07 ***

n  0.50808    0.01330  38.195  < 2e-16 ***

k -0.03333    0.01471  -2.265   0.0342 *  

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



Residual standard error: 2.7 on 21 degrees of freedom



Number of iterations to convergence: 16 

Achieved convergence tolerance: 1.913e-07

Получаем несколько иные коэффициенты в пределах доверительного интервала. В обоих случаях толщина у нас не при делах.

Код:

> (Nu3<-0.78664*Re^0.50808*Sd^-0.03333)

 [1] 139.23948 126.92047 113.31649 103.29129  86.17087  72.62981 136.05956 124.02189

 [9] 110.72860 100.93235  84.20292  70.97111 139.23948 126.92047 113.31649 103.29129

[17]  86.17087  72.62981 136.05956 124.02189 110.72860 100.93235  84.20292  70.97111

> Nu3-Nu

 [1] -3.06051548  5.12046869 -1.48350586 -2.50870593 -2.22912754 -1.77019073 -2.94043534

 [8] -1.37811227  0.72859768 -1.16764914  3.20292157 -2.22889389 -0.76051548  6.92046869

[15] -1.38350586 -0.70870593 -0.82912754  0.62980927 -0.94043534  0.72188773  1.72859768

[22] -0.06764914  4.10292157 -0.02889389



> confint(nlm2)

         2.5%        97.5%

A  0.59839995  1.031156973

n  0.48048555  0.535954713

k -0.06394209 -0.002731841

Относительная ошибка аппроксимации составляет для первого и второго решения соответственно (в процентах):

Код:

> mean(abs(Nu2-Nu)/Nu)*100

[1] 1.825308

> mean(abs(Nu3-Nu)/Nu)*100

[1] 1.877169

Поиск решения в MS Excel дает два варианта.
Если логарифмировать, то А=0,92, n=0,50, k=-0,042, m=-0,023
Если не логарифмировать, то А=0,82, n=0,50, k=-0,029, m=0,017.

В обоих случаях я бы не сказал, что совпадение расчетных величин и исходных имеет место. L^m и (s/d)^k меня, как бывшего ракетчика, очень смущают, особенно L^m. Вы в критериальное (безразмерное!) уравнение толкаете размерную величину в миллиметрах, и что-то хотите получить. Я думаю, этот паровоз не полетит...

Файл с деталями "поиска решения" в аттаче.

avz, m не при делах в модели. Гляньте t-статистику.

Hogfather, супер-прога! спасибо, что показали! а разбираться в ней можно только по r-documentation, или что-то типа самоучителя есть?

Linka, сейчас, добрые люди дали ссылку, как раз начался курс на курсере
https://www.coursera.org/course/econometrics

А вот тут список учебников
http://r-analytics.blogspot.ru/p/blog-page_20.html

А вот даже тема есть на Портале по данному вопросу
http://www.aspirantura.spb.ru/forum/...ad.php?t=10501

Хог, благодарю, а скиньте сам файл, может будет время разберусь в проге.

Вы в критериальное (безразмерное!) уравнение толкаете размерную величину в миллиметрах, и что-то хотите получить. Я думаю, этот паровоз не полетит...

avz, спасибо, вот это очень полезная мысль, про которую я то позабыл, надо будет пересмотреть параметры.

Цитата:

Сообщение от 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

Цитата:

Сообщение от Hogfather (Сообщение 511670)

еще и на русском ))) отвал башки просто)