算術平均

　python では numpy.mean を使って配列に含まれる値の算術平均を取得できます。

numpy.mean
https://docs.scipy.org/doc/numpy/reference/generated/numpy.mean.html

>>> array
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
>>> numpy.mean(array)
5.5

　ruby では同様のメソッドがみつからなかったので、配列内の値を合計したものを要素数で割って平均を取得します。

irb(main):020:0* array                              
=> [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
irb(main):021:0> array.inject(:+) / array.size.to_f
=> 5.5

分散

　python で分散を算出するには numpy.var を使用します。

numpy.var
https://docs.scipy.org/doc/numpy/reference/generated/numpy.var.html

>>> array
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
>>> numpy.var(array)
8.25

　こちらも ruby ではメソッドがみつからなかったので、下記サイトを参考にさせていただきました。

Rubyで標本の分散・標準偏差・変動係数算出まで | WEBサービス創造記

　平均、偏差平方和を算出し、要素数で割って分散を算出しています。

irb(main):044:0* array
=> [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
irb(main):045:0> mean = array.inject(:+) / array.size.to_f                         
=> 5.5
irb(main):046:0> sum_of_squares = array.inject(0) {|sum, i| sum + (i - mean) ** 2} 
=> 82.5
irb(main):047:0> sum_of_squares / array.size.to_f                                  
=> 8.25

等間隔の数値配列生成

　ある一定範囲内で、指定したステップごとの数値の配列の取得は、 python では numpy.arrange を利用します。

numpy.arange
https://docs.scipy.org/doc/numpy/reference/generated/numpy.arange.html

>>> numpy.arange(-10, 10, 2)                             
array([-10,  -8,  -6,  -4,  -2,   0,   2,   4,   6,   8])

　ruby では範囲オブジェクトと step メソッドを利用して生成します。

irb(main):059:0* s = 2
=> 2
irb(main):060:0> (-10..10-s).step(s).to_a
=> [-10, -8, -6, -4, -2, 0, 2, 4, 6, 8]

確率密度関数

　python では確率密度関数の利用は scipy.stats.norm から行います。

scipy.stats.norm
https://docs.scipy.org/doc/scipy-0.16.1/reference/generated/scipy.stats.norm.html

>>> from scipy.stats import norm
>>> orig = norm(loc=0, scale=1)
>>> orig                                                              
<scipy.stats._distn_infrastructure.rv_frozen object at 0x7f095ae1b950>
>>> linex = numpy.arange(-10, 10, 2)                     
>>> linex                                                
array([-10,  -8,  -6,  -4,  -2,   0,   2,   4,   6,   8])
>>> orig.pdf(linex)                                         
array([  7.69459863e-23,   5.05227108e-15,   6.07588285e-09,
         1.33830226e-04,   5.39909665e-02,   3.98942280e-01,
         5.39909665e-02,   1.33830226e-04,   6.07588285e-09,
         5.05227108e-15])

　ruby では rubystats という gem があったので、そちらを使ってみます。

github.com

Distribution::Normal というのもあるのですが、こちらは平均と標準偏差がそれぞれ 0, 1 のケースにしか適用できなかったので、 rubystats を使うことにしました。

Distribution http://www.rubydoc.info/github/sciruby/distribution

irb(main):002:0* require 'rubystats/normal_distribution'
=> true
irb(main):003:0> norm = Rubystats::NormalDistribution.new(0, 1)
=> #<Rubystats::NormalDistribution:0x007f6dfa935930 @mean=0, @stdev=1, @variance=1, @pdf_denominator=2.5066282746310007, @cdf_denominator=1.4142135623730951, @use_last=nil>
irb(main):004:0> s = 2
=> 2
irb(main):005:0> linex = (-10..10-s).step(s).to_a
=> [-10, -8, -6, -4, -2, 0, 2, 4, 6, 8]
irb(main):006:0> linex.map {|x| norm.pdf(x) }
=> [7.694598626706419e-23, 5.052271083536892e-15, 6.075882849823285e-09, 0.00013383022576488534, 0.05399096651318805, 0.39894228040143265, 0.05399096651318805, 0.00013383022576488534, 6.075882849823285e-09, 5.052271083536892e-15]

スクリプト全体

　ここまでの内容を踏まえて、スクリプト全体としては下記のように実装しました。

require 'nyaplot'
require 'rubystats/normal_distribution'

def normal_rand(mu = 0, sigma = 1.0)
  random = Random.new
  (Math.sqrt(-2 * Math.log(random.rand)) * Math.sin(2 * Math::PI * random.rand) * sigma) + mu
end

fig = Nyaplot::Frame.new

[2,4,10,100].each do |datapoints|
  ds = datapoints.times.map { normal_rand }

  mu = ds.inject(:+) / ds.size.to_f
  sum_of_squares = ds.inject(0) {|sum, i| sum + (i - mu) ** 2} 
  var = sum_of_squares / ds.size.to_f
  sigma = Math.sqrt(var)

  plot = Nyaplot::Plot.new

  s = 0.1
  linex = (-10..10-s).step(s).to_a
  
  # 真の曲線を表示
  orig = Rubystats::NormalDistribution.new(0, 1)  
  orig_pdfs = linex.map {|x| orig.pdf(x) }
  line = plot.add(:line, linex, orig_pdfs)
  line.color('green')
  line.title('Original')

  # 推定した曲線を表示
  est = Rubystats::NormalDistribution.new(mu, Math.sqrt(sigma))
  est_pdfs = linex.map {|x| est.pdf(x)}
  line = plot.add(:line, linex, est_pdfs)
  line.color('red')
  line.title("Sigma=#{sprintf("%.2f", sigma)}")
  
  # サンプルの表示
  df = Nyaplot::DataFrame.new({x: ds,y: ds.map {|x| orig.pdf(x)}})
  scatter = plot.add_with_df(df, :scatter, :x, :y)
  scatter.color('blue')
  scatter.title('Sample')

  fig.add(plot)
  
  plot.configure do
    x_label("N=#{datapoints}")
    y_label('')
    xrange([-4, 4])
    legend(true)
    height(300)
    width(490)
  end
end

fig.show