#1. TYÖARKKI
 #datan syöttö näppäimistöltä
x <- c(1.2,2.3,5.1)
x
#toinen tapa on scan-functio (kokeile help(scan))
x <- scan()
0.1 2.2 4.0


          #tyhjä rivi datan lopetusmerkki
x

 #merkkijonon syöttö
x <- scan(,"")
 a b c
 a c b
 b a c


x
#data <- scan("A:\\Tmallit\\Datat\\Weld.txt",list(x=0,y=0))
# data-tiedosto "Weld.txt" hakemistossa "A:\\Datat\\Weld.txt"
data <- scan("C:\\Kurssit\\TMallit\\Datat\\Weld.txt",list(x=0,y=0))
data

names(data)
x <- data$x; y <- data$y
x; y
 # x=virran voimakkuus (ampeereina) hitsauksessa, y=syntyvän hitsisauman minimiläpimitta.
 # The Welding Institute, Abingdon (Altham)
summary(x); summary(y)
#
hist(x)  #x:n histogramma
X <- matrix(scan("A:\\TMallit\\Datat\\Weld.txt"),ncol=2,byrow=T) # T on "true"
X
XX <- t(X)%*%X   # Matriisien t(X)  ja  X  kertolasku
XX
t(X)
weldtable <- read.table("A:\\TMallit\\Datat\\Weld.txt",header=F) # F on "false"
weldtable
# Seuraavassa tehdään lin. regressianalyysi, jossa x:llä selitetään y:tä
plot(x,y)
plot(x,y,xlab="Virrankulutus (amp)",
ylab="Sauman minimileveys (mm)")
title("Sauman minimileveyden riippuvuus\n virrankulutuksesta")
teeny <- lm(y~x)
teeny   # teeny on S-terminologian (sekä S-Plus ja R) mukaan "objekti"
summary(teeny)
anova(teeny)
names(teeny) 
fv1 <- teeny$fitted.values
par(mfrow=c(2,1)) # plottaus kahdelle riville yhteen sarakkeeseen
plot(x,fv1)
plot(x,y)
abline(teeny,x)
par(mfrow=c(1,1) )# yksi plottaus näytölle
#
# Korrelaatiodiagramma tulostetaan postscript-tiedostoon"foo.eps"
#
postscript("foo.eps")
plot(x,y)
dev.off() # Laite pois päältä
#plot(teeny,ask=T)

Y <- cbind(y,fv1); # "cbind" = "columnbind"
Y
matplot(x,Y,type="pl")  #help("matplot")
res <- teeny$residuals   
 plot(x,res)  # Kvadraattinen trendi
 xx <- x*x
xx
teeny2 <- lm(y ~ x + xx) # Kvadraattinen regressio
teeny2
summary(teeny2)
summary(teeny2,correlation=TRUE)
fv2 <- teeny2$fitted.values
plot(x,y)
lines(x,fv2,lty=2)    # lty=line type, p. 354
#
yks <- rep(1,21)
yks
xx <- x*x
X <- cbind(yks,x,xx) # Mallimatriisi
XX <- t(X)%*%X  # Lähes singulaarinen matriisi
eigen(XX)
eigen(XX, only.values=TRUE)

#
# Ortogonaaliset polynomit
#
tpoly2 <- poly(x, degree=2)
#
tpoly2
#
#########			   
teenyp2 <- lm(y ~ poly(x,2))
teenyp2
summary(teenyp2,correlation=TRUE)
anova(teenyp2)
teenyp1 <- lm(y ~ poly(x,1))
summary(teenyp1)
tp2up <- update(teenyp2,.~.-Run)
teenyp1
summary(tp2up)
anova(tp2up)


#########
# Ortogonaalisista polynomeista
#########
tpoly2
a <- c(8.680952, 8.790475) # Ortogonaalisen polynomin alpha-kertoimet
dim(a) <- c(2,1)
P0 <- rep(1,21)
a1 <- sum(x)/sum(P0^2)
P1 <- 2*(x-a1)
b1 <- sum(P1^2)/sum(P0^2)
a2 <- sum(x*P1^2)/sum(P1^2)
P2 <- 2*(x-a2)*P1-b1
# Normeerataan polynomit
tpoly2 <- cbind(P1,P2)
t(tpoly2)%*%tpoly2
#########
Lasketaan ortog. polyn. mallin kertoimet suoraan matriisilausekkeina
#########
Z <- cbind(P0,P1,P2)
a <- solve(t(Z)%*%Z)%*%t(Z)%*%y
#
# Toinen tapa, lsfit-funktio
#
Z <- cbind(P1,P2)
a <- lsfit(Z,y)




#
q()   # Poistuttaessa "yes" tallentaa työtilan