library(maptools)
# downloaded via http://www.globalmap.org/download/data/shp/1_1/shape_jp_1_1.zip
boundp <- readShapePoly("shape_jp_1_1/bnda_1_1.shp")
transl <- readShapeLines("shape_jp_1_1/transl_1_1.shp")
# Plot shapefile
plot(boundp)
plot(transl)
# plot shapefile within specified boundary
plot(boundp, xlim=c(135, 136), ylim=c(34,35))
lines(transl, xlim=c(135, 136), ylim=c(34,35), col='red')
# plot polygons which name equals 'OSAKA'
boundp_hyogo <- shp[(shp$nam == 'HYOGO'),]
x_range <- bbox(shp_tokyo)[1,]
y_range <- bbox(shp_tokyo)[2,]
plot(boundp_hyogo, xlim=x_range, ylim=y_range)
lines(transl, xlim=x_range, ylim=y_range, col='green')
Showing posts with label R. Show all posts
Showing posts with label R. Show all posts
2012/10/26
R tips: read and visualize shapefile with R
Below sample code shows how to read shapefiles with R and plot them.
2012/08/23
R tips: operating probability density function in R
There are several probability density functions (PDFs) defined in R. You can find defined PDFs with "library(help='stats')" in R-console. Now, let's plot Gaussian distribution N(0,1) with *norm function.
# Probability density function > curve(dnorm(x, mean=0, sd=1, log=FALSE), -3, 3, n=1001, col = "blue")
# Cumulative probability density function > curve(pnorm(x, mean=0, sd=1, log=FALSE), -3, 3, n=1001, col = "blue")
# Probabilistic generative model
> rnorm(n=100, mean=0, sd=1)
[1] 0.76291389 0.32249468 1.43329768 0.39455640 -0.96082434 -0.03498441
[7] -0.98776053 -0.69866914 1.31514002 -1.42647303 -0.46658770 -1.32506735
[13] -1.52074716 0.45169290 -1.11993007 -0.27547321 0.86232800 -2.36687435
[19] -1.34228574 -0.72800613 -0.41979553 -2.24104445 1.86935264 -0.88038071
[25] -0.71816560 -0.53788946 0.50974550 -1.23530999 0.77035374 0.06483997
[31] -0.62251296 1.60636131 0.36695288 0.76654543 0.86672801 -0.75831284
[37] -1.48413034 -0.37661399 1.34613939 -1.20912506 -0.22231654 2.66835415
[43] 0.49528498 0.26435606 -0.74862141 -0.41349885 -0.94850251 0.83038086
[49] 0.29994578 -0.15778865 -1.16049616 2.24109521 0.86357121 1.28825845
[55] 0.50549591 0.08668917 -0.34392882 2.08620620 -0.02514741 -0.38974372
[61] -0.07943300 1.01109304 1.26749487 1.65848787 1.40626520 0.13476383
[67] -0.81541581 -0.23397894 -1.54089193 2.69271063 1.02584790 0.34583390
[73] 0.35808556 -2.23486925 -0.97331144 -0.75742258 -1.49065873 -0.24702416
[79] 1.46008300 0.26804260 1.25457656 0.55169794 -0.34509144 -1.24602872
[85] 0.65070530 -0.69068100 0.98802706 -1.76179226 -0.30577751 0.30197949
[91] -0.18236187 -0.86313716 1.36986087 -1.23400546 -0.78780711 -2.14942716
[97] 1.17743473 -0.95787888 -0.27079632 1.79351280
2012/08/21
R tips: first contact with linear pattern recognition in R.
ラベル:
machine learning,
R,
tips
This is just test code for linear pattern recognition. Below code is introduced in
http://www.slideshare.net/sleipnir002/05-12739580.
> index <- sample(nrow(data), nrow(data)*0.2) #sample data
> index
[1] 83 47 146 59 111 94 121 14 15 122 33 22 87 78 66 133 18 79 19
[20] 9 88 32 116 60 7 104 67 20 86 1
> data[index,]
Sepal.Length Sepal.Width Petal.Length Petal.Width Species Species2
83 5.8 2.7 3.9 1.2 versicolor vervir
47 5.1 3.8 1.6 0.2 setosa setosa
146 6.7 3.0 5.2 2.3 virginica vervir
59 6.6 2.9 4.6 1.3 versicolor vervir
111 6.5 3.2 5.1 2.0 virginica vervir
94 5.0 2.3 3.3 1.0 versicolor vervir
121 6.9 3.2 5.7 2.3 virginica vervir
14 4.3 3.0 1.1 0.1 setosa setosa
15 5.8 4.0 1.2 0.2 setosa setosa
122 5.6 2.8 4.9 2.0 virginica vervir
33 5.2 4.1 1.5 0.1 setosa setosa
22 5.1 3.7 1.5 0.4 setosa setosa
87 6.7 3.1 4.7 1.5 versicolor vervir
78 6.7 3.0 5.0 1.7 versicolor vervir
66 6.7 3.1 4.4 1.4 versicolor vervir
133 6.4 2.8 5.6 2.2 virginica vervir
18 5.1 3.5 1.4 0.3 setosa setosa
79 6.0 2.9 4.5 1.5 versicolor vervir
19 5.7 3.8 1.7 0.3 setosa setosa
9 4.4 2.9 1.4 0.2 setosa setosa
88 6.3 2.3 4.4 1.3 versicolor vervir
32 5.4 3.4 1.5 0.4 setosa setosa
116 6.4 3.2 5.3 2.3 virginica vervir
60 5.2 2.7 3.9 1.4 versicolor vervir
7 4.6 3.4 1.4 0.3 setosa setosa
104 6.3 2.9 5.6 1.8 virginica vervir
67 5.6 3.0 4.5 1.5 versicolor vervir
20 5.1 3.8 1.5 0.3 setosa setosa
86 6.0 3.4 4.5 1.6 versicolor vervir
1 5.1 3.5 1.4 0.2 setosa setosa
> res <- glm(Species~.,iris[index,], family=binomial) #linear regression
警告メッセージ：
glm.fit: 数値的に 0 か 1 である確率が生じました
> res$coefficient
(Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
6.587948 -7.617771 -13.925058 36.005598 -31.809827
> logLik(res) #log-scale likelyhood
'log Lik.' -1.785923e-10 (df=5)
> predict(res,iris,type="response") # predictions from fitting functions
1 2 3 4 5 6
2.220446e-16 4.328295e-13 2.220446e-16 3.870520e-11 2.220446e-16 2.220446e-16
7 8 9 10 11 12
2.220446e-16 2.220446e-16 7.856922e-11 9.478634e-11 2.220446e-16 4.737275e-12
13 14 15 16 17 18
2.231699e-11 2.220446e-16 2.220446e-16 2.220446e-16 2.220446e-16 2.220446e-16
19 20 21 22 23 24
2.220446e-16 2.220446e-16 1.795636e-12 2.220446e-16 2.220446e-16 2.220446e-16
25 26 27 28 29 30
2.326153e-07 2.709491e-10 2.220446e-16 2.220446e-16 2.220446e-16 1.643915e-10
31 32 33 34 35 36
3.088882e-10 2.220446e-16 2.220446e-16 2.220446e-16 3.937881e-12 2.220446e-16
37 38 39 40 41 42
2.220446e-16 2.220446e-16 5.330791e-13 2.220446e-16 2.220446e-16 1.769196e-10
43 44 45 46 47 48
2.220446e-16 2.220446e-16 1.556449e-13 2.220446e-16 2.220446e-16 2.626083e-13
49 50 51 52 53 54
2.220446e-16 2.220446e-16 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
55 56 57 58 59 60
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
61 62 63 64 65 66
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 9.999926e-01 1.000000e+00
67 68 69 70 71 72
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
73 74 75 76 77 78
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
79 80 81 82 83 84
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
85 86 87 88 89 90
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
91 92 93 94 95 96
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
97 98 99 100 101 102
1.000000e+00 1.000000e+00 9.974026e-01 1.000000e+00 1.000000e+00 1.000000e+00
103 104 105 106 107 108
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
109 110 111 112 113 114
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
115 116 117 118 119 120
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
121 122 123 124 125 126
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
127 128 129 130 131 132
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
133 134 135 136 137 138
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
139 140 141 142 143 144
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
145 146 147 148 149 150
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
> levels(data$Species2)
[1] "setosa" "vervir"
> yhat <- levels(data$Species2)[(predict(res,iris,type="response")>0.5)+1]
> yhat
[1] "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa"
[9] "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa"
[17] "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa"
[25] "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa"
[33] "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa"
[41] "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa"
[49] "setosa" "setosa" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[57] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[65] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[73] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[81] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[89] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[97] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[105] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[113] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[121] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[129] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[137] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[145] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
> mean(yhat[index] != data$Species2[index]) #training error
[1] 0
> mean(yhat[-index] != data$Species2[-index]) #estimation error
[1] 0
> table(true=data$Species2,prediction=yhat) #show result
prediction
true setosa vervir
setosa 50 0
vervir 0 100
R tips: inserting factor value into data frame
Here is sample R script for inserting factor value into data frame.
> data <- data.frame(iris)
> cind <- length(data) + 1
> data[data$Species == 'setosa', cind] <- c('setosa')
> data[data$Species != 'setosa', cind] <- c('vervir')
> data[,cind] <- as.factor(data[,cind])
> names(data)
[1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" "Species"
[6] "V6"
> names(data)[6]
[1] "V6"
> names(data)[6] <- c('Species2')
2011/06/05
Hello, ess (1st contact)
I tried to use emacs + ess (Emacs Speaks Statistics) to program R code quickly.
INSTALL
ESS COMMAND REFERENCE (more ess command is introduced in RjpWiki.)
INSTALL
- aptitude install ess
- echo "(require 'ess-site)" >> ~/.emacs.d/init.el
- boot emacs
- M-x 2 => divide emacs buffer
- M-x R => boot R-console
- C-c C-j: execute current line
- C-c M-j: execute current line and move to end of R-console
- C-c C-b: execute current buffer
- C-c M-b: execute current buffer and move to end of R-console
- C-c C-r: execute current region
- C-c M-r: execute current region and move to end of R-console
ESS COMMAND REFERENCE (more ess command is introduced in RjpWiki.)
 |
| Command of ess (by RjpWiki) |
Subscribe to:
Posts (Atom)