Laboratorium numer 1. Numeryczne Algorytmy Inżynierskie, Tomasz Owerko KGIiB

Regresja nieliniowa, aproksymacja krzywych dla potrzeb obliczeń geodezyjnych i
inżynierskich:

Wstęp

The R Project for Statistical Computing –

http://www.r-project.org/

Literatura:
D. M. Bates and D. G. Watts. Nonlinear Regression Analysis and
Its Applications. Wiley, 1988.
S. E. Maxwell and H. D. Delaney (1990), "Designing experiments and
analyzing data: A model comparison perspective".

W „R” Równania które nie mogą być zlinearyzowane lub ich linearyzowana postać nie jest znana mogą
być wpasowywane za pomocą metody nls
Postać wstępna funkcji powinna być znana z teorii np. z fizyki zjawiska (termowizja, konstrukcje:
maszty, kominy, wiadukty, kładki itp.)

Zadanie 1 - wprawka:

Stworzyć ramkę jednolitych danych (o liczebności 24) wraz z odpowiadającą zmienną zależną o
pewnym losowym zaszumieniu. Narysować dane oraz krzywą teoretyczną.

Założenie:
a) Radomizujemy pomiary funkcją rnorm
b) Długość wektora opracowania = 24 (pomiar dyskretny)

len=24
x=runif(len)

// jednolicie rozłożona zmienna

y=x^2+rnorm(len,0,0.04)

// wprowadzmy szum losowy do wartości funkcji

ds=data.frame(x=x,y=y)

// tworzymy ramkę danych

ds

// wyświetlając ramkę

x y
1 0.89630372 0.82353177
2 0.02321050 -0.04986645
3 0.41016706 0.08474166
4 0.33757243 0.10262701
5 0.86798392 0.80236827
6 0.34340661 0.10700315
7 0.08304545 0.04435323
8 0.60756028 0.45617045
9 0.21540084 0.10360606
10 0.23889671 0.04242461
11 0.98617209 0.94627259
12 0.98083996 0.96310611
13 0.92455574 0.80253030
14 0.09107546 -0.00885334
15 0.61064812 0.45171589
16 0.35651122 0.11266692
17 0.96864061 0.88724460
18 0.58030270 0.28182619
19 0.42512862 0.20558007
20 0.70878168 0.48703572
21 0.92171906 0.83941377

Laboratorium numer 1. Numeryczne Algorytmy Inżynierskie, Tomasz Owerko KGIiB

22 0.98873698 1.01740812
23 0.63951224 0.43364385
24 0.69595835 0.51077657

plot(y~x)

//rysujemy "dane"

s=seq(0,1,by=0.01)
s

[1] 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14
[16] 0.15 0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29
[31] 0.30 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.40 0.41 0.42 0.43 0.44
[46] 0.45 0.46 0.47 0.48 0.49 0.50 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59
[61] 0.60 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.70 0.71 0.72 0.73 0.74
[76] 0.75 0.76 0.77 0.78 0.79 0.80 0.81 0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89
[91] 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1.00

lines(s,s^2,lty=2,col="green")

// funkcja "modelowa" [w praktyce nie znana!]

Laboratorium numer 1. Numeryczne Algorytmy Inżynierskie, Tomasz Owerko KGIiB

m=nls(y~I(x^potegi),data=ds,start=list(potegi=1),trace=T) )//

model wykładniczy - estymujemy (tylko) wykładnik - proces iteracyjny!

0.696303 : 1
0.08328449 : 1.647024
0.04389877 : 1.922860
0.04367553 : 1.947776
0.04367551 : 1.947541
0.04367551 : 1.947545

summary(m)

Formula: y ~ I(x^potegi)
Parameters:
Estimate Std. Error t value Pr(>|t|)
potegi 1.94755 0.07317 26.62 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.04358 on 23 degrees of freedom
Number of iterations to convergence: 5
Achieved convergence tolerance: 2.020e-07

lines(s,predict(m,list(x=s),lty=1,col="blue"))

// rysyjemy linię z wygenerowanego

modelu

Laboratorium numer 1. Numeryczne Algorytmy Inżynierskie, Tomasz Owerko KGIiB

Na nasze potrzeby:
Def 1.
RSS: residual sum of squaers - świadczy o wariancji w modelowaniu błędów
Def 2.
TSS Total sum of squares - świadczy o wariancji w modelowanych danych
Def 3.
ESS: Explained sum of squares - parametr stosowany do wyrażenia jak dobrze model (regresji) modeluje
zadane dane

Ogólnie:
TSS=ESS+RSS

Szacujemy jakość wpasowania

RSS.p=sum(residuals(m)^2)

// liczymy RSS

RSS.p

[1] 0.04367551

TSS=sum((y-mean(y))^2)

// liczymy TSS

TSS

Laboratorium numer 1. Numeryczne Algorytmy Inżynierskie, Tomasz Owerko KGIiB

[1] 3.047115

1-(RSS.p/TSS)

// jakość modelu wygenerowanego na podstawie nls

[1] 0.9856666

1-sum(x^2-y)^2/TSS

// jakość znanej funkcji (zależy od szumu)

[1] 0.9996179

Może lepszy model?

kwadracik=function(x,w,b)
+ {x^w+b}
m.2=nls(y~kwadracik(x,w,b), data=ds, start=list(w=1,b=0),trace=T)

0.696303 : 1 0
0.07191987 : 1.58621651 -0.01789465
0.04278197 : 1.85001676 -0.01101228
0.04252747 : 1.86518399 -0.01244997
0.04252742 : 1.86456054 -0.01252547
0.04252742 : 1.86458905 -0.01252242

summary(m.2)

Formula: y ~ kwadracik(x, w, b)

Parameters:
Estimate Std. Error t value Pr(>|t|)
w 1.86459 0.12930 14.420 1.08e-12 ***
b -0.01252 0.01656 -0.756 0.458
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04397 on 22 degrees of freedom

Number of iterations to convergence: 5
Achieved convergence tolerance: 2.163e-06

RSS.pb=sum(residuals(m.2)^2)

[1] 0.04252742

1-(RSS.pb/TSS)

[1] 0.9860434

1-(RSS.p/TSS)

[1] 0.9856666

1-sum(x^2-y)^2/TSS

[1] 0.9996179

Anova

anova(m.2,m)

Analysis of Variance Table
Model 1: y ~ kwadracik(x, w, b)
Model 2: y ~ I(x^potegi)
Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1 22 0.052105
2 23 0.053218 -1 -0.0011123 0.4696

0.5003

Komentarz:
Pr (>F) to jest prawdopodobieństwo że odrzucenie hipotezy zerowej
(bardziej skomplikowany model nie pasuje lepiej niż prostszy)
byłoby błędem.
My wiemy że nie powinno być "b" (parametru przecięcia - bo

Laboratorium numer 1. Numeryczne Algorytmy Inżynierskie, Tomasz Owerko KGIiB

zaburzaliśmy czyste x^2!) - więc liczymy na to że wartość Pr(>F)
będzie wysoka - i tak jest.
Dalmierz:
c(t)=Csin(фt)
m(t)=sin(wt+P)
y=(C+Msin(wt+P))*sin(фt)

Przykłady do wykonania:

1. Modulacja Amplitudowa (dalmierze EDM)
x=seq(0,60,by=0.01)
y=(0.05+1*sin(0.1*x+2))*sin(3*x)

2. Funkcja okresowa z trendem liniowym
sinusik=function(x,a,a2,w,b2){a*x+a2*sin(w*x)+b2}