#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 = ""))