#Autor: Krzysztof Pytka (K.Pytka@ci.edu.pl) #Laboratorium 1 z RD2 #Data: 16/02/2012 #POMOC help("plot") ?plot help.search('linear model') ??'linear model' RSiteSearch('linear model') #Wektory x<-c(1,2,3,4,5) c(1,2,3,4,5)->x x2<-1:4 x7<-c(x,x2) x7+2 x3<-seq(from=1, to=11, by=2) x7+x3 #bezsens 1:3+x3 #caveat!!! seq(1,100)^2 2^seq(1,20) x3[1] x3[1:5] x3[c(1,5)] #Zadanie: kolejne nieparzyste potêgi 3: #brzydki sposób: x4<-3^seq(1,20, 1) x4[seq(1,20,by=2)] c(seq(1,20,by=2), x4[seq(1,20,by=2)]) #bezsens cbind(seq(1,20,by=2), x4[seq(1,20,by=2)]) rbind(seq(1,20,by=2), x4[seq(1,20,by=2)]) #kolejno¶æ wykonywania dzia³añ 1:7-1 1:(7-1) 1:7*2 z<-c(seq(1:10), seq(20, 0, by=-2)) sort(z) length(z) unique(z) length(unique(z)) sort(unique(z)) a<-c('tekst', 2) a[2]+1 as.numeric(a[2])+1 #CZYNNIKI plec <- c('dz','ch','dz','dz','ch','ch','dz','ch','dz','dz') plec plec <- factor(plec) plec plec2 <- factor(c('dz','dz','dz','dz','dz', 'ch'),levels=c('dz','ch')) plec2 table(plec) table(plec2) length(plec) wiek <- factor(c('>18','<18','<18','<18','<18','>18','>18','>18','>18','>18')) wiek length(wiek) t <- table(plec,wiek) t prop.table(t) margin.table(t,1) margin.table(t,2) prop.table(t,1) prop.table(t,2) #MACIERZE A<-matrix( nrow=10, ncol=8) B<-matrix(seq(1:70), nrow=10, ncol=7) t(B) B[2,3] B[,4] B[1,] B%*%B t(B)%*%B 1:10%*%B t(1:10)%*%B #data frame uczestnicy<-data.frame(name=c("Madzia", "Szyma", "Brzez", "Borys", "Krzysiek", "Grzegorz"), gender=factor(c("f","f", rep("m",4)))) uczestnicy[1,1] uczestnicy$name[1] uczestnicy$gender[uczestnicy$name=="Krzysiek"] uczestnicy$name[uczestnicy$gender=="m"] ############################################# #Funkcje #Zadanie: #Autor: M. Ramsza # Dany jest uk³ad dynamiczny postaci $x_{n+1}=x_n+1$, gdzie $x_0\in \mathbb{R}$ jest dane. # Napisaæ skrypt generj±cy pierwsze $K\in \mathbb{N}_{++}$ wyrazów ci±gu. # Proszê nie u¿ywaæ ¿adnych pêtli. #Rozwi±zanie z pêtlami dyn1<-function(x0, k) { result<-numeric(k+1) result[1]<-x0 for(i in 2:(k+1) ) { result[i]<-result[i-1]+1 } return(result[2:(k+1)]) } #Rozwi±zanie bez pêtli dyn2<-function(x0,k) {return(rep(x0,k)+1:k)} #Funkcja licz±ca pole i obwód ko³a o zadanym promieniu pole<-function(r) { return(pi*r^2) } obwod<-function(r) { return(2*pi*r) } kolo<-function(x, a="pole") { if(a=="pole") { return(pole(x)) } else if(a=="obwod") { return(obwod(x)) } else ramka<-data.frame(argument=c("Pole", "Obwód"), wynik=c(pole(x), obwod(x))) return(ramka) } kolo(2) #Pole ko³a dla pierwszych 10 liczb naturalnych wynik<-numeric(10) for(i in 1:10) { wynik[i]<-kolo(i) } wynik kolo(seq(1:10), a="pole") #GENEROWANIE LICZB PSEUDOLOSOWYCH runif(10) uczestnicy<-data.frame(name=factor(c("Madzia", "Szyma", "Brzez", "Borys", "Krzysiek", "Grzegorz")), gender=factor(c("f","f", rep("m",4)))) losuj<-function(x) { losowanie<-runif(1, min=0, max=length(x$name)) zaokr<-ceiling(losowanie) if(x$gender[zaokr]=="f") { gender_wynik<-c(" dziewczyna.") } else gender_wynik<-c(" ch³opak.") paste("Obecn± rozgrywkê rozpoczyna: ", x$name[zaokr], " i p³eæ tego gracza to:", gender_wynik) } losuj(uczestnicy) #MONTE CARLO #Liczenie pola ko³a pole_MC<-function(r, N) { punkty<-cbind(runif(N, min=-r, max=r),runif(N, min=-r, max=r)) kolo<-0 for(i in 1:N) { if( (punkty[i,1])^2+(punkty[i,2])^2<=r^2 ){kolo<-kolo+1} } return(kolo/N*(2*r)^2) } pole_MC(1,1000) Pole<-data.frame(liczebnosc=seq(25, 10000, by=25), pole=0) par(mfrow=c(1,1)) for(i in seq(25, 10000, by=25) ) { Pole$pole[i/25]<-pole_MC(1,i) } plot(Pole$pole~Pole$liczebnosc, xlab="liczebno¶æ", ylab="Pole", pch=8, col="orange") abline(pole(1), 0, col="red", lwd=2, lty=3) grid() #Porównanie z rozwi±zaniem bez pêtli mateusz <- function(r,N) return(4 * r ^ 2 * mean(runif(N, min=-r, max=r) ^ 2 + runif(N, min=-r, max=r) ^ 2 < (r ^ 2))) cat(format("N", width = 10), format("mateusz", width = 5), format("pole_MC", width = 5)) for (N in 10 ^ (1:6)) cat(format(N, scientific = F, width = 10), format(system.time(mateusz(1, N))[3], nsmal=2, width = 5), format(system.time(pole_MC(1, N))[3], nsmall=2, width = 5),"\n") ########Regresja dane<-read.csv2(file.choose()) dim(dane) dane<-dane[,2:9] attach(dane) plot(q_vodka~p_vodka, ylab="Zrealizowany popyt", xlab="Cena", main="Wykres zale¿nosci popytu od ceny produktu", col="mediumpurple") grid() pair(dane) plot(q_vodka~p_vodka, ylab="Zrealizowany popyt", xlab="Cena", main="Wykres zale¿nosci popytu od ceny produktu", col="mediumpurple") grid() model<-lm(q_vodka~p_vodka, data=dane) lm(q_vodka~1+p_vodka, data=dane) lm(q_vodka~0+p_vodka, data=dane) summary(model) abline(model, col="green", lwd=3, lty=2) model2<-lm(q_vodka~., data=dane) model3<-lm(q_vodka~.-wrzesien-p_Nesquik, data=dane) ############W£ASNO¦CI ESTYMATORA KMNK - PODEJ¦CIE MONTE CARLO #Autorzy: Krzysztof Pytka (K.Pytka@ci.edu.pl) # oraz Mateusz Zawisza (mat.zawisza@gmail.com) #Data: 17/09/2010 ##Instalacja potrzebnych pakietów install.packages("quantreg", repos = getOption("repos")) ############ set.seed(38166, kind = "default", normal.kind = "default") library(MASS) library(quantreg) alpha<-25 beta<-12 slope_mnk<-c() slope_mnk_2<-c() slope_orr<-c() slope_LAD<-c() n<-25 ratio<-5 x<-runif(n,min=0,max=10) X2<-ginv(t(cbind(1,x))%*%cbind(1,x))[2,] XOR<-ginv(t(cbind(1,x)) %*% cbind(1,x) + 1.5*diag(2))[2,] for(i in 1:10000) { eps<-rnorm(n, mean=0, sd=5) y<-alpha+beta*x+eps slope_mnk[i]<-X2%*%(rbind(1,x)%*%y) slope_orr[i]<-XOR%*%(rbind(1,x)%*%y) slope_LAD[i]<-rq(y~1+x)$coefficients[2] } # close(pb) #OBCI¡¯ONO¦Æ: Klasyczna Metoda Najmniejszych Kwadrat?w vs. Regresja Grzbietowa (ORR) par(mfrow=c(2,1)) hist(slope_mnk, col="orange", main=paste("Estymator nieobci±¿ony (KMNK), 'prawdziwa' beta=",beta), xlim=c(min(slope_mnk, slope_orr), max(slope_orr, slope_mnk))) hist(slope_orr, col="orange", main=paste("Estymator obci±¿ony (ORR), 'prawdziwa' beta=",beta), xlim=c(min(slope_mnk, slope_orr), max(slope_orr, slope_mnk))) ######################################### #EFEKTYWNO¦Æ: Klasyczna Metoda Najmniejszych Kwadratów vs. Najmniejszych Odchyleñ Bezwzglêdnych (LAD) par(mfrow=c(2,1)) hist(slope_mnk, xlim=c(min(slope_mnk, slope_LAD), max(slope_mnk, slope_LAD)), freq=TRUE, col="orange", breaks=10, ylim=c(0,4000), main=paste("Estymator nieobci±¿ony i EFEKTYWNY (KMNK), 'prawdziwa' beta=", beta)) hist(slope_LAD, xlim=c(min(slope_mnk, slope_LAD), max(slope_mnk, slope_LAD)), freq=TRUE, col="orange", breaks=30, ylim=c(0,4000), main=paste("Estymator nieobci±¿ony i NIEEFEKTYWNY (LAD), 'prawdziwa' beta=",beta)) ######################################### #ZGODNO¦Æ #################Próba n x<-runif(n,min=0,max=10) X2<-ginv(t(cbind(1,x))%*%cbind(1,x))[2,] for(i in 1:10000) { eps<-rnorm(n, mean=0, sd=5) y<-alpha+beta*x+eps slope_mnk[i]<-X2%*%(rbind(1,x)%*%y) } slope_mnk_10<-c() x10<-runif(10*n,min=0,max=10) X2_10<-ginv(rbind(1,x10)%*%cbind(1,x10))[2,] #############Próba 10*n for(i in 1:10000) { eps<-rnorm(10*n, mean=0, sd=5) y<-alpha+beta*x10+eps slope_mnk_10[i]<-X2_10%*%(rbind(1,x10)%*%y) } #############Próba 100*n slope_mnk_100<-c() x100<-runif(100*n,min=0,max=10) X2_100<-ginv(t(cbind(1,x100))%*%cbind(1,x100))[2,] # pb <- tkProgressBar(title = "Symulacja w toku...", min = 0, # max = n*1000, width = 300) for(i in 1:10000) { eps<-rnorm(n*100, mean=0, sd=5) eps2<-ratio*eps y<-alpha+beta*x100+eps slope_mnk_100[i]<-X2_100%*%(rbind(1,x100)%*%y) # setTkProgressBar(pb, i, label=paste( round(i/(n*10), 0), # "% wykonane")) } # close(pb) par(mfrow=c(3,1)) hist(slope_mnk_100, xlim=c(min(slope_mnk, slope_mnk_10, slope_mnk_100), max(slope_mnk, slope_mnk_10, slope_mnk_100)), freq=TRUE, col="orange", main=paste("Liczebno¶æ próby jednej symulacji n=",100*n, ", 'prawdziwa' beta=", beta)) hist(slope_mnk_10, xlim=c(min(slope_mnk, slope_mnk_10, slope_mnk_100), max(slope_mnk, slope_mnk_10, slope_mnk_100)), freq=TRUE, col="orange", main=paste("Liczebno¶æ próby jednej symulacji n=",10*n, ", 'prawdziwa' beta=",beta)) hist(slope_mnk, xlim=c(min(slope_mnk, slope_mnk_10, slope_mnk_100), max(slope_mnk, slope_mnk_10, slope_mnk_100)), freq=TRUE, col="orange", main=paste("Liczebno¶æ próby jednej symulacji n=",n, ", 'prawdziwa' beta=",beta)) #B£¡D OSZACOWANIA set.seed(38166, kind = "default", normal.kind = "default") library(MASS) library(quantreg) alpha<-25 beta<-8 slope_mnk<-c() slope_mnk_2<-c() slope_orr<-c() slope_LAD<-c() n<-25 ratio<-5 x<-runif(n,min=0,max=10) X2<-ginv(t(cbind(1,x))%*%cbind(1,x))[2,] XOR<-ginv(t(cbind(1,x)) %*% cbind(1,x) + 1.5*diag(2))[2,] slope_mnk<-c() slope_mnk_2<-c() slope_mnk_3<-c() for(i in 1:10000) { eps<-rnorm(n, mean=0, sd=1) eps2<-ratio*eps eps3<-2*ratio*eps y<-alpha+beta*x+eps slope_mnk[i]<-X2%*%(rbind(1,x)%*%y) y2<-alpha+beta*x+eps2 slope_mnk_2[i]<-X2%*%(rbind(1,x)%*%y2) y3<-alpha+beta*x+eps3 slope_mnk_3[i]<-X2%*%(rbind(1,x)%*%y3) } par(mfrow=c(3,1)) hist(slope_mnk, col="orange", xlim=c(min(slope_mnk_3),max(slope_mnk_3)), main=paste("Odch. stand. skl. losowego = 1, beta=",beta)) hist(slope_mnk_2, col="orange", xlim=c(min(slope_mnk_3),max(slope_mnk_3)), main=paste("Odch. stand. skl. losowego = 5, beta=", beta)) hist(slope_mnk_3, col="orange", xlim=c(min(slope_mnk_3),max(slope_mnk_3)), main=paste("Odch. stand. skl. losowego = 10, beta",beta)) #Autorzy: B. Kamiñski i Mateusz Zawisza #Wspó³czynnik determinacji R^2 jest wyznaczany na podstawie próby, # a wiêc ma okre¶lony rozk³ad, który miêdzy innymi # zale¿y od jej liczebno¶ci. Pos³uguj±c siê symulacjami Monte Carlo, # wyznaczyæ warto¶æ oczekiwan± i 90\%-przedzia³ ufno¶ci tego rozk³adu # dla liczebno¶ci prób bêd±cych wielokrotno¶ciami liczby 10, # ale nie wiêkszymi ni¿ 200. # Przyjmij, ¿e zmienn± obja¶niaj±c± jest X o rozk³adzie N(0,1), # a obja¶nian± Y=X+epsilon, gdzie epsilon jest b³êdem losowym # o rozk³adzie N(0,1) niezale¿nym od X. set.seed(1) sizes <- seq(from = 10, to = 200, by = 10) reps <- 2^10 sim.r.squared <- function(n) { x <- rnorm(n) y <- x + rnorm(n) model <- lm(y ~ x) return(summary(model)$r.squared) } r.squared.mean <- numeric(length(sizes)) r.squared.q5 <- numeric(length(sizes)) r.squared.q95 <- numeric(length(sizes)) for (i in 1:length(sizes)) { result <- replicate(reps, sim.r.squared(sizes[i])) r.squared.mean[i] <- mean(result) r.squared.q5[i] <- quantile(result, 0.05) r.squared.q95[i] <- quantile(result, 0.95) } plot(r.squared.mean~sizes, xlab=sizes, ylim=c(min(r.squared.q5), max(r.squared.q95)), col="red", ylab=) grid() lines(r.squared.q5~sizes, lwd=2, col="green") lines(r.squared.q95~sizes, lwd=2, col="green") # ZADANIE # Copyright by Bogumi³ Kamiñski, Mateusz Zawisza # Podczas szacowania regresji liniowej klasyczn± metod± najmniejszych kwadratów # napotykamy na problem niepewno¶ci oszacowañ warto¶ci parametrów regresji. # Przedstaw graficznie wp³yw liczebno¶ci próby na t± niepewno¶æ oszacowañ. # Wynik zapisz do pliku PDF i wy¶wietl go. # Badanie przeprowad? symulacyjnie dla wielko?ci pr?by # nale¿±cych do zbioru {8,16,32,64}. # Za³ó¿ model postaci: Y = beta_0 + beta_1 * X + epsilon, # gdzie epsilon ~ N(0,1), zmienna X jest jednostajnie roz³o¿ona na przedziale [0;2], # a prawdziwe warto¶ci parametrów to (beta_0, beta_1)=(0,1). library(graphics) set.seed(1) file.name <- "testRegresji.pdf" counts <- 2^(3:6) pdf(file.name) par(mfrow = c(2, 2)) for (count in counts) { count<-8 x <- seq(0, 2, length.out = count) y <- x + rnorm(count) model <- lm(y ~ x) xlab.label <- paste("y =", format(model$coeff[1], digits = 4), "+", format(model$coeff[2], digits = 4), "* x + e") plot(x, y, xlab = xlab.label, ylab = "", main = paste("n =",count)) matlines(x, predict(model, interval = "confidence"), type = "l", lty = c(1, 2, 2), col = "red") abline(0, 1 , col = "blue", lwd = 2) } dev.off() shell.exec(paste(getwd(), "/", file.name, sep = ""))