r_lab_rozw2

2812 days ago by aniag

%r rucircle <- function(n, r=1.0){ estymN<-ceiling(n * 4/pi) #estymacja liczby potrzebnych losowan (2r)^2/(pi*r^2) = 4/pi randData<-runif(2*estymN, max=2*r) #Losujemy 2*estymN jednostajnych punktow z przedzialu (0,2r) points<-matrix(randData,ncol=2) #Ukladamy w pary - jako wspolrzedne kwadratu (0,2r)x(0,2r) colnames(points)<-c("x","y") #Kosmetyka - zmieniamy nazwy kolumn points<-points - r #Przesuwamy i skalujemy kwadrat, by jego srodek byl w p-cie 0.0 points2<-points*points #To nie jest mnozenie macierzy tylko wyraz po wyrazie #Mozna tez points2<-points^2 r2s<-points2[,1]+points2[,2] good<- r2s <= r^2 okN<-sum(good) if (okN<n) return( rbind(points[good,], rucircle(n-okN)) ) #Losujemy brakujace jesli potrzeba else return( points[good,][1:n,] ) #Biezemy tylko pierwsze n dobrych elementow } 
       
%r #Testujemy funkcje v<-rucircle(100000) print("Dim(wynik):") dim(v) print("Summary:") summary(v) plot(v[,"x"],v[,"y"],pch=".") #Dzieki "kosmetyce" mozemy odwolac sie do kolumn przez nazwe png_output=dev.off() 
       
[1] "Dim(wynik):"
[1] 100000      2
[1] "Summary:"
       x                   y             
 Min.   :-0.999660   Min.   :-0.9996166  
 1st Qu.:-0.400968   1st Qu.:-0.4047041  
 Median : 0.002157   Median : 0.0002983  
 Mean   : 0.001124   Mean   :-0.0010346  
 3rd Qu.: 0.406001   3rd Qu.: 0.4012702  
 Max.   : 0.999088   Max.   : 0.9992486  
%r v2<-v*v r2<-v2[,1]+v2[,2] png(width=1000,height=500) par(mfcol=c(1,2)) #Ta instrukcja tworzy "grid" 1x2 do rysowania kolejnych wykresow rhist<-hist(sqrt(r2)) r2hist<-hist(r2) png_output=dev.off() 
       
%r #I jeszcze test na jednostajnosc r^2 #uzyjemy testu X^2-Pearsona i wyniku kubełkowania wykonanego przez hist: r2hist$counts 
       
 [1] 4927 4959 4983 4984 4988 5035 4961 4973 5129 4994 4943 4978
5002 5002 5024 5062 4983 5001 5074
[20] 4998
%r #http://en.wikipedia.org/wiki/Pearson's_chi-square_test#Discrete_uniform_distribution n<-length(r2hist$counts) binProbabilities <- rep(1/n,n) chisq.test(r2hist$counts, p=binProbabilities) 
       
	Chi-squared test for given probabilities

data:  r2hist$counts 
X-squared = 8.3564, df = 19, p-value = 0.9827
%r #Podobnie przyglądamy się, czy kat w kole rozklada sie jednostajnie. atanHist<-hist(atan2(v[,1],v[,2])) png_output=dev.off() #Tutaj widac problemy numeryczne gdy v[,1]/v[,2] bliskie +/- 0 czyli atan ~= +/- pi #Dla ciekawych: https://stat.ethz.ch/pipermail/r-help/2009-January/183776.html #Nalezaloby to inaczej testowac, moze macie jakies pomysly chisq.test(atanHist$counts) chisq.test(atanHist$counts[-c(1,length(atanHist$counts))]) 
       
	Chi-squared test for given probabilities

data:  atanHist$counts 
X-squared = 7848.767, df = 13, p-value < 2.2e-16


	Chi-squared test for given probabilities

data:  atanHist$counts[-c(1, length(atanHist$counts))] 
X-squared = 10.7471, df = 11, p-value = 0.4647