自分用メモも兼ねて初投稿です。
Rを使ったモンテカルロシミュレーションで、ベルヌーイ試行から二項分布と幾何分布を作りました。
自分でシミュレーションを考えたり、パラメーターをいじったりすることで、確率分布の理解がより深まると思います。

追記

コメントでご指摘をいただいたので記事を書き直しました。
ありがとうございました。

参考

こちらのスライドが作るきっかけになりました
指数分布とポアソン分布のいけない関係

ベルヌーイ試行、二項分布、幾何分布の関係

一番わかりやすいベルヌーイ試行の例として、コイン投げを想定しましょう。
表が出たら成功とします。このとき
「成功回数」が従うのが二項分布、
「初めて成功した回数」が従うのが幾何分布になります。

シミュレーションを実行

まずはシミュレーション回数とパラメータを決めます。

N <- 10000 # シミュレーション回数
n <- 10 #　コイン投げの回数
p <- 0.5 # 表の出る確率

シミュレーションを行う関数を作ります。

bernoulli <- function() {
	k.sim <- integer(N)
	j.sim <- integer(N)

	for (i in 1:N) {
	
		k <- 0 # 成功回数
		j = Inf # 初めて成功した回
		counter <- 0 # 試行回数のカウンター		
			
		## n回に満たない,または成功していない場合は処理を繰り返す	
		while (counter < n | j == Inf) {

			counter <- counter + 1

			##　ベルヌーイ試行のシミュレーション
			sim <- sample(0:1, size = 1, replace = TRUE, prob = c(1-p, p))
			
			##　試行回数n回目までは二項分布のkをカウント
			if (counter <= n) {
				k <- k + sim
			}
			
			 ##　成功して、かつ　試行回数が j より小さいとき、j　を更新
			if (counter < j && sim == 1) {
				 j <- counter
			}

		}
	k.sim[i] <- k
	j.sim[i] <- j
	}

	return(list(k = k.sim, j = j.sim))
}

シミュレーションを実行して結果を保存。

sim <- bernoulli()
k.sim <- sim$k
j.sim <- sim$j

二項分布

シミュレーション結果のデータを整理。

### 二項分布

# シミュレーション結果
binom.result <- as.data.frame(table(k.sim)/N)

# シミュレーション結果を入れるデータフレームの作成
binom.sim <- data.frame(k = 0:n, freq = integer(n+1))

# 結果を入れる
for (j in 1:(n+1)) {
	for (i in 1:nrow(binom.result)) {
		if (binom.result[i, 1] == binom.sim$k[j]) {
			binom.sim$freq[j] <- binom.result[i, 2]
		}
	}
}

# 二項分布の理論値
binom.theo <- as.data.frame(dbinom(0:n, size = n, prob = p))
binom.theo$k <- 0:n
colnames(binom.theo) <- c("freq", "k")

## 平均の計算
# 理論値
mean.binom.theo <- n*p
# シミュレーション
mean.binom.sim <- sum(k.sim)/N

## 分散の計算
# 理論値
var.binom.theo <- n*p*(1-p)
# シミュレーション
var.binom.sim <- var(k.sim)

幾何分布

同じ要領で幾何分布もやってみます。

# シミュレーション結果
geo.result <- as.data.frame(table(j.sim)/N)

geo.result$j.sim <- as.numeric(as.character(geo.result$j.sim))

nj <- max(geo.result$j.sim)

geo.sim <- data.frame(freq = rep(0, nj), j = 1:nj)

for (j in 1:nj) {
	for (i in 1:nrow(geo.result)) {
		if (geo.result[i, 1] == geo.sim$j[j]) {
			geo.sim$freq[j] <- geo.result[i, 2]
		}
	}
}
geo.sim <- temp2

# 幾何分布の理論値
geo.theo <- as.data.frame(dgeom(1:nj -1, prob = p)) 
# dgeom(x, p) ではx は失敗回数であることに注意

geo.theo$j <- (1:nj)
colnames(geo.theo) <- c("freq", "j")

## 平均の計算
# 理論値
mean.geo.theo <- 1/p
# シミュレーション
mean.geo.sim <- sum(j.sim)/N

## 分散の計算
# 理論値
var.geo.theo <- (1-p)/p^2
# シミュレーション
var.geo.sim <- var(j.sim)

結果の出力

binom.merge <- merge(binom.sim, binom.theo, by = "k")
binom.plot <- rbind(binom.merge$freq.x, binom.merge$freq.y)
colnames(binom.plot) <- c(0:n)

geo.merge <- merge(geo.sim, geo.theo, by = "j")
geo.plot <- rbind(geo.merge$freq.x, geo.merge$freq.y)
colnames(geo.plot) <- c(1:nj)

par(mfrow = c(2,1)) # 上下に画面を分割

barplot(binom.plot,
	col = c("violetred1", "slateblue4"),
	legend.text = c("シミュレーション", "理論値"),
	args.legend = list(x = "topright"),
	beside = TRUE,
	main = "二項分布",
	xlab = "成功回数", ylab = "確率")


barplot(geo.plot,
	col = c("violetred1", "slateblue4"),
	legend.text = c("シミュレーション", "理論値"),
	args.legend = list(x = "topright"),
	beside = TRUE,
	main = "幾何分布",
	xlab = "初めて成功した回", ylab = "確率")


## 結果の出力
result <- function() {
	cat("二項分布\n シミュレーション\n  平均", mean.binom.sim,
		"\n  分散", var.binom.sim,

		"\n 理論値 \n  平均", mean.binom.theo,
		"\n  分散", var.binom.theo,

		"\n\n幾何分布\n シミュレーション\n  平均", mean.geo.sim,
		"\n  分散", var.geo.sim,

		"\n 理論値\n  平均", mean.geo.theo,
		"\n  分散", var.geo.theo)
}
result()

> result()
二項分布
 シミュレーション
  平均 5.0004 
  分散 2.573257 
 理論値 
  平均 5 
  分散 2.5 

幾何分布
 シミュレーション
  平均 1.9983 
  分散 2.053903 
 理論値
  平均 2 
  分散 2>