「R」の編集履歴（バックアップ）一覧はこちら
「R」(2008/10/06 (月) 21:29:59) の最新版変更点
追加された行は緑色になります。
削除された行は赤色になります。
#contents
* Conditional simulation
** Example
bore <- read.table("borehole.txt", header=T)
v <- variogram(log(lambda)~1, ~x+y, bore)
m <- vgm(.003, "Exp", 30)
plot(v, model=m)
sim <- krige(formula=log(lambda)~1, locations=~x+y, model=m, data=bore, newdata=dom, nmax=15, beta=1.16, nsim=4)
levelplot(exp(z)~x+y|name, map.to.lev(sim, z=c(3:6)), aspect="iso")
levelplot(z~x+y|name, map.to.lev(sim, z=c(3:6)), aspect="iso")
** 3D SGSIM
library(gstat)
library(lattice)
range1d <- seq(from=1, to=100, length=50)
grid3d <- expand.grid(x=range1d, y=range1d, z=range1d)
coordinates(grid3d) <- ~x+y+z
bore3d <- expand.grid(x=c(1,100), y=50, z=range1d)
bore3d <- cbind(bore3d, lambda=3.2)
coordinates(bore3d) <- ~x+y+z
res3D <- krige(formula = log(lambda) ~ 1, bore3d, grid3d, model = vgm(0.03, "Exp", 30), beta=1.16, nmax=15, nsim=1, debug.level=-1)
levelplot(exp(sim1) ~ x + y | z, as.data.frame(res3D))
summary(exp(as.data.frame(res3D)$sim1))
* Histogram
library(lattice)
histogram(meuse$lead, nint=12)
hist(mhw$grain, breaks=seq(2.6, 5.2, by=0.1), col="lightblue", border="red", main="Grain yield")
h <- hist(mhw$grain, breaks = seq(2.6, 5.2, by = 0.2), plot = F)
plot(h, col = heat.colors(length(h$mids))[length(h$count) -
rank(h$count) + 1], ylim = c(0, max(h$count) + 5),
main = "Frequency histogram, Mercer & Hall grain yield",
sub = "Counts shown above bar, actual values shown with rug plot",
xlab = "Grain yield, lbs. per plot")
text(h$mids, h$count, h$count, pos = 3)
hist(mhw$grain, breaks = 20, col = "lightblue", border = "red", freq = F, main = "Grain yield, lbs per plot")
rug(mhw$grain)
lines(density(mhw$grain))
* Empirical cumulative distribution function
plot(ecdf(mhw$grain), pch = 1, xlab = "Mercer & Hall, Grain yield, lbs. per plot",
ylab = "Cumulative proportion of plots", main = "Empirical CDF",
sub = "Quantiles shown with vertical lines")
* Cumulative distribution function
x<-rnorm(1000)
h<-hist(x)
h$counts<-cumsum(h$counts)
plot.histogram(h)
* Graph
** Line
lines(density(mhw$grain))
** Straight line
abline(v = median(grain), col = "blue")
** Title
title("Straw yield predicted by grain yield")
* Plot
plot(mhw$grain, mhw$straw)
pts <- identify(mhw$grain, mhw$straw)
* Variogram
** 3D
#highlight(R){{
library(gstat)
demo(gstat3D)
# $Id: gstat3D.R,v 1.4 2006-02-10 19:05:02 edzer Exp $
# simple demo of 3D interpolation of 50 points with random normal values,
# randomly located in the unit cube
n <- 50
data3D <- data.frame(x = runif(n), y = runif(n), z = runif(n), v = rnorm(n))
coordinates(data3D) = ~x+y+z
range1D <- seq(from = 0, to = 1, length = 20)
grid3D <- expand.grid(x = range1D, y = range1D, z = range1D)
gridded(grid3D) = ~x+y+z
res3D <- krige(formula = v ~ 1, data3D, grid3D, model = vgm(1, "Exp", .2))
library(lattice)
levelplot(var1.pred ~ x + y | z, as.data.frame(res3D))
}}
** Show variogram
show.vgms()
show.vgms(models = "Exp", sill=0.03, range=30, max=100)
* Grid
xy <- expand.grid(x=1:100, y=1:100)
xy <- expand.grid(x=1:100, y=1:100, z=1:100)
* Utility commands
** Generate random numbers from uniform distribution
runif(12, -1, 1)
** Round numbers to two decimal place
round(sample, 2)
** Version of R
R.version
** Find out if the method is available in any package
help.search("lqs")
** Show current directly
getwd()
** Change current directly
setwd("E:\tmp\user")
** Show list of files on the current directly
list.files()
** remove all variables
rm(list = ls())
** Define as a spatial object
coordinates(sim) <- ~x+y
** Min and max values in data
range(mhw$grain)
which.min(grain)
mhw[which.min(grain), "straw"]
** Subset
subset(mhw, grain < 3)
#contents
* Conditional simulation
** Example
bore <- read.table("borehole.txt", header=T)
v <- variogram(log(lambda)~1, ~x+y, bore)
m <- vgm(.003, "Exp", 30)
plot(v, model=m)
sim <- krige(formula=log(lambda)~1, locations=~x+y, model=m, data=bore, newdata=dom, nmax=15, beta=1.16, nsim=4)
levelplot(exp(z)~x+y|name, map.to.lev(sim, z=c(3:6)), aspect="iso")
levelplot(z~x+y|name, map.to.lev(sim, z=c(3:6)), aspect="iso")
** 3D SGSIM
library(gstat)
library(lattice)
range1d <- seq(from=1, to=100, length=50)
grid3d <- expand.grid(x=range1d, y=range1d, z=range1d)
coordinates(grid3d) <- ~x+y+z
bore3d <- expand.grid(x=c(1,100), y=50, z=range1d)
bore3d <- cbind(bore3d, lambda=3.2)
coordinates(bore3d) <- ~x+y+z
res3D <- krige(formula = log(lambda) ~ 1, bore3d, grid3d, model = vgm(0.03, "Exp", 30), beta=1.16, nmax=15, nsim=1, debug.level=-1)
levelplot(exp(sim1) ~ x + y | z, as.data.frame(res3D))
summary(exp(as.data.frame(res3D)$sim1))
* Histogram
library(lattice)
histogram(meuse$lead, nint=12)
hist(mhw$grain, breaks=seq(2.6, 5.2, by=0.1), col="lightblue", border="red", main="Grain yield")
h <- hist(mhw$grain, breaks = seq(2.6, 5.2, by = 0.2), plot = F)
plot(h, col = heat.colors(length(h$mids))[length(h$count) -
rank(h$count) + 1], ylim = c(0, max(h$count) + 5),
main = "Frequency histogram, Mercer & Hall grain yield",
sub = "Counts shown above bar, actual values shown with rug plot",
xlab = "Grain yield, lbs. per plot")
text(h$mids, h$count, h$count, pos = 3)
hist(mhw$grain, breaks = 20, col = "lightblue", border = "red", freq = F, main = "Grain yield, lbs per plot")
rug(mhw$grain)
lines(density(mhw$grain))
* Empirical cumulative distribution function
plot(ecdf(mhw$grain), pch = 1, xlab = "Mercer & Hall, Grain yield, lbs. per plot",
ylab = "Cumulative proportion of plots", main = "Empirical CDF",
sub = "Quantiles shown with vertical lines")
* Cumulative distribution function
x<-rnorm(1000)
h<-hist(x)
h$counts<-cumsum(h$counts)
plot.histogram(h)
* Graph
** Line
lines(density(mhw$grain))
** Straight line
abline(v = median(grain), col = "blue")
** Title
title("Straw yield predicted by grain yield")
* Plot
plot(mhw$grain, mhw$straw)
pts <- identify(mhw$grain, mhw$straw)
* Variogram
** 3D
#highlight(R){{
library(gstat)
demo(gstat3D)
# $Id: gstat3D.R,v 1.4 2006-02-10 19:05:02 edzer Exp $
# simple demo of 3D interpolation of 50 points with random normal values,
# randomly located in the unit cube
n <- 50
data3D <- data.frame(x = runif(n), y = runif(n), z = runif(n), v = rnorm(n))
coordinates(data3D) = ~x+y+z
range1D <- seq(from = 0, to = 1, length = 20)
grid3D <- expand.grid(x = range1D, y = range1D, z = range1D)
gridded(grid3D) = ~x+y+z
res3D <- krige(formula = v ~ 1, data3D, grid3D, model = vgm(1, "Exp", .2))
library(lattice)
levelplot(var1.pred ~ x + y | z, as.data.frame(res3D))
}}
** Show variogram
show.vgms()
show.vgms(models = "Exp", sill=0.03, range=30, max=100)
* Grid
xy <- expand.grid(x=1:100, y=1:100)
xy <- expand.grid(x=1:100, y=1:100, z=1:100)
* Utility commands
** Interpolation
n <- seq(0.001, 0.2, by=0.001)
k <- (234.0*(n)^1.25+2094.0*(10*n)^3.88)*1e-18
f_kn <- approxfun(k, n)
** Generate random numbers from uniform distribution
runif(12, -1, 1)
** Round numbers to two decimal place
round(sample, 2)
** Version of R
R.version
** Find out if the method is available in any package
help.search("lqs")
** Show current directly
getwd()
** Change current directly
setwd("E:\tmp\user")
** Show list of files on the current directly
list.files()
** remove all variables
rm(list = ls())
** Define as a spatial object
coordinates(sim) <- ~x+y
** Min and max values in data
range(mhw$grain)
which.min(grain)
mhw[which.min(grain), "straw"]
** Subset
subset(mhw, grain < 3)
