Rで等比数列(幾何数列)を求めたいと思った。
等比数列を求める既存のRの関数を見つけられなかったので、作ってみた。
でも下記の関数を使わなくてもわりと簡単に求められることも関数作っていてわかった。
{seq.geometric}
等比級数列を求めるRの関数を含む(というかそれ単一の関数しか入っていない)パッケージ{seq.geometric}
を作ってGithubに公開した。
インストールは
install.packages("remotes") remotes::install_github("indenkun/seq.geometric")
でできる。
使い方は、
library(seq.geometric)
でロード。
関数は唯一つ、seq_geometric()
のみ。
from
で初項を指定し、to
で最終項か最終項に最も近い値を指定し、by.rate
で公比を指定することで、"初項*公比n“の等比数列を”from“から”to"に最も近い値まで求めることができる。
seq_geometric(from = 1, to = 128, by.ratio = 2)
## [1] 1 2 4 8 16 32 64 128
求める数列の解が収束しない場合(収束しないパターンの値を指定した場合)はエラーが出る。
seq_geometric(from = 1, to = 100, by.ratio = 0.1)
## Error: The calculation results do not converge.
求める数列の長さをlength.out
で指定することもできる。
seq_geometric(from = 1, to = 128, by.ratio = 2, length.out = 3)
## [1] 1 2 4
{base}
のseq()
に近い挙動をするように調整した。seq()
のalong.with
は挙動がイマイチよくわからなかったので搭載していない。
等比数列を上記の関数を使わずに求める方法
等比数列を求める既存の関数を見つけることができなかった、実は結構簡単に求めることができることがわかった。だから多分関数がないんだと思う。そんなに頻繁に等比数列を求めるわけじゃないだろうし、自分でサクッと計算しろということだと思う。
簡単なバージョンを書く前に、初項に公比をかけて、第二項を求めて、第二項に公比をかけて第三項を求めて……というのをぱっと思いつくのでこれを実現すると次のようになる。
# 1を初項として1000000までの公比10の等比数列を上記の処理で求める場合 # 第二項以降の答えをansに格納することにする ans <- NULL # ループ処理用のカウンタをiとしておく i <- 1 # n-1項をaとしておく a <- 1 while(a < 1000000){ # 定義?通りn-1項に公比をかけて第n項を求める処理を値が1000000になるまでループする。 a <- a * 10 ans[i] <- a i <- i + 1 } # 初項とそれ以降の項が格納されているansをくっつけると求める数列がでる c(1, ans)
## [1] 1e+00 1e+01 1e+02 1e+03 1e+04 1e+05 1e+06
いろいろな条件分岐のパターンがありえるので、実際にこの処理で関数化する場合にはもう少しif文やブレイクポイントが増えるが基本はこんな感じ。
ループ処理をしているのでちょっぴり遅そうだが、実際のところはパケージに採用した処理と比較も計算内容が単純なのか微妙にしか劣らないが、ちょっと長い(がこれしか思いつかないときはまぁいいかなぁと思っていた)。
ベクトルで処理すると、簡単だということがわかった。つまり簡単なバージョンとしては次のようになる。
# 1を初項として1000000までの公比10の等比数列を上記の処理で求める場合 # 求める数列の長さを計算する # 最終項/初項が公比の何乗になるのかを求める # これが求める数列の長さ-1(初項が公比の0乗となるので)になる length.out <- log(1000000 / 1, base = 10) # 求める数列は初項*公比^(0:length.out)で計算できる。 1 * 10 ^ (0:length.out)
## [1] 1e+00 1e+01 1e+02 1e+03 1e+04 1e+05 1e+06
わずか二行である。
数列の長さがわかっていれば最後の一行だけでいける。しかもベクトルだから、ループ処理するよりもわずかに早い。
上のseq_geometric()
で採用しているのもこの処理で、あとは解が収束しない場合などで処理をストップしたり、任意の長さで数列を切り上げるようにしただけ。