私は風力タービンからのデータを分析しています。通常、これは私がExcelで行うようなものですが、データの量にはかなりの負荷がかかります。私はこれまでRを使用したことがないので、いくつかのポインターを探しています。
データはWindSpeedとPowerの2つの列で構成されています。これまでのところ、CSVファイルからデータをインポートし、2つを相互に分散プロットしました。
次にやりたいのは、データを範囲に並べ替えることです。たとえば、WindSpeedがxとyの間にあるすべてのデータを検索し、各範囲で生成された電力の平均を求めて、形成された曲線をグラフ化します。
この平均から、平均の2つの標準偏差のいずれかに該当するデータに基づいて平均を再計算します（基本的に外れ値を無視します）。
どんなポインタでも大歓迎です。
興味のある方のために、これに似たグラフを作成しようとしています。これはかなり標準的なタイプのグラフですが、前述したように、データのせん断量には、Excelよりも重いものが必要です。
もう Excel を使用していないので、データの大雑把なビニングを必要としない最新の統計手法と、外れ値を削除するためのアドホックな方法を使用してみませんか? loess によって実装されているように、局所的に滑らかな回帰です。
csgilespie のサンプル データを少し変更すると、次のようになります。
w_sp <- sample(seq(0, 100, 0.01), 1000)
power <- 1/(1+exp(-(w_sp -40)/5)) + rnorm(1000, sd = 0.1)
plot(w_sp, power)
x_grid <- seq(0, 100, length = 100)
lines(x_grid, predict(loess(power ~ w_sp), x_grid), col = "red", lwd = 3)
まず、問題を具体化するためにいくつかのサンプル データを作成します。
w_sp = sample(seq(0, 100, 0.01), 1000)
power = 1/(1+exp(-(rnorm(1000, mean=w_sp, sd=5) -40)/5))
power[0,5)、[5,10) などの値をビンに入れたいとします。
bin_incr = 5
bins = seq(0, 95, bin_incr)
y_mean = sapply(bins, function(x) mean(power[w_sp >= x & w_sp < (x+bin_incr)]))
これで、対象範囲間の平均値が作成されました。中央値が必要な場合は、 に変更meanしてmedianください。あとは、それらをプロットするだけです。
plot(w_sp, power)
points(seq(2.5, 97.5, 5), y_mean, col=3, pch=16)
平均の 2 標準偏差以内に収まるデータに基づいて平均を取得するには、もう少し複雑な関数を作成する必要があります。
noOutliers = function(x, power, w_sp, bin_incr) {
d = power[w_sp >= x & w_sp < (x + bin_incr)]
m_d = mean(d)
d_trim = mean(d[d > (m_d - 2*sd(d)) & (d < m_d + 2*sd(d))])
return(mean(d_trim))
}
y_no_outliers = sapply(bins, noOutliers, power, w_sp, bin_incr)
@hadley のものと動機が似ているこのバージョンを、 package を使用して適応スムーザーを備えた加法的モデルを使用してミックスに投入しますmgcv。
@hadley が使用するダミーデータを最初に使用
w_sp <- sample(seq(0, 100, 0.01), 1000)
power <- 1/(1+exp(-(w_sp -40)/5)) + rnorm(1000, sd = 0.1)
df <- data.frame(power = power, w_sp = w_sp)
を使用して加法的モデルをgam()適合させ、REML による適応スムーザーと滑らかさの選択を使用します。
require(mgcv)
mod <- gam(power ~ s(w_sp, bs = "ad", k = 20), data = df, method = "REML")
summary(mod)
モデルから予測して適合の標準誤差を取得し、後者を使用して約 95% の信頼区間を生成します
x_grid <- with(df, data.frame(w_sp = seq(min(w_sp), max(w_sp), length = 100)))
pred <- predict(mod, x_grid, se.fit = TRUE)
x_grid <- within(x_grid, fit <- pred$fit)
x_grid <- within(x_grid, upr <- fit + 2 * pred$se.fit)
x_grid <- within(x_grid, lwr <- fit - 2 * pred$se.fit)
すべてをプロットし、黄土を比較のために当てはめます
plot(power ~ w_sp, data = df, col = "grey")
lines(fit ~ w_sp, data = x_grid, col = "red", lwd = 3)
## upper and lower confidence intervals ~95%
lines(upr ~ w_sp, data = x_grid, col = "red", lwd = 2, lty = "dashed")
lines(lwr ~ w_sp, data = x_grid, col = "red", lwd = 2, lty = "dashed")
## add loess fit from @hadley's answer
lines(x_grid$w_sp, predict(loess(power ~ w_sp, data = df), x_grid), col = "blue",
lwd = 3)
商用タービンの近似曲線 (ワイブル分析) の例を次に示します。
http://www.inl.gov/wind/software/
Hadley 自身の ggplot2 も試してみることをお勧めします。彼の Web サイトは素晴らしいリソースです: http://had.co.nz/ggplot2/ .
# If you haven't already installed ggplot2:
install.pacakges("ggplot2", dependencies = T)
# Load the ggplot2 package
require(ggplot2)
# csgillespie's example data
w_sp <- sample(seq(0, 100, 0.01), 1000)
power <- 1/(1+exp(-(w_sp -40)/5)) + rnorm(1000, sd = 0.1)
# Bind the two variables into a data frame, which ggplot prefers
wind <- data.frame(w_sp = w_sp, power = power)
# Take a look at how the first few rows look, just for fun
head(wind)
# Create a simple plot
ggplot(data = wind, aes(x = w_sp, y = power)) + geom_point() + geom_smooth()
# Create a slightly more complicated plot as an example of how to fine tune
# plots in ggplot
p1 <- ggplot(data = wind, aes(x = w_sp, y = power))
p2 <- p1 + geom_point(colour = "darkblue", size = 1, shape = "dot")
p3 <- p2 + geom_smooth(method = "loess", se = TRUE, colour = "purple")
p3 + scale_x_continuous(name = "mph") +
scale_y_continuous(name = "power") +
opts(title = "Wind speed and power")