次のような関数があるとします。
fun <- function(x) 2.3*exp(-x/2)
2 から 20 までの間隔でこの関数の平均値を取得したいと思います。
平均を得るために、最初に頭に浮かぶのはこれです:
mean(fun(2:20))
関数に値を与えて平均を計算するのと同じくらい簡単です。
しかし、これを取得するためのより正確な方法が他にあるのではないかと思います。何か案が?
分析的に、次を使用して区間 [a,b] の関数の平均値を決定できます。
したがって、積分を行った後、関数を 2 点で評価し、平均値を解析的に取得できます。-4.6 * exp(0.5 * x)
あなたの場合、これは の積分と の平均値につながります1/(20-2) * (-4.6 * exp(-0.5 * 20) + 4.6 * exp(-0.5 * 2)) = 0.09400203
。
ここで、間隔に沿ってサンプリングし、そのような平均を計算することに焦点を当てます。
get_sample_mean_from_function = function(func, interval, n = 1000) {
interval_samples = seq(interval[1], interval[2], length = n)
function_values = sapply(interval_samples, func)
return(mean(function_values))
}
fun <- function(x) 2.3*exp(-x/2)
get_sample_mean_from_function(fun, interval = c(2,20))
数n
(取得したサンプル数) を増やすことで、回答の精度を上げることができます。これは、サンプルサイズが大きくなるにつれて平均値がどのように発展するかです。
n_list = c(1,4,10,15,25,50,100,500,1000,10e3,100e3,100e4,100e5)
mean_list = sapply(n_list,
function(x) get_sample_mean_from_function(fun,
interval = c(2,20), n = x))
library(ggplot2)
qplot(n_list, mean_list, geom = "point", log = "x")
収束を得るには、少なくとも 1000 個のサンプルが必要であることに注意してください。この数値解を解析値と比較すると、次のようになります。
mean_list - real_value
[1] 7.521207e-01 1.286106e-01 3.984653e-02 2.494165e-02 1.421951e-02
[6] 6.841070e-03 3.355199e-03 6.607662e-04 3.297467e-04 3.291750e-05
[11] 3.291179e-06 3.291122e-07 3.291116e-08
100e5
サンプルの場合でも、倍精度の浮動小数点精度と比較して、解析解と数値解の差が依然として大きいことがわかります。
非常に高い精度がどうしても必要な場合は、分析ソリューションを試してみます. ただし、実際には、妥当な精度を得るには 5000 サンプルで十分です。